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*/