1
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>
7
main(int argc,char ** args)8 int main(int argc,char **args)
9 {
10 Mat A,Atrans,sA,*submatA,*submatsA;
11 PetscErrorCode ierr;
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
22
23 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
25 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
26
27 ierr = PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);CHKERRQ(ierr);
28 ierr = PetscOptionsGetInt(NULL,NULL,"-mat_mbs",&mbs,NULL);CHKERRQ(ierr);
29 ierr = PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);CHKERRQ(ierr);
30 ierr = PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);CHKERRQ(ierr);
31 ierr = PetscOptionsGetInt(NULL,NULL,"-view_id",&vid,NULL);CHKERRQ(ierr);
32 ierr = PetscOptionsHasName(NULL,NULL, "-test_overlap", &TestOverlap);CHKERRQ(ierr);
33 ierr = PetscOptionsHasName(NULL,NULL, "-test_submat", &TestSubMat);CHKERRQ(ierr);
34 ierr = PetscOptionsHasName(NULL,NULL, "-test_allcols", &TestAllcols);CHKERRQ(ierr);
35 ierr = PetscOptionsGetBool(NULL,NULL,"-test_sorted",&test_sorted,NULL);CHKERRQ(ierr);
36
37 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
38 ierr = MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
39 ierr = MatSetType(A,MATBAIJ);CHKERRQ(ierr);
40 ierr = MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);CHKERRQ(ierr);
41 ierr = MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);CHKERRQ(ierr);
42
43 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
44 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
45
46 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
47 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
48 Mbs = M/bs;
49
50 ierr = PetscMalloc1(bs,&rows);CHKERRQ(ierr);
51 ierr = PetscMalloc1(bs,&cols);CHKERRQ(ierr);
52 ierr = PetscMalloc1(bs*bs,&vals);CHKERRQ(ierr);
53 ierr = PetscMalloc1(M,&idx);CHKERRQ(ierr);
54
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 ierr = MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);CHKERRQ(ierr);
64 }
65 /* second, add random blocks */
66 for (i=0; i<20*bs; i++) {
67 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr);
68 cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
69 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr);
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 }
75
76 for (j=0; j<bs*bs; j++) {
77 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr);
78 vals[j] = rval;
79 }
80 ierr = MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);CHKERRQ(ierr);
81 }
82
83 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
84 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
85
86 /* make A a symmetric matrix: A <- A^T + A */
87 ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);CHKERRQ(ierr);
88 ierr = MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
89 ierr = MatDestroy(&Atrans);CHKERRQ(ierr);
90 ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);CHKERRQ(ierr);
91 ierr = MatEqual(A, Atrans, &flg);CHKERRQ(ierr);
92 if (flg) {
93 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
94 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A+A^T is non-symmetric");
95 ierr = MatDestroy(&Atrans);CHKERRQ(ierr);
96
97 /* create a SeqSBAIJ matrix sA (= A) */
98 ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
99 if (vid >= 0 && vid < size) {
100 if (!rank) printf("A: \n");
101 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
102 if (!rank) printf("sA: \n");
103 ierr = MatView(sA,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
104 }
105
106 /* Test sA==A through MatMult() */
107 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
108 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Error in MatConvert(): A != sA");
109
110 /* Test MatIncreaseOverlap() */
111 ierr = PetscMalloc1(nd,&is1);CHKERRQ(ierr);
112 ierr = PetscMalloc1(nd,&is2);CHKERRQ(ierr);
113
114 for (i=0; i<nd; i++) {
115 if (!TestAllcols) {
116 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr);
117 sz = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
118
119 for (j=0; j<sz; j++) {
120 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);CHKERRQ(ierr);
125 ierr = ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);CHKERRQ(ierr);
126 if (rank == vid) {
127 ierr = PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);CHKERRQ(ierr);
128 ierr = ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
129 }
130 } else { /* Test all rows and colums */
131 sz = M;
132 ierr = ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);CHKERRQ(ierr);
133 ierr = ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);CHKERRQ(ierr);
134
135 if (rank == vid) {
136 PetscBool colflag;
137 ierr = ISIdentity(is2[i],&colflag);CHKERRQ(ierr);
138 printf("[%d] is2[%d], colflag %d\n",rank,(int)i,(int)colflag);
139 ierr = ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
140 }
141 }
142 }
143
144 ierr = PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);CHKERRQ(ierr);
145 ierr = PetscLogStageRegister("MatOv_BAIJ",&stages[1]);CHKERRQ(ierr);
146
147 /* Test MatIncreaseOverlap */
148 if (TestOverlap) {
149 ierr = PetscLogStagePush(stages[0]);CHKERRQ(ierr);
150 ierr = MatIncreaseOverlap(sA,nd,is2,ov);CHKERRQ(ierr);
151 ierr = PetscLogStagePop();CHKERRQ(ierr);
152
153 ierr = PetscLogStagePush(stages[1]);CHKERRQ(ierr);
154 ierr = MatIncreaseOverlap(A,nd,is1,ov);CHKERRQ(ierr);
155 ierr = PetscLogStagePop();CHKERRQ(ierr);
156
157 if (rank == vid) {
158 printf("\n[%d] IS from BAIJ:\n",rank);
159 ierr = ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
160 printf("\n[%d] IS from SBAIJ:\n",rank);
161 ierr = ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
162 }
163
164 for (i=0; i<nd; ++i) {
165 ierr = ISEqual(is1[i],is2[i],&flg);CHKERRQ(ierr);
166 if (!flg) {
167 if (!rank) {
168 ierr = ISSort(is1[i]);CHKERRQ(ierr);
169 /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); */
170 ierr = ISSort(is2[i]);CHKERRQ(ierr);
171 /* ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); */
172 }
173 SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
174 }
175 }
176 }
177
178 /* Test MatCreateSubmatrices */
179 if (TestSubMat) {
180 if (test_sorted) {
181 for (i = 0; i < nd; ++i) {
182 ierr = ISSort(is1[i]);CHKERRQ(ierr);
183 }
184 }
185 ierr = MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);CHKERRQ(ierr);
186 ierr = MatCreateSubMatrices(sA,nd,is1,is1,MAT_INITIAL_MATRIX,&submatsA);CHKERRQ(ierr);
187
188 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
189 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");
190
191 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
192 ierr = MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);CHKERRQ(ierr);
193 ierr = MatCreateSubMatrices(sA,nd,is1,is1,MAT_REUSE_MATRIX,&submatsA);CHKERRQ(ierr);
194 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
195 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatCreateSubmatrices(): A != sA");
196
197 ierr = MatDestroySubMatrices(nd,&submatA);CHKERRQ(ierr);
198 ierr = MatDestroySubMatrices(nd,&submatsA);CHKERRQ(ierr);
199 }
200
201 /* Free allocated memory */
202 for (i=0; i<nd; ++i) {
203 ierr = ISDestroy(&is1[i]);CHKERRQ(ierr);
204 ierr = ISDestroy(&is2[i]);CHKERRQ(ierr);
205 }
206 ierr = PetscFree(is1);CHKERRQ(ierr);
207 ierr = PetscFree(is2);CHKERRQ(ierr);
208 ierr = PetscFree(idx);CHKERRQ(ierr);
209 ierr = PetscFree(rows);CHKERRQ(ierr);
210 ierr = PetscFree(cols);CHKERRQ(ierr);
211 ierr = PetscFree(vals);CHKERRQ(ierr);
212 ierr = MatDestroy(&A);CHKERRQ(ierr);
213 ierr = MatDestroy(&sA);CHKERRQ(ierr);
214 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
215 ierr = PetscFinalize();
216 return ierr;
217 }
218
219
220 /*TEST
221
222 test:
223 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
224 output_file: output/ex92_1.out
225
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
231
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
237
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
243
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
249
250 TEST*/
251