Actual source code: ex92.c


  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;
 11:   PetscMPIInt    size,rank;
 12:   PetscInt       bs=1,mbs=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs;
 13:   PetscScalar    *vals,rval,one=1.0;
 14:   IS             *is1,*is2;
 15:   PetscRandom    rand;
 16:   PetscBool      flg,TestOverlap,TestSubMat,TestAllcols,test_sorted=PETSC_FALSE;
 17:   PetscInt       vid = -1;
 18: #if defined(PETSC_USE_LOG)
 19:   PetscLogStage  stages[2];
 20: #endif

 22:   PetscInitialize(&argc,&args,(char*)0,help);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 24:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 26:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-mat_mbs",&mbs,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
 29:   PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
 30:   PetscOptionsGetInt(NULL,NULL,"-view_id",&vid,NULL);
 31:   PetscOptionsHasName(NULL,NULL, "-test_overlap", &TestOverlap);
 32:   PetscOptionsHasName(NULL,NULL, "-test_submat", &TestSubMat);
 33:   PetscOptionsHasName(NULL,NULL, "-test_allcols", &TestAllcols);
 34:   PetscOptionsGetBool(NULL,NULL,"-test_sorted",&test_sorted,NULL);

 36:   MatCreate(PETSC_COMM_WORLD,&A);
 37:   MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
 38:   MatSetType(A,MATBAIJ);
 39:   MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
 40:   MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);

 42:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 43:   PetscRandomSetFromOptions(rand);

 45:   MatGetOwnershipRange(A,&rstart,&rend);
 46:   MatGetSize(A,&M,&N);
 47:   Mbs  = M/bs;

 49:   PetscMalloc1(bs,&rows);
 50:   PetscMalloc1(bs,&cols);
 51:   PetscMalloc1(bs*bs,&vals);
 52:   PetscMalloc1(M,&idx);

 54:   /* Now set blocks of values */
 55:   for (j=0; j<bs*bs; j++) vals[j] = 0.0;
 56:   for (i=0; i<Mbs; i++) {
 57:     cols[0] = i*bs; rows[0] = i*bs;
 58:     for (j=1; j<bs; j++) {
 59:       rows[j] = rows[j-1]+1;
 60:       cols[j] = cols[j-1]+1;
 61:     }
 62:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 63:   }
 64:   /* second, add random blocks */
 65:   for (i=0; i<20*bs; i++) {
 66:     PetscRandomGetValue(rand,&rval);
 67:     cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
 68:     PetscRandomGetValue(rand,&rval);
 69:     rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
 70:     for (j=1; j<bs; j++) {
 71:       rows[j] = rows[j-1]+1;
 72:       cols[j] = cols[j-1]+1;
 73:     }

 75:     for (j=0; j<bs*bs; j++) {
 76:       PetscRandomGetValue(rand,&rval);
 77:       vals[j] = rval;
 78:     }
 79:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 80:   }

 82:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 85:   /* make A a symmetric matrix: A <- A^T + A */
 86:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
 87:   MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
 88:   MatDestroy(&Atrans);
 89:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
 90:   MatEqual(A, Atrans, &flg);
 91:   if (flg) {
 92:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 93:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A+A^T is non-symmetric");
 94:   MatDestroy(&Atrans);

 96:   /* create a SeqSBAIJ matrix sA (= A) */
 97:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
 98:   if (vid >= 0 && vid < size) {
 99:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"A:\n");
100:     MatView(A,PETSC_VIEWER_STDOUT_WORLD);
101:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"sA:\n");
102:     MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
103:   }

105:   /* Test sA==A through MatMult() */
106:   MatMultEqual(A,sA,10,&flg);

109:   /* Test MatIncreaseOverlap() */
110:   PetscMalloc1(nd,&is1);
111:   PetscMalloc1(nd,&is2);

