Actual source code: ex92.c


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

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

 23:   PetscInitialize(&argc, &args, (char *)0, help);
 24:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 25:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

 37:   MatCreate(PETSC_COMM_WORLD, &A);
 38:   MatSetSizes(A, mbs * bs, mbs * bs, PETSC_DECIDE, PETSC_DECIDE);
 39:   MatSetType(A, MATBAIJ);
 40:   MatSeqBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL);
 41:   MatMPIBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);

 43:   PetscRandomCreate(PETSC_COMM_WORLD, &rand);
 44:   PetscRandomSetFromOptions(rand);

 46:   MatGetOwnershipRange(A, &rstart, &rend);
 47:   MatGetSize(A, &M, &N);
 48:   Mbs = M / bs;

 50:   PetscMalloc1(bs, &rows);
 51:   PetscMalloc1(bs, &cols);
 52:   PetscMalloc1(bs * bs, &vals);
 53:   PetscMalloc1(M, &idx);

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

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

 84:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 85:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

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

107:   /* Test sA==A through MatMult() */
108:   MatMultEqual(A, sA, 10, &flg);

111:   /* Test MatIncreaseOverlap() */
112:   PetscMalloc1(nd, &is1);
113:   PetscMalloc1(nd, &is2);

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

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

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

145:   PetscLogStageRegister("MatOv_SBAIJ", &stages[0]);
146:   PetscLogStageRegister("MatOv_BAIJ", &stages[1]);

148:   /* Test MatIncreaseOverlap */
149:   if (TestOverlap) {
150:     PetscLogStagePush(stages[0]);
151:     MatIncreaseOverlap(sA, nd, is2, ov);
152:     PetscLogStagePop();

154:     PetscLogStagePush(stages[1]);
155:     MatIncreaseOverlap(A, nd, is1, ov);
156:     PetscLogStagePop();

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

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

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

185:     MatMultEqual(A, sA, 10, &flg);

188:     /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
189:     MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
190:     MatCreateSubMatrices(sA, nd, is1, is1, MAT_REUSE_MATRIX, &submatsA);
191:     MatMultEqual(A, sA, 10, &flg);

194:     MatDestroySubMatrices(nd, &submatA);
195:     MatDestroySubMatrices(nd, &submatsA);
196:   }

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

216: /*TEST

218:    test:
219:       args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
220:       output_file: output/ex92_1.out

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

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

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

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

246: TEST*/