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