Actual source code: ex92.c
petsc-3.9.4 2018-09-11
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() 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 -test_overlap -test_submat
5: */
6: #include <petscmat.h>
8: int main(int argc,char **args)
9: {
10: Mat A,Atrans,sA,*submatA,*submatsA;
12: PetscMPIInt size,rank;
13: PetscInt bs=1,mbs=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs;
14: PetscScalar *vals,rval,one=1.0;
15: IS *is1,*is2;
16: PetscRandom rand;
17: PetscBool flg,TestOverlap,TestSubMat,TestAllcols,test_sorted=PETSC_FALSE;
18: PetscInt vid = -1;
19: #if defined(PETSC_USE_LOG)
20: PetscLogStage stages[2];
21: #endif
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-mat_mbs",&mbs,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-view_id",&vid,NULL);
32: PetscOptionsHasName(NULL,NULL, "-test_overlap", &TestOverlap);
33: PetscOptionsHasName(NULL,NULL, "-test_submat", &TestSubMat);
34: PetscOptionsHasName(NULL,NULL, "-test_allcols", &TestAllcols);
35: PetscOptionsGetBool(NULL,NULL,"-test_sorted",&test_sorted,NULL);
37: MatCreate(PETSC_COMM_WORLD,&A);
38: MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
39: MatSetType(A,MATBAIJ);
40: MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
41: MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
43: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
44: PetscRandomSetFromOptions(rand);
46: MatGetOwnershipRange(A,&rstart,&rend);
47: MatGetSize(A,&M,&N);
48: Mbs = M/bs;
50: PetscMalloc1(bs,&rows);
51: PetscMalloc1(bs,&cols);
52: PetscMalloc1(bs*bs,&vals);
53: PetscMalloc1(M,&idx);
55: /* Now set blocks of values */
56: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
57: for (i=0; i<Mbs; i++) {
58: cols[0] = i*bs; rows[0] = i*bs;
59: for (j=1; j<bs; j++) {
60: rows[j] = rows[j-1]+1;
61: cols[j] = cols[j-1]+1;
62: }
63: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
64: }
65: /* second, add random blocks */
66: for (i=0; i<20*bs; i++) {
67: PetscRandomGetValue(rand,&rval);
68: cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
69: PetscRandomGetValue(rand,&rval);
70: rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
71: for (j=1; j<bs; j++) {
72: rows[j] = rows[j-1]+1;
73: cols[j] = cols[j-1]+1;
74: }
76: for (j=0; j<bs*bs; j++) {
77: PetscRandomGetValue(rand,&rval);
78: vals[j] = rval;
79: }
80: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
81: }
83: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
86: /* make A a symmetric matrix: A <- A^T + A */
87: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
88: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
89: MatDestroy(&Atrans);
90: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
91: MatEqual(A, Atrans, &flg);
92: if (flg) {
93: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
94: } else SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
95: MatDestroy(&Atrans);
97: /* create a SeqSBAIJ matrix sA (= A) */
98: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
99: if (vid >= 0 && vid < size) {
100: if (!rank) printf("A: \n");
101: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
102: if (!rank) printf("sA: \n");
103: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
104: }
106: /* Test sA==A through MatMult() */
107: MatMultEqual(A,sA,10,&flg);
108: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Error in MatConvert(): A != sA");
110: /* Test MatIncreaseOverlap() */
111: PetscMalloc1(nd,&is1);
112: PetscMalloc1(nd,&is2);
114: for (i=0; i<nd; i++) {
115: if (!TestAllcols) {
116: PetscRandomGetValue(rand,&rval);
117: sz = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
119: for (j=0; j<sz; j++) {
120: PetscRandomGetValue(rand,&rval);
121: idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
122: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
123: }
124: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
125: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
126: if (rank == vid) {
127: PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
128: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
129: }
130: } else { /* Test all rows and colums */
131: sz = M;
132: ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);
133: ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);
135: if (rank == vid) {
136: PetscBool colflag;
137: ISIdentity(is2[i],&colflag);
138: printf("[%d] is2[%d], colflag %d\n",rank,(int)i,(int)colflag);
139: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
140: }
141: }
142: }
144: PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
145: PetscLogStageRegister("MatOv_BAIJ",&stages[1]);
147: /* Test MatIncreaseOverlap */
148: if (TestOverlap) {
149: PetscLogStagePush(stages[0]);
150: MatIncreaseOverlap(sA,nd,is2,ov);
151: PetscLogStagePop();
153: PetscLogStagePush(stages[1]);
154: MatIncreaseOverlap(A,nd,is1,ov);
155: PetscLogStagePop();
157: if (rank == vid) {
158: printf("\n[%d] IS from BAIJ:\n",rank);
159: ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
160: printf("\n[%d] IS from SBAIJ:\n",rank);
161: ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
162: }
164: for (i=0; i<nd; ++i) {
165: ISEqual(is1[i],is2[i],&flg);
166: if (!flg) {
167: if (!rank) {
168: ISSort(is1[i]);
169: /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF); */
170: ISSort(is2[i]);
171: /* ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
172: }
173: SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
174: }
175: }
176: }
178: /* Test MatCreateSubmatrices */
179: if (TestSubMat) {
180: if (test_sorted) {
181: for (i = 0; i < nd; ++i) {
182: ISSort(is1[i]);
183: }
184: }
185: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
186: MatCreateSubMatrices(sA,nd,is1,is1,MAT_INITIAL_MATRIX,&submatsA);
188: MatMultEqual(A,sA,10,&flg);
189: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");
191: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
192: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
193: MatCreateSubMatrices(sA,nd,is1,is1,MAT_REUSE_MATRIX,&submatsA);
194: MatMultEqual(A,sA,10,&flg);
195: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatCreateSubmatrices(): A != sA");
197: MatDestroySubMatrices(nd,&submatA);
198: MatDestroySubMatrices(nd,&submatsA);
199: }
201: /* Free allocated memory */
202: for (i=0; i<nd; ++i) {
203: ISDestroy(&is1[i]);
204: ISDestroy(&is2[i]);
205: }
206: PetscFree(is1);
207: PetscFree(is2);
208: PetscFree(idx);
209: PetscFree(rows);
210: PetscFree(cols);
211: PetscFree(vals);
212: MatDestroy(&A);
213: MatDestroy(&sA);
214: PetscRandomDestroy(&rand);
215: PetscFinalize();
216: return ierr;
217: }
220: /*TEST
222: test:
223: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
224: output_file: output/ex92_1.out
226: test:
227: suffix: 2
228: nsize: {{3 4}}
229: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
230: output_file: output/ex92_1.out
232: test:
233: suffix: 3
234: nsize: {{3 4}}
235: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols
236: output_file: output/ex92_1.out
238: test:
239: suffix: 3_sorted
240: nsize: {{3 4}}
241: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted
242: output_file: output/ex92_1.out
244: test:
245: suffix: 4
246: nsize: {{3 4}}
247: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols
248: output_file: output/ex92_1.out
250: TEST*/