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) {
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));
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"));
121: /* convert C to other formats */
122: for (i = 0; i < ntypes; i++) {
123: PetscBool ismpisbaij, isseqsbaij;
125: PetscCall(PetscStrcmp(type[i], MATMPISBAIJ, &ismpisbaij));
126: PetscCall(PetscStrcmp(type[i], MATMPISBAIJ, &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], MATMPISBAIJ, &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]));
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]);
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]);
156: PetscCall(MatDestroy(&D));
157: }
158: PetscCall(MatDestroy(&B));
159: PetscCall(MatDestroy(&D));
160: }
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], MATMPISBAIJ, &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: }
172: PetscCall(MatDestroy(&A));
173: }
175: /* test BAIJ to MATIS */
176: if (size > 1) {
177: MatType ctype;
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;
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;
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;
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));
229: PetscCall(PetscFinalize());
230: return 0;
231: }
233: /*TEST
235: test:
237: test:
238: suffix: 2
239: nsize: 3
241: testset:
242: requires: parmetis
243: output_file: output/ex55_1.out
244: nsize: 3
245: args: -mat_is_disassemble_l2g_type nd -mat_partitioning_type parmetis
246: test:
247: suffix: matis_baij_parmetis_nd
248: test:
249: suffix: matis_aij_parmetis_nd
250: args: -testseqaij
251: test:
252: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
253: suffix: matis_poisson1_parmetis_nd
254: args: -f ${DATAFILESPATH}/matrices/poisson1
256: testset:
257: requires: ptscotch defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
258: output_file: output/ex55_1.out
259: nsize: 4
260: args: -mat_is_disassemble_l2g_type nd -mat_partitioning_type ptscotch
261: test:
262: suffix: matis_baij_ptscotch_nd
263: test:
264: suffix: matis_aij_ptscotch_nd
265: args: -testseqaij
266: test:
267: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
268: suffix: matis_poisson1_ptscotch_nd
269: args: -f ${DATAFILESPATH}/matrices/poisson1
271: TEST*/