Actual source code: ex92.c
petsc-3.7.3 2016-08-01
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 -test_overlap -test_submat
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=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs;
16: PetscScalar *vals,rval,one=1.0;
17: IS *is1,*is2;
18: PetscRandom rand;
19: PetscBool flg,TestOverlap,TestSubMat,TestAllcols;
20: PetscInt vid = -1;
21: #if defined(PETSC_USE_LOG)
22: PetscLogStage stages[2];
23: #endif
25: PetscInitialize(&argc,&args,(char*)0,help);
26: MPI_Comm_size(PETSC_COMM_WORLD,&size);
27: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
29: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-mat_mbs",&mbs,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-view_id",&vid,NULL);
34: PetscOptionsHasName(NULL,NULL, "-test_overlap", &TestOverlap);
35: PetscOptionsHasName(NULL,NULL, "-test_submat", &TestSubMat);
36: PetscOptionsHasName(NULL,NULL, "-test_allcols", &TestAllcols);
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
40: MatSetType(A,MATBAIJ);
41: MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
42: MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
44: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
45: PetscRandomSetFromOptions(rand);
47: MatGetOwnershipRange(A,&rstart,&rend);
48: MatGetSize(A,&M,&N);
49: Mbs = M/bs;
51: PetscMalloc1(bs,&rows);
52: PetscMalloc1(bs,&cols);
53: PetscMalloc1(bs*bs,&vals);
54: PetscMalloc1(M,&idx);
56: /* Now set blocks of values */
57: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
58: for (i=0; i<Mbs; i++) {
59: cols[0] = i*bs; rows[0] = i*bs;
60: for (j=1; j<bs; j++) {
61: rows[j] = rows[j-1]+1;
62: cols[j] = cols[j-1]+1;
63: }
64: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
65: }
66: /* second, add random blocks */
67: for (i=0; i<20*bs; i++) {
68: PetscRandomGetValue(rand,&rval);
69: cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
70: PetscRandomGetValue(rand,&rval);
71: rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
72: for (j=1; j<bs; j++) {
73: rows[j] = rows[j-1]+1;
74: cols[j] = cols[j-1]+1;
75: }
77: for (j=0; j<bs*bs; j++) {
78: PetscRandomGetValue(rand,&rval);
79: vals[j] = rval;
80: }
81: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
82: }
84: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
87: /* make A a symmetric matrix: A <- A^T + A */
88: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
89: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
90: MatDestroy(&Atrans);
91: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
92: MatEqual(A, Atrans, &flg);
93: if (flg) {
94: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
95: } else SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
96: MatDestroy(&Atrans);
98: /* create a SeqSBAIJ matrix sA (= A) */
99: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
100: if (vid >= 0 && vid < size) {
101: if (!rank) printf("A: \n");
102: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
103: if (!rank) printf("sA: \n");
104: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
105: }
107: /* Test sA==A through MatMult() */
108: MatMultEqual(A,sA,10,&flg);
109: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Error in MatConvert(): A != sA");
111: /* Test MatIncreaseOverlap() */
112: PetscMalloc1(nd,&is1);
113: PetscMalloc1(nd,&is2);
115: for (i=0; i<nd; i++) {
116: if (!TestAllcols) {
117: PetscRandomGetValue(rand,&rval);
118: sz = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
120: for (j=0; j<sz; j++) {
121: PetscRandomGetValue(rand,&rval);
122: idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
123: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
124: }
125: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
126: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
127: if (rank == vid) {
128: PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
129: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
130: }
131: } else { /* Test all rows and colums */
132: sz = M;
133: ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);
134: ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);
136: if (rank == vid) {
137: PetscBool colflag;
138: ISIdentity(is2[i],&colflag);
139: printf("[%d] is2[%d], colflag %d\n",rank,(int)i,(int)colflag);
140: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
141: }
142: }
143: }
145: PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
146: PetscLogStageRegister("MatOv_BAIJ",&stages[1]);
148: /* Test MatIncreaseOverlap */
149: if (TestOverlap) {
150: PetscLogStagePush(stages[0]);
151: MatIncreaseOverlap(sA,nd,is2,ov);
152: PetscLogStagePop();
154: PetscLogStagePush(stages[1]);
155: MatIncreaseOverlap(A,nd,is1,ov);
156: PetscLogStagePop();
158: if (rank == vid) {
159: printf("\n[%d] IS from BAIJ:\n",rank);
160: ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
161: printf("\n[%d] IS from SBAIJ:\n",rank);
162: ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
163: }
165: for (i=0; i<nd; ++i) {
166: ISEqual(is1[i],is2[i],&flg);
167: if (!flg) {
168: if (!rank) {
169: ISSort(is1[i]);
170: /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF); */
171: ISSort(is2[i]);
172: /* ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
173: }
174: SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
175: }
176: }
177: }
179: /* Test MatGetSubmatrices */
180: if (TestSubMat) {
181: for (i = 0; i < nd; ++i) {
182: ISSort(is1[i]);
183: ISSort(is2[i]);
184: }
185: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
186: MatGetSubMatrices(sA,nd,is2,is2,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 MatGetSubmatrices with MAT_REUSE_MATRIX option */
192: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
193: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
194: MatMultEqual(A,sA,10,&flg);
195: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatGetSubmatrices(): A != sA");
197: for (i=0; i<nd; ++i) {
198: MatDestroy(&submatA[i]);
199: MatDestroy(&submatsA[i]);
200: }
201: PetscFree(submatA);
202: PetscFree(submatsA);
203: }
205: /* Free allocated memory */
206: for (i=0; i<nd; ++i) {
207: ISDestroy(&is1[i]);
208: ISDestroy(&is2[i]);
209: }
210: PetscFree(is1);
211: PetscFree(is2);
212: PetscFree(idx);
213: PetscFree(rows);
214: PetscFree(cols);
215: PetscFree(vals);
216: MatDestroy(&A);
217: MatDestroy(&sA);
218: PetscRandomDestroy(&rand);
219: PetscFinalize();
220: return 0;
221: }