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, (char *)0, 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*/