xref: /petsc/src/mat/tests/ex55.c (revision bef158480efac06de457f7a665168877ab3c2fd7) !
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   PetscErrorCode ierr;
11   PetscInt       i,j,ntypes,bs,mbs,m,block,d_nz=6, o_nz=3,col[3],row,verbose=0;
12   PetscMPIInt    size,rank;
13   MatType        type[9];
14   char           file[PETSC_MAX_PATH_LEN];
15   PetscViewer    fd;
16   PetscBool      equal,flg_loadmat,flg,issymmetric;
17   PetscScalar    value[3];
18 
19   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20   ierr = PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL);CHKERRQ(ierr);
21   ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg_loadmat);CHKERRQ(ierr);
22   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
23   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
24 
25   ierr = PetscOptionsHasName(NULL,NULL,"-testseqaij",&flg);CHKERRQ(ierr);
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     ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
49 
50     /* Load a baij matrix, then destroy the viewer. */
51     ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
52     if (size == 1) {
53       ierr = MatSetType(C,MATSEQBAIJ);CHKERRQ(ierr);
54     } else {
55       ierr = MatSetType(C,MATMPIBAIJ);CHKERRQ(ierr);
56     }
57     ierr = MatSetFromOptions(C);CHKERRQ(ierr);
58     ierr = MatLoad(C,fd);CHKERRQ(ierr);
59     ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
60   } else { /* Create a baij mat with bs>1  */
61     bs   = 2; mbs=8;
62     ierr = PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr);
63     ierr = PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);CHKERRQ(ierr);
64     if (bs <= 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
65     m    = mbs*bs;
66     ierr = MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,NULL,o_nz,NULL,&C);CHKERRQ(ierr);
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         ierr   = MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
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       ierr    = MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
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       ierr    = MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
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       ierr  = MatSetValues(C,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
87       col[0]=i; row=i+bs;
88       ierr  = MatSetValues(C,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
89     }
90     ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
91     ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
92   }
93   {
94     /* Check the symmetry of C because it will be converted to a sbaij matrix */
95     Mat Ctrans;
96     ierr = MatTranspose(C,MAT_INITIAL_MATRIX,&Ctrans);CHKERRQ(ierr);
97     ierr = MatEqual(C,Ctrans,&flg);CHKERRQ(ierr);
98 /*
99     {
100       ierr = MatAXPY(C,1.,Ctrans,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
101       flg  = PETSC_TRUE;
102     }
103 */
104     ierr = MatSetOption(C,MAT_SYMMETRIC,flg);CHKERRQ(ierr);
105     ierr = MatDestroy(&Ctrans);CHKERRQ(ierr);
106   }
107   ierr = MatIsSymmetric(C,0.0,&issymmetric);CHKERRQ(ierr);
108   ierr = MatViewFromOptions(C,NULL,"-view_mat");CHKERRQ(ierr);
109 
110   /* convert C to other formats */
111   for (i=0; i<ntypes; i++) {
112     PetscBool ismpisbaij,isseqsbaij;
113 
114     ierr = PetscStrcmp(type[i],MATMPISBAIJ,&ismpisbaij);CHKERRQ(ierr);
115     ierr = PetscStrcmp(type[i],MATMPISBAIJ,&isseqsbaij);CHKERRQ(ierr);
116     if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
117     ierr = MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
118     ierr = MatMultEqual(A,C,10,&equal);CHKERRQ(ierr);
119     if (!equal) SETERRQ1(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       ierr = PetscStrcmp(type[j],MATMPISBAIJ,&ismpisbaij);CHKERRQ(ierr);
122       ierr = PetscStrcmp(type[j],MATMPISBAIJ,&isseqsbaij);CHKERRQ(ierr);
123       if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
124       if (verbose>0) {
125         ierr = PetscPrintf(PETSC_COMM_WORLD," \n[%d] test conversion between %s and %s\n",rank,type[i],type[j]);CHKERRQ(ierr);
126       }
127 
128       if (!rank && verbose) printf("Convert %s A to %s B\n",type[i],type[j]);
129       ierr = MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
130       /*
131       if (j == 2) {
132         ierr = PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);CHKERRQ(ierr);
133         ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
134         ierr = PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);CHKERRQ(ierr);
135         ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
136       }
137        */
138       ierr = MatMultEqual(A,B,10,&equal);CHKERRQ(ierr);
139       if (!equal) SETERRQ2(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 && verbose) printf("Convert %s B to %s D\n",type[j],type[i]);
143         ierr = MatConvert(B,type[i],MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr);
144         ierr = MatMultEqual(B,D,10,&equal);CHKERRQ(ierr);
145         if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[j],type[i]);
146 
147         ierr = MatDestroy(&D);CHKERRQ(ierr);
148       }
149       ierr = MatDestroy(&B);CHKERRQ(ierr);
150       ierr = MatDestroy(&D);CHKERRQ(ierr);
151     }
152 
153     /* Test in-place convert */
154     if (size == 1) { /* size > 1 is not working yet! */
155       j = (i+1)%ntypes;
156       ierr = PetscStrcmp(type[j],MATMPISBAIJ,&ismpisbaij);CHKERRQ(ierr);
157       ierr = PetscStrcmp(type[j],MATMPISBAIJ,&isseqsbaij);CHKERRQ(ierr);
158       if (!issymmetric && (ismpisbaij || isseqsbaij)) continue;
159       /* printf("[%d] i: %d, j: %d\n",rank,i,j); */
160       ierr = MatConvert(A,type[j],MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr);
161     }
162 
163     ierr = MatDestroy(&A);CHKERRQ(ierr);
164   }
165 
166   /* test BAIJ to MATIS */
167   if (size > 1) {
168     MatType ctype;
169 
170     ierr = MatGetType(C,&ctype);CHKERRQ(ierr);
171     ierr = MatConvert(C,MATIS,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
172     ierr = MatMultEqual(A,C,10,&equal);CHKERRQ(ierr);
173     ierr = MatViewFromOptions(A,NULL,"-view_conv");CHKERRQ(ierr);
174     if (!equal) {
175       Mat C2;
176 
177       ierr = MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2);CHKERRQ(ierr);
178       ierr = MatViewFromOptions(C2,NULL,"-view_conv_assembled");CHKERRQ(ierr);
179       ierr = MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
180       ierr = MatChop(C2,PETSC_SMALL);CHKERRQ(ierr);
181       ierr = MatViewFromOptions(C2,NULL,"-view_err");CHKERRQ(ierr);
182       ierr = MatDestroy(&C2);CHKERRQ(ierr);
183       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
184     }
185     ierr = MatConvert(C,MATIS,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr);
186     ierr = MatMultEqual(A,C,10,&equal);CHKERRQ(ierr);
187     ierr = MatViewFromOptions(A,NULL,"-view_conv");CHKERRQ(ierr);
188     if (!equal) {
189       Mat C2;
190 
191       ierr = MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2);CHKERRQ(ierr);
192       ierr = MatViewFromOptions(C2,NULL,"-view_conv_assembled");CHKERRQ(ierr);
193       ierr = MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
194       ierr = MatChop(C2,PETSC_SMALL);CHKERRQ(ierr);
195       ierr = MatViewFromOptions(C2,NULL,"-view_err");CHKERRQ(ierr);
196       ierr = MatDestroy(&C2);CHKERRQ(ierr);
197       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
198     }
199     ierr = MatDestroy(&A);CHKERRQ(ierr);
200     ierr = MatDuplicate(C,MAT_COPY_VALUES,&A);CHKERRQ(ierr);
201     ierr = MatConvert(A,MATIS,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr);
202     ierr = MatViewFromOptions(A,NULL,"-view_conv");CHKERRQ(ierr);
203     ierr = MatMultEqual(A,C,10,&equal);CHKERRQ(ierr);
204     if (!equal) {
205       Mat C2;
206 
207       ierr = MatViewFromOptions(A,NULL,"-view_conv");CHKERRQ(ierr);
208       ierr = MatConvert(A,ctype,MAT_INITIAL_MATRIX,&C2);CHKERRQ(ierr);
209       ierr = MatViewFromOptions(C2,NULL,"-view_conv_assembled");CHKERRQ(ierr);
210       ierr = MatAXPY(C2,-1.,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
211       ierr = MatChop(C2,PETSC_SMALL);CHKERRQ(ierr);
212       ierr = MatViewFromOptions(C2,NULL,"-view_err");CHKERRQ(ierr);
213       ierr = MatDestroy(&C2);CHKERRQ(ierr);
214       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in conversion from BAIJ to MATIS");
215     }
216     ierr = MatDestroy(&A);CHKERRQ(ierr);
217   }
218   ierr = MatDestroy(&C);CHKERRQ(ierr);
219 
220   ierr = PetscFinalize();
221   return ierr;
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 !define(PETSC_USE_64BIT_INDICES)
244         suffix: matis_poisson1_parmetis_nd
245         args: -f ${DATAFILESPATH}/matrices/poisson1
246 
247    testset:
248       requires: ptscotch define(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 !define(PETSC_USE_64BIT_INDICES)
259         suffix: matis_poisson1_ptscotch_nd
260         args: -f ${DATAFILESPATH}/matrices/poisson1
261 
262 TEST*/
263