Actual source code: ex92.c

  1: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel MatSBAIJ format.\n";
  2: /* Example of usage:
  3:       mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat
  4: */
  5: #include <petscmat.h>

  7: int main(int argc, char **args)
  8: {
  9:   Mat           A, Atrans, sA, *submatA, *submatsA;
 10:   PetscMPIInt   size, rank;
 11:   PetscInt      bs = 1, mbs = 10, ov = 1, i, j, k, *rows, *cols, nd = 2, *idx, rstart, rend, sz, M, N, Mbs;
 12:   PetscScalar  *vals, rval, one = 1.0;
 13:   IS           *is1, *is2;
 14:   PetscRandom   rand;
 15:   PetscBool     flg, TestOverlap, TestSubMat, TestAllcols, test_sorted = PETSC_FALSE;
 16:   PetscInt      vid = -1;
 17:   PetscLogStage stages[2];

 19:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 22:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 24:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_mbs", &mbs, NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
 27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
 28:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-view_id", &vid, NULL));
 29:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_overlap", &TestOverlap));
 30:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_submat", &TestSubMat));
 31:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_allcols", &TestAllcols));
 32:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sorted", &test_sorted, NULL));

 34:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 35:   PetscCall(MatSetSizes(A, mbs * bs, mbs * bs, PETSC_DECIDE, PETSC_DECIDE));
 36:   PetscCall(MatSetType(A, MATBAIJ));
 37:   PetscCall(MatSeqBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL));
 38:   PetscCall(MatMPIBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));

 40:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 41:   PetscCall(PetscRandomSetFromOptions(rand));

 43:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 44:   PetscCall(MatGetSize(A, &M, &N));
 45:   Mbs = M / bs;

 47:   PetscCall(PetscMalloc1(bs, &rows));
 48:   PetscCall(PetscMalloc1(bs, &cols));
 49:   PetscCall(PetscMalloc1(bs * bs, &vals));
 50:   PetscCall(PetscMalloc1(M, &idx));

 52:   /* Now set blocks of values */
 53:   for (j = 0; j < bs * bs; j++) vals[j] = 0.0;
 54:   for (i = 0; i < Mbs; i++) {
 55:     cols[0] = i * bs;
 56:     rows[0] = i * bs;
 57:     for (j = 1; j < bs; j++) {
 58:       rows[j] = rows[j - 1] + 1;
 59:       cols[j] = cols[j - 1] + 1;
 60:     }
 61:     PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
 62:   }
 63:   /* second, add random blocks */
 64:   for (i = 0; i < 20 * bs; i++) {
 65:     PetscCall(PetscRandomGetValue(rand, &rval));
 66:     cols[0] = bs * (PetscInt)(PetscRealPart(rval) * Mbs);
 67:     PetscCall(PetscRandomGetValue(rand, &rval));
 68:     rows[0] = rstart + bs * (PetscInt)(PetscRealPart(rval) * mbs);
 69:     for (j = 1; j < bs; j++) {
 70:       rows[j] = rows[j - 1] + 1;
 71:       cols[j] = cols[j - 1] + 1;
 72:     }

 74:     for (j = 0; j < bs * bs; j++) {
 75:       PetscCall(PetscRandomGetValue(rand, &rval));
 76:       vals[j] = rval;
 77:     }
 78:     PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
 79:   }

 81:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 82:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 84:   /* make A a symmetric matrix: A <- A^T + A */
 85:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
 86:   PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN));
 87:   PetscCall(MatDestroy(&Atrans));
 88:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
 89:   PetscCall(MatEqual(A, Atrans, &flg));
 90:   if (flg) {
 91:     PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
 92:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric");
 93:   PetscCall(MatDestroy(&Atrans));

 95:   /* create a SeqSBAIJ matrix sA (= A) */
 96:   PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
 97:   if (vid >= 0 && vid < size) {
 98:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "A:\n"));
 99:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
100:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "sA:\n"));
101:     PetscCall(MatView(sA, PETSC_VIEWER_STDOUT_WORLD));
102:   }

104:   /* Test sA==A through MatMult() */
105:   PetscCall(MatMultEqual(A, sA, 10, &flg));
106:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Error in MatConvert(): A != sA");

108:   /* Test MatIncreaseOverlap() */
109:   PetscCall(PetscMalloc1(nd, &is1));
110:   PetscCall(PetscMalloc1(nd, &is2));

