Actual source code: ex92.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatSBAIJ format.\n";
3: /* Example of usage:
4: mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0
5: */
6: #include <petscmat.h>
10: int main(int argc,char **args)
11: {
12: Mat A,Atrans,sA,*submatA,*submatsA;
14: PetscMPIInt size,rank;
15: PetscInt bs=1,mbs=1,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,M,N,Mbs;
16: PetscScalar *vals,rval,one=1.0;
17: IS *is1,*is2;
18: PetscRandom rand;
19: PetscBool flg;
20: PETSC_UNUSED PetscLogStage stages[2];
21: PetscInt vid = -1;
23: PetscInitialize(&argc,&args,(char *)0,help);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
28: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&mbs,PETSC_NULL);
29: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
30: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
31: PetscOptionsGetInt(PETSC_NULL,"-view_id",&vid,PETSC_NULL);
32:
33: MatCreate(PETSC_COMM_WORLD,&A);
34: MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
35: MatSetType(A,MATBAIJ);
36: MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL);
37: MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL);
39: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
40: PetscRandomSetFromOptions(rand);
42: MatGetOwnershipRange(A,&rstart,&rend);
43: MatGetSize(A,&M,&N);
44: Mbs = M/bs;
46: PetscMalloc(bs*sizeof(PetscInt),&rows);
47: PetscMalloc(bs*sizeof(PetscInt),&cols);
48: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
49: PetscMalloc(M*sizeof(PetscScalar),&idx);
51: /* Now set blocks of values */
52: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
53: for (i=0; i<Mbs; i++){
54: cols[0] = i*bs; rows[0] = i*bs;
55: for (j=1; j<bs; j++) {
56: rows[j] = rows[j-1]+1;
57: cols[j] = cols[j-1]+1;
58: }
59: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
60: }
61: /* second, add random blocks */
62: for (i=0; i<20*bs; i++) {
63: PetscRandomGetValue(rand,&rval);
64: cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
65: PetscRandomGetValue(rand,&rval);
66: rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
67: for (j=1; j<bs; j++) {
68: rows[j] = rows[j-1]+1;
69: cols[j] = cols[j-1]+1;
70: }
72: for (j=0; j<bs*bs; j++) {
73: PetscRandomGetValue(rand,&rval);
74: vals[j] = rval;
75: }
76: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
77: }
79: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
81:
82: /* make A a symmetric matrix: A <- A^T + A */
83: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
84: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
85: MatDestroy(&Atrans);
86: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
87: MatEqual(A, Atrans, &flg);
88: if (flg) {
89: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
90: } else {
91: SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
92: }
93: MatDestroy(&Atrans);
95: /* create a SeqSBAIJ matrix sA (= A) */
96: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
97: if (vid >= 0 && vid < size){
98: if (!rank) printf("A: \n");
99: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
100: if (!rank) printf("sA: \n");
101: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
102: }
104: /* Test sA==A through MatMult() */
105: MatMultEqual(A,sA,10,&flg);
106: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG ,"Error in MatConvert(): A != sA");
107:
108: /* Test MatIncreaseOverlap() */
109: PetscMalloc(nd*sizeof(IS **),&is1);
110: PetscMalloc(nd*sizeof(IS **),&is2);
112: for (i=0; i<nd; i++) {
113: PetscRandomGetValue(rand,&rval);
114: sz = (PetscInt)((0.5 + 0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
115: //sz /= (size*nd*10);
116:
117: if (rank == vid){
118: PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
119: }
120:
121: for (j=0; j<sz; j++) {
122: PetscRandomGetValue(rand,&rval);
123: idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
124: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
125: }
127: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
128: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
130: if (rank == vid){
131: printf("[%d] is2[%d]\n",rank,i);
132: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
133: }
134: }
136: PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
137: PetscLogStageRegister("MatOv_BAIJ",&stages[1]);
139: /* Test MatIncreaseOverlap */
140: PetscLogStagePush(stages[0]);
141: MatIncreaseOverlap(sA,nd,is2,ov);
142: PetscLogStagePop();
144: PetscLogStagePush(stages[1]);
145: MatIncreaseOverlap(A,nd,is1,ov);
146: PetscLogStagePop();
148: if (rank == vid){
149: printf("\n[%d] IS from BAIJ:\n",rank);
150: ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
151: printf("\n[%d] IS from SBAIJ:\n",rank);
152: ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
153: }
155: for (i=0; i<nd; ++i) {
156: ISEqual(is1[i],is2[i],&flg);
157: if (!flg ){
158: if (rank == 0){
159: ISSort(is1[i]);
160: ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
161: ISSort(is2[i]);
162: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
163: }
164: SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
165: }
166: }
167:
168: /* Test MatGetSubmatrices */
169: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
170: MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
172: MatMultEqual(A,sA,10,&flg);
173: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");
175: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
176: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
177: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
178: MatMultEqual(A,sA,10,&flg);
179: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");
181: /* Free allocated memory */
182: for (i=0; i<nd; ++i) {
183: ISDestroy(&is1[i]);
184: ISDestroy(&is2[i]);
185: MatDestroy(&submatA[i]);
186: MatDestroy(&submatsA[i]);
187: }
188: PetscFree(submatA);
189: PetscFree(submatsA);
190: PetscFree(is1);
191: PetscFree(is2);
192: PetscFree(idx);
193: PetscFree(rows);
194: PetscFree(cols);
195: PetscFree(vals);
196: MatDestroy(&A);
197: MatDestroy(&sA);
198: PetscRandomDestroy(&rand);
199: PetscFinalize();
200: return 0;
201: }