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