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*/