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