112:   for (i = 0; i < nd; i++) {
113:     if (!TestAllcols) {
114:       PetscCall(PetscRandomGetValue(rand, &rval));
115:       sz = (PetscInt)((0.5 + 0.2 * PetscRealPart(rval)) * mbs); /* 0.5*mbs < sz < 0.7*mbs */

117:       for (j = 0; j < sz; j++) {
118:         PetscCall(PetscRandomGetValue(rand, &rval));
119:         idx[j * bs] = bs * (PetscInt)(PetscRealPart(rval) * Mbs);
120:         for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
121:       }
122:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
123:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
124:       if (rank == vid) {
125:         PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n", rank, i, sz));
126:         PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
127:       }
128:     } else { /* Test all rows and columns */
129:       sz = M;
130:       PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is1 + i));
131:       PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is2 + i));

133:       if (rank == vid) {
134:         PetscBool colflag;
135:         PetscCall(ISIdentity(is2[i], &colflag));
136:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] is2[%" PetscInt_FMT "], colflag %d\n", rank, i, colflag));
137:         PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
138:       }
139:     }
140:   }

142:   PetscCall(PetscLogStageRegister("MatOv_SBAIJ", &stages[0]));
143:   PetscCall(PetscLogStageRegister("MatOv_BAIJ", &stages[1]));

145:   /* Test MatIncreaseOverlap */
146:   if (TestOverlap) {
147:     PetscCall(PetscLogStagePush(stages[0]));
148:     PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));
149:     PetscCall(PetscLogStagePop());

151:     PetscCall(PetscLogStagePush(stages[1]));
152:     PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
153:     PetscCall(PetscLogStagePop());

155:     if (rank == vid) {
156:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from BAIJ:\n", rank));
157:       PetscCall(ISView(is1[0], PETSC_VIEWER_STDOUT_SELF));
158:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from SBAIJ:\n", rank));
159:       PetscCall(ISView(is2[0], PETSC_VIEWER_STDOUT_SELF));
160:     }

162:     for (i = 0; i < nd; ++i) {
163:       PetscCall(ISEqual(is1[i], is2[i], &flg));
164:       if (!flg) {
165:         if (rank == 0) {
166:           PetscCall(ISSort(is1[i]));
167:           PetscCall(ISSort(is2[i]));
168:         }
169:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i);
170:       }
171:     }
172:   }

174:   /* Test MatCreateSubmatrices */
175:   if (TestSubMat) {
176:     if (test_sorted) {
177:       for (i = 0; i < nd; ++i) PetscCall(ISSort(is1[i]));
178:     }
179:     PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
180:     PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_INITIAL_MATRIX, &submatsA));

182:     PetscCall(MatMultEqual(A, sA, 10, &flg));
183:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A != sA");

185:     /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
186:     PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
187:     PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_REUSE_MATRIX, &submatsA));
188:     PetscCall(MatMultEqual(A, sA, 10, &flg));
189:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatCreateSubmatrices(): A != sA");

191:     PetscCall(MatDestroySubMatrices(nd, &submatA));
192:     PetscCall(MatDestroySubMatrices(nd, &submatsA));
193:   }

195:   /* Free allocated memory */
196:   for (i = 0; i < nd; ++i) {
197:     PetscCall(ISDestroy(&is1[i]));
198:     PetscCall(ISDestroy(&is2[i]));
199:   }
200:   PetscCall(PetscFree(is1));
201:   PetscCall(PetscFree(is2));
202:   PetscCall(PetscFree(idx));
203:   PetscCall(PetscFree(rows));
204:   PetscCall(PetscFree(cols));
205:   PetscCall(PetscFree(vals));
206:   PetscCall(MatDestroy(&A));
207:   PetscCall(MatDestroy(&sA));
208:   PetscCall(PetscRandomDestroy(&rand));
209:   PetscCall(PetscFinalize());
210:   return 0;
211: }

213: /*TEST

215:    test:
216:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
217:       output_file: output/ex92_1.out

219:    test:
220:       suffix: 2
221:       nsize: {{3 4}}
222:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
223:       output_file: output/ex92_1.out

225:    test:
226:       suffix: 3
227:       nsize: {{3 4}}
228:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols
229:       output_file: output/ex92_1.out

231:    test:
232:       suffix: 3_sorted
233:       nsize: {{3 4}}
234:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted
235:       output_file: output/ex92_1.out

237:    test:
238:       suffix: 4
239:       nsize: {{3 4}}
240:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols
241:       output_file: output/ex92_1.out

243: TEST*/