xref: /petsc/src/mat/tests/ex55.c (revision ebead697dbf761eb322f829370bbe90b3bd93fa3)
1 
2 static char help[] = "Tests converting a matrix to another format with MatConvert().\n\n";
3 
4 #include <petscmat.h>
5 /* Usage: mpiexec -n <np> ex55 -verbose <0 or 1> */
6 
7 int main(int argc,char **args)
8 {
9   Mat            C,A,B,D;
10   PetscInt       i,j,ntypes,bs,mbs,m,block,d_nz=6, o_nz=3,col[3],row,verbose=0;
11   PetscMPIInt    size,rank;
12   MatType        type[9];
13   char           file[PETSC_MAX_PATH_LEN];
14   PetscViewer    fd;
15   PetscBool      equal,flg_loadmat,flg,issymmetric;
16   PetscScalar    value[3];
17 
18   PetscFunctionBeginUser;
19   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
20   PetscCall(PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL));
21   PetscCall(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg_loadmat));
22   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
23   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
24 
25   PetscCall(PetscOptionsHasName(NULL,NULL,"-testseqaij",&flg));
26   if (flg) {
27     if (size == 1) {
28       type[0] = MATSEQAIJ;
29     } else {
30       type[0] = MATMPIAIJ;
31     }
32   } else {
33     type[0] = MATAIJ;
34   }
35   if (size == 1) {
36     ntypes  = 3;
37     type[1] = MATSEQBAIJ;
38     type[2] = MATSEQSBAIJ;
39   } else {
40     ntypes  = 3;
41     type[1] = MATMPIBAIJ;
42     type[2] = MATMPISBAIJ;
43   }
44 
45   /* input matrix C */
46   if (flg_loadmat) {
47     /* Open binary file. */
48     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd));
49 
50     /* Load a baij matrix, then destroy the viewer. */
51     PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
52     if (size == 1) {
53       PetscCall(MatSetType(C,MATSEQBAIJ));
54     } else {
55       PetscCall(MatSetType(C,MATMPIBAIJ));
56     }
57     PetscCall(MatSetFromOptions(C));
58     PetscCall(MatLoad(C,fd));
59     PetscCall(PetscViewerDestroy(&fd));
60   } else { /* Create a baij mat with bs>1  */
61     bs   = 2; mbs=8;
62     PetscCall(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL));
63     PetscCall(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL));
64     PetscCheck(bs > 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
65     m    = mbs*bs;
66     PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,NULL,o_nz,NULL,&C));
67     for (block=0; block<mbs; block++) {
68       /* diagonal blocks */
69       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
70       for (i=1+block*bs; i<bs-1+block*bs; i++) {
71         col[0] = i-1; col[1] = i; col[2] = i+1;
72         PetscCall(MatSetValues(C,1,&i,3,col,value,INSERT_VALUES));
73       }
74       i       = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
75       value[0]=-1.0; value[1]=4.0;
76       PetscCall(MatSetValues(C,1,&i,2,col,value,INSERT_VALUES));
77 
78       i       = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
79       value[0]=4.0; value[1] = -1.0;
80       PetscCall(MatSetValues(C,1,&i,2,col,value,INSERT_VALUES));
81     }
82     /* off-diagonal blocks */
83     value[0]=-1.0;
84     for (i=0; i<(mbs-1)*bs; i++) {
85       col[0]=i+bs;
86       PetscCall(MatSetValues(C,1,&i,1,col,value,INSERT_VALUES));
87       col[0]=i; row=i+bs;
88       PetscCall(MatSetValues(C,1,&row,1,col,value,INSERT_VALUES));
89     }
90     PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
91     PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
92   }
93   {
94     /* Check the symmetry of C because it will be converted to a sbaij matrix */
95     Mat Ctrans;
96     PetscCall(MatTranspose(C,MAT_INITIAL_MATRIX,&Ctrans));
97     PetscCall(MatEqual(C,Ctrans,&flg));
98 /*
99     {
100       PetscCall(MatAXPY(C,1.,Ctrans,DIFFERENT_NONZERO_PATTERN));
101       flg  = PETSC_TRUE;
102     }
103 */
104     PetscCall(MatSetOption(C,MAT_SYMMETRIC,flg));
105     PetscCall(MatDestroy(&Ctrans));
106   }
107   PetscCall(MatIsSymmetric(C,0.0,&issymmetric));
108   PetscCall(MatViewFromOptions(C,NULL,"-view_mat"));
109 
110   /* convert C to other formats */
111   for (i=0; i<ntypes; i++) {
112     PetscBool ismpisbaij,isseqsbaij;
113 
114     PetscCall(PetscStrcmp(type[i],MATMPISBAIJ,&ismpisbaij));
115     PetscCall(PetscStrcmp(type[i],MATMPISBAIJ,&isseqsbaij));
116     if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
117     PetscCall(MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A));
118     PetscCall(MatMultEqual(A,C,10,&equal));
119     PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from BAIJ to %s",type[i]);
120     for (j=i+1; j<ntypes; j++) {
121       PetscCall(PetscStrcmp(type[j],MATMPISBAIJ,&ismpisbaij));
122       PetscCall(PetscStrcmp(type[j],MATMPISBAIJ,&isseqsbaij));
123       if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
124       if (verbose>0) {
125         PetscCall(PetscPrintf(PETSC_COMM_WORLD," \n[%d] test conversion between %s and %s\n",rank,type[i],type[j]));
126       }
127 
128       if (rank == 0 && verbose) printf("Convert %s A to %s B\n",type[i],type[j]);
129       PetscCall(MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B));
130       /*
131       if (j == 2) {
132         PetscCall(PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]));
133         PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD));
134         PetscCall(PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]));
135         PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD));
136       }
137        */
138       PetscCall(MatMultEqual(A,B,10,&equal));
139       PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[i],type[j]);
140 
141       if (size == 1 || j != 2) { /* Matconvert from mpisbaij mat to other formats are not supported */
142         if (rank == 0 && verbose) printf("Convert %s B to %s D\n",type[j],type[i]);
143         PetscCall(MatConvert(B,type[i],MAT_INITIAL_MATRIX,&D));
144         PetscCall(MatMultEqual(B,D,10,&equal));
145         PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[j],type[i]);
146 
147         PetscCall(MatDestroy(&D));
148       }
149       PetscCall(MatDestroy(&B));
150       PetscCall(MatDestroy(&D));
151     }
152 
153     /* Test in-place convert */
154     if (size == 1) { /* size > 1 is not working yet! */
155       j = (i+1)%ntypes;
156       PetscCall(PetscStrcmp(type[j],MATMPISBAIJ,&ismpisbaij));
157       PetscCall(PetscStrcmp(type[j],MATMPISBAIJ,&isseqsbaij));
158       if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
159       /* printf("[%d] i: %d, j: %d\n",rank,i,j); */
160       PetscCall(MatConvert(A,type[j],MAT_INPLACE_MATRIX,&A));
161     }
162 
163     PetscCall(MatDestroy(&A));
164   }
165 
166   /* test BAIJ to MATIS */
167   if (size > 1) {
168     MatType ctype;
169 
170     PetscCall(MatGetType(C,&ctype));
171     PetscCall(MatConvert(C,MATIS,MAT_INITIAL_MATRIX,&A));
172     PetscCall(MatMultEqual(A,C,10,&equal));
173     PetscCall(MatViewFromOptions(A,NULL,"-view_conv"));
174     if (!equal) {
175       Mat C2;
176 
177       PetscCall(MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2));
178       PetscCall(MatViewFromOptions(C2,NULL,"-view_conv_assembled"));
179       PetscCall(MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN));
180       PetscCall(MatChop(C2,PETSC_SMALL));
181       PetscCall(MatViewFromOptions(C2,NULL,"-view_err"));
182       PetscCall(MatDestroy(&C2));
183       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
184     }
185     PetscCall(MatConvert(C,MATIS,MAT_REUSE_MATRIX,&A));
186     PetscCall(MatMultEqual(A,C,10,&equal));
187     PetscCall(MatViewFromOptions(A,NULL,"-view_conv"));
188     if (!equal) {
189       Mat C2;
190 
191       PetscCall(MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2));
192       PetscCall(MatViewFromOptions(C2,NULL,"-view_conv_assembled"));
193       PetscCall(MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN));
194       PetscCall(MatChop(C2,PETSC_SMALL));
195       PetscCall(MatViewFromOptions(C2,NULL,"-view_err"));
196       PetscCall(MatDestroy(&C2));
197       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
198     }
199     PetscCall(MatDestroy(&A));
200     PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&A));
201     PetscCall(MatConvert(A,MATIS,MAT_INPLACE_MATRIX,&A));
202     PetscCall(MatViewFromOptions(A,NULL,"-view_conv"));
203     PetscCall(MatMultEqual(A,C,10,&equal));
204     if (!equal) {
205       Mat C2;
206 
207       PetscCall(MatViewFromOptions(A,NULL,"-view_conv"));
208       PetscCall(MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2));
209       PetscCall(MatViewFromOptions(C2,NULL,"-view_conv_assembled"));
210       PetscCall(MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN));
211       PetscCall(MatChop(C2,PETSC_SMALL));
212       PetscCall(MatViewFromOptions(C2,NULL,"-view_err"));
213       PetscCall(MatDestroy(&C2));
214       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
215     }
216     PetscCall(MatDestroy(&A));
217   }
218   PetscCall(MatDestroy(&C));
219 
220   PetscCall(PetscFinalize());
221   return 0;
222 }
223 
224 /*TEST
225 
226    test:
227 
228    test:
229       suffix: 2
230       nsize: 3
231 
232    testset:
233       requires: parmetis
234       output_file: output/ex55_1.out
235       nsize: 3
236       args: -mat_is_disassemble_l2g_type nd -mat_partitioning_type parmetis
237       test:
238         suffix: matis_baij_parmetis_nd
239       test:
240         suffix: matis_aij_parmetis_nd
241         args: -testseqaij
242       test:
243         requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
244         suffix: matis_poisson1_parmetis_nd
245         args: -f ${DATAFILESPATH}/matrices/poisson1
246 
247    testset:
248       requires: ptscotch defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
249       output_file: output/ex55_1.out
250       nsize: 4
251       args: -mat_is_disassemble_l2g_type nd -mat_partitioning_type ptscotch
252       test:
253         suffix: matis_baij_ptscotch_nd
254       test:
255         suffix: matis_aij_ptscotch_nd
256         args: -testseqaij
257       test:
258         requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
259         suffix: matis_poisson1_ptscotch_nd
260         args: -f ${DATAFILESPATH}/matrices/poisson1
261 
262 TEST*/
263