113:   for (i=0; i<nd; i++) {
114:     if (!TestAllcols) {
115:       PetscRandomGetValue(rand,&rval);
116:       sz   = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */

118:       for (j=0; j<sz; j++) {
119:         PetscRandomGetValue(rand,&rval);
120:         idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
121:         for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
122:       }
123:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
124:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
125:       if (rank == vid) {
126:         PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n",rank,i,sz);
127:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
128:       }
129:     } else { /* Test all rows and columns */
130:       sz   = M;
131:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);
132:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);

134:       if (rank == vid) {
135:         PetscBool colflag;
136:         ISIdentity(is2[i],&colflag);
137:         PetscPrintf(PETSC_COMM_SELF,"[%d] is2[%" PetscInt_FMT "], colflag %d\n",rank,i,colflag);
138:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
139:       }
140:     }
141:   }

143:   PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
144:   PetscLogStageRegister("MatOv_BAIJ",&stages[1]);

146:   /* Test MatIncreaseOverlap */
147:   if (TestOverlap) {
148:     PetscLogStagePush(stages[0]);
149:     MatIncreaseOverlap(sA,nd,is2,ov);
150:     PetscLogStagePop();

152:     PetscLogStagePush(stages[1]);
153:     MatIncreaseOverlap(A,nd,is1,ov);
154:     PetscLogStagePop();

156:     if (rank == vid) {
157:       PetscPrintf(PETSC_COMM_SELF,"\n[%d] IS from BAIJ:\n",rank);
158:       ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
159:       PetscPrintf(PETSC_COMM_SELF,"\n[%d] IS from SBAIJ:\n",rank);
160:       ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
161:     }

163:     for (i=0; i<nd; ++i) {
164:       ISEqual(is1[i],is2[i],&flg);
165:       if (!flg) {
166:         if (rank == 0) {
167:           ISSort(is1[i]);
168:           ISSort(is2[i]);
169:         }
170:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%" PetscInt_FMT ", is1 != is2",i);
171:       }
172:     }
173:   }

175:   /* Test MatCreateSubmatrices */
176:   if (TestSubMat) {
177:     if (test_sorted) {
178:       for (i = 0; i < nd; ++i) {
179:         ISSort(is1[i]);
180:       }
181:     }
182:     MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
183:     MatCreateSubMatrices(sA,nd,is1,is1,MAT_INITIAL_MATRIX,&submatsA);

185:     MatMultEqual(A,sA,10,&flg);

188:     /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
189:     MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
190:     MatCreateSubMatrices(sA,nd,is1,is1,MAT_REUSE_MATRIX,&submatsA);
191:     MatMultEqual(A,sA,10,&flg);

194:     MatDestroySubMatrices(nd,&submatA);
195:     MatDestroySubMatrices(nd,&submatsA);
196:   }

198:   /* Free allocated memory */
199:   for (i=0; i<nd; ++i) {
200:     ISDestroy(&is1[i]);
201:     ISDestroy(&is2[i]);
202:   }
203:   PetscFree(is1);
204:   PetscFree(is2);
205:   PetscFree(idx);
206:   PetscFree(rows);
207:   PetscFree(cols);
208:   PetscFree(vals);
209:   MatDestroy(&A);
210:   MatDestroy(&sA);
211:   PetscRandomDestroy(&rand);
212:   PetscFinalize();
213:   return 0;
214: }

216: /*TEST

218:    test:
219:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
220:       output_file: output/ex92_1.out

222:    test:
223:       suffix: 2
224:       nsize: {{3 4}}
225:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
226:       output_file: output/ex92_1.out

228:    test:
229:       suffix: 3
230:       nsize: {{3 4}}
231:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols
232:       output_file: output/ex92_1.out

234:    test:
235:       suffix: 3_sorted
236:       nsize: {{3 4}}
237:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted
238:       output_file: output/ex92_1.out

240:    test:
241:       suffix: 4
242:       nsize: {{3 4}}
243:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols
244:       output_file: output/ex92_1.out

246: TEST*/