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