Actual source code: ex4.c

  1: static char help[] = "Creates a matrix, inserts some values, and tests MatCreateSubMatrices() and MatZeroEntries().\n\n";

  3: #include <petscmat.h>
  4: #include <petsc/private/petscimpl.h>

  6: static PetscErrorCode TestSingleISStructureOnly(MPI_Comm comm)
  7: {
  8:   Mat             A;
  9:   Mat            *reuse = NULL, *fresh;
 10:   IS              rows, cols;
 11:   PetscInt        N, rstart, rend, nreuse, nfresh;
 12:   const PetscInt *ia, *ja, *iref, *jref;
 13:   PetscMPIInt     size;
 14:   PetscBool       done, equal;

 16:   PetscFunctionBeginUser;
 17:   PetscCallMPI(MPI_Comm_size(comm, &size));
 18:   N = 2 * size;
 19:   PetscCall(MatCreate(comm, &A));
 20:   PetscCall(MatSetSizes(A, 2, 2, N, N));
 21:   PetscCall(MatSetType(A, MATMPIAIJ));
 22:   PetscCall(MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE));
 23:   PetscCall(MatMPIAIJSetPreallocation(A, 2, NULL, N - 2, NULL));
 24:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 25:   for (PetscInt i = rstart; i < rend; ++i) {
 26:     PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
 27:     PetscCall(MatSetValue(A, i, (i + 1) % N, 1.0, INSERT_VALUES));
 28:   }
 29:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 30:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 31:   PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &rows));
 32:   PetscCall(ISCreateStride(PETSC_COMM_SELF, N / 2, 0, 2, &cols));
 33:   PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
 34:   PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &reuse));
 35:   for (PetscInt iteration = 0; iteration < 2; ++iteration) {
 36:     if (iteration) PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_REUSE_MATRIX, &reuse));
 37:     PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_FALSE));
 38:     PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &fresh));
 39:     PetscCall(MatGetRowIJ(reuse[0], 0, PETSC_FALSE, PETSC_FALSE, &nreuse, &ia, &ja, &done));
 40:     PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot access reused submatrix graph");
 41:     PetscCall(MatGetRowIJ(fresh[0], 0, PETSC_FALSE, PETSC_FALSE, &nfresh, &iref, &jref, &done));
 42:     PetscCheck(done && nreuse == nfresh, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Structure-only submatrix row counts differ");
 43:     PetscCall(PetscArraycmp(ia, iref, nreuse + 1, &equal));
 44:     PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Structure-only submatrix row offsets differ");
 45:     PetscCall(PetscArraycmp(ja, jref, ia[nreuse], &equal));
 46:     PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Structure-only submatrix column indices differ");
 47:     PetscCall(MatRestoreRowIJ(fresh[0], 0, PETSC_FALSE, PETSC_FALSE, &nfresh, &iref, &jref, &done));
 48:     PetscCall(MatRestoreRowIJ(reuse[0], 0, PETSC_FALSE, PETSC_FALSE, &nreuse, &ia, &ja, &done));
 49:     PetscCall(MatDestroySubMatrices(1, &fresh));
 50:   }
 51:   PetscCall(MatDestroySubMatrices(1, &reuse));
 52:   PetscCall(ISDestroy(&rows));
 53:   PetscCall(ISDestroy(&cols));
 54:   PetscCall(MatDestroy(&A));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: static PetscErrorCode TestSingleISOutputGraph(MPI_Comm comm)
 59: {
 60:   Mat              A, difference;
 61:   Mat             *reuse, *fresh;
 62:   IS               rows, cols;
 63:   PetscInt         rstart, rend, row = 0;
 64:   PetscMPIInt      rank;
 65:   PetscBool        equal;
 66:   PetscReal        norm;
 67:   PetscObjectState parentstate, newparentstate, nzstate, newnzstate;

 69:   PetscFunctionBeginUser;
 70:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 71:   PetscCall(MatCreate(comm, &A));
 72:   PetscCall(MatSetSizes(A, 6, 6, PETSC_DETERMINE, PETSC_DETERMINE));
 73:   PetscCall(MatSetType(A, MATMPIAIJ));
 74:   PetscCall(MatMPIAIJSetPreallocation(A, 2, NULL, 0, NULL));
 75:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 76:   for (PetscInt i = rstart; i < rend; ++i) {
 77:     PetscCall(MatSetValue(A, i, i, 10 + i, INSERT_VALUES));
 78:     if (i + 1 < rend) PetscCall(MatSetValue(A, i, i + 1, 20 + i, INSERT_VALUES));
 79:   }
 80:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 81:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 82:   PetscCall(MatGetNonzeroState(A, &parentstate));
 83:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 4, rstart, 1, &rows));
 84:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 4, rstart, 1, &cols));
 85:   for (PetscInt mode = 0; mode < 4; ++mode) {
 86:     PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
 87:     PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &reuse));
 88:     // Exercise graph changes both before the first reuse and after its maps are built.
 89:     if (mode >= 2) PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_REUSE_MATRIX, &reuse));
 90:     PetscCall(MatSetOption(reuse[0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 91:     PetscCall(MatGetNonzeroState(reuse[0], &nzstate));
 92:     // Change only rank 0 so other ranks can build or continue using their cached positions.
 93:     if (!rank) {
 94:       if (!(mode % 2)) PetscCall(MatZeroRows(reuse[0], 1, &row, 0.0, NULL, NULL));
 95:       else PetscCall(MatSetValue(reuse[0], 0, 3, 0.0, INSERT_VALUES));
 96:       PetscCall(MatAssemblyBegin(reuse[0], MAT_FINAL_ASSEMBLY));
 97:       PetscCall(MatAssemblyEnd(reuse[0], MAT_FINAL_ASSEMBLY));
 98:     }
 99:     PetscCall(MatGetNonzeroState(reuse[0], &newnzstate));
100:     PetscCheck(rank ? newnzstate == nzstate : newnzstate != nzstate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Submatrix graph must change only on rank 0");
101:     for (PetscInt iteration = 0; iteration < 2; ++iteration) {
102:       PetscCall(MatScale(A, 2.0));
103:       PetscCall(MatGetNonzeroState(A, &newparentstate));
104:       PetscCheck(newparentstate == parentstate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Parent graph must remain unchanged");
105:       PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_REUSE_MATRIX, &reuse));
106:       PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_FALSE));
107:       PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &fresh));
108:       if (!(mode % 2)) {
109:         PetscCall(MatEqual(reuse[0], fresh[0], &equal));
110:         PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reused submatrix differs after removing a row, iteration %" PetscInt_FMT, iteration);
111:       }
112:       // An inserted explicit zero changes the graph but must not change the numerical result.
113:       PetscCall(MatDuplicate(reuse[0], MAT_COPY_VALUES, &difference));
114:       PetscCall(MatSetOption(difference, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
115:       PetscCall(MatAXPY(difference, -1.0, fresh[0], DIFFERENT_NONZERO_PATTERN));
116:       PetscCall(MatNorm(difference, NORM_INFINITY, &norm));
117:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reused submatrix differs after output graph mutation, mode %" PetscInt_FMT ", iteration %" PetscInt_FMT ": %g", mode, iteration, (double)norm);
118:       PetscCall(MatDestroy(&difference));
119:       PetscCall(MatDestroySubMatrices(1, &fresh));
120:     }
121:     PetscCall(MatDestroySubMatrices(1, &reuse));
122:   }
123:   PetscCall(ISDestroy(&rows));
124:   PetscCall(ISDestroy(&cols));
125:   PetscCall(MatDestroy(&A));
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: static PetscErrorCode TestSingleISReuse(MPI_Comm comm)
130: {
131:   Mat              A, retained, duplicate;
132:   Mat             *reuse, *fresh;
133:   IS               rows, cols;
134:   Vec              x, y, reference;
135:   PetscInt         mlocal = 6, nlocal, M, N, rstart, rend, nr, nc;
136:   PetscInt        *rowidx, *colidx;
137:   PetscMPIInt      rank, size;
138:   PetscBool        empty_rank = PETSC_FALSE, rectangular = PETSC_FALSE, equal;
139:   PetscReal        norm;
140:   PetscObjectState state, newstate, nzstate, newnzstate;

142:   PetscFunctionBeginUser;
143:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
144:   PetscCallMPI(MPI_Comm_size(comm, &size));
145:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_empty_rank", &empty_rank, NULL));
146:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_rectangular", &rectangular, NULL));
147:   if (empty_rank && size > 1 && rank == size - 1) mlocal = 0;
148:   nlocal = rectangular && mlocal ? 4 : mlocal;
149:   M      = 6 * (empty_rank && size > 1 ? size - 1 : size);
150:   N      = (rectangular ? 4 : 6) * (empty_rank && size > 1 ? size - 1 : size);
151:   PetscCall(MatCreate(comm, &A));
152:   PetscCall(MatSetSizes(A, mlocal, nlocal, M, N));
153:   PetscCall(MatSetType(A, MATMPIAIJ));
154:   PetscCall(MatSetOptionsPrefix(A, "reuse_"));
155:   PetscCall(MatSetFromOptions(A));
156:   PetscCall(MatMPIAIJSetPreallocation(A, nlocal, NULL, N - nlocal, NULL));
157:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
158:   for (PetscInt i = rstart; i < rend; ++i) {
159:     for (PetscInt j = 0; j < N; ++j) {
160:       PetscScalar value = 5 * (i + 1) + 3 * (j + 1);

162: #if PetscDefined(USE_COMPLEX)
163:       value += PETSC_i * (i - 2 * j + 1);
164: #endif
165:       if ((i + j) % 3 != 1 || i == j) PetscCall(MatSetValue(A, i, j, value, INSERT_VALUES));
166:     }
167:   }
168:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
169:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
170:   PetscCall(PetscMalloc2(M, &rowidx, N, &colidx));

172:   // Include fallback selections, empty graphs, remote-only rows, and a proper subset of owned rows.
173:   for (PetscInt scenario = 0; scenario < 9; ++scenario) {
174:     nr = M;
175:     nc = N / 2;
176:     for (PetscInt i = 0; i < nr; ++i) rowidx[i] = scenario == 1 ? M - 1 - i : i;
177:     for (PetscInt i = 0; i < nc; ++i) colidx[i] = 2 * (scenario == 2 ? nc - 1 - i : i);
178:     if (scenario == 3) {
179:       nr = M / 2;
180:       nc = N;
181:       for (PetscInt i = 0; i < nr; ++i) rowidx[i] = 2 * i;
182:       for (PetscInt i = 0; i < nc; ++i) colidx[i] = i;
183:     }
184:     if (scenario == 4 && !rank) nr = nc = 0;
185:     if (scenario == 5) {
186:       nr = 0;
187:       for (PetscInt i = 0; i < M; ++i) {
188:         if (i < rstart || i >= rend) rowidx[nr++] = i;
189:       }
190:     }
191:     if (scenario == 6) {
192:       nr = M / 2;
193:       for (PetscInt i = 0; i < nr; ++i) rowidx[i] = 2 * i;
194:     }
195:     if (scenario == 7) nc = 0;
196:     if (scenario == 8) {
197:       nr = nc   = 1;
198:       rowidx[0] = 0;
199:       colidx[0] = 1; // This entry is absent from the parent graph.
200:     }
201:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nr, rowidx, PETSC_COPY_VALUES, &rows));
202:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nc, colidx, PETSC_COPY_VALUES, &cols));
203:     PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
204:     PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &reuse));
205:     PetscCall(MatCreateVecs(reuse[0], &x, &y));
206:     PetscCall(VecDuplicate(y, &reference));
207:     PetscCall(VecSet(x, 1.0));
208:     for (PetscInt iteration = 0; iteration < 3; ++iteration) {
209:       // Exercise cached/device values before changing the parent.
210:       PetscCall(MatMult(reuse[0], x, y));
211:       PetscCall(PetscObjectStateGet((PetscObject)reuse[0], &state));
212:       PetscCall(MatGetNonzeroState(reuse[0], &nzstate));
213:       PetscCall(MatScale(A, -0.5));
214:       if (!rectangular) PetscCall(MatShift(A, iteration + 1));
215:       PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_REUSE_MATRIX, &reuse));
216:       PetscCall(PetscObjectStateGet((PetscObject)reuse[0], &newstate));
217:       PetscCall(MatGetNonzeroState(reuse[0], &newnzstate));
218:       PetscCheck(newstate > state && newnzstate == nzstate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reusing a submatrix must update its values without changing its graph");
219:       // A fresh extraction uses the general path, independently of the cached SingleIS maps.
220:       PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_FALSE));
221:       PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &fresh));
222:       PetscCall(MatEqual(reuse[0], fresh[0], &equal));
223:       PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reused submatrix differs from fresh extraction in scenario %" PetscInt_FMT ", iteration %" PetscInt_FMT, scenario, iteration);
224:       PetscCall(MatMult(reuse[0], x, y));
225:       PetscCall(MatMult(fresh[0], x, reference));
226:       PetscCall(VecAXPY(y, -1.0, reference));
227:       PetscCall(VecNorm(y, NORM_INFINITY, &norm));
228:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reused submatrix product differs from fresh extraction: %g", (double)norm);
229:       PetscCall(MatDestroySubMatrices(1, &fresh));
230:     }
231:     retained = reuse[0];
232:     PetscCall(PetscObjectReference((PetscObject)retained));
233:     PetscCall(MatDuplicate(retained, MAT_COPY_VALUES, &duplicate));
234:     PetscCall(MatDestroySubMatrices(1, &reuse));
235:     PetscCall(MatEqual(retained, duplicate, &equal));
236:     PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Retained submatrix values changed after destroying its containing array");
237:     PetscCall(MatMult(retained, x, y));
238:     PetscCall(MatDestroy(&duplicate));
239:     PetscCall(MatDestroy(&retained));
240:     PetscCall(VecDestroy(&x));
241:     PetscCall(VecDestroy(&y));
242:     PetscCall(VecDestroy(&reference));
243:     PetscCall(ISDestroy(&rows));
244:     PetscCall(ISDestroy(&cols));
245:   }
246:   PetscCall(PetscFree2(rowidx, colidx));
247:   PetscCall(MatDestroy(&A));
248:   PetscCall(TestSingleISStructureOnly(comm));
249:   PetscCall(TestSingleISOutputGraph(comm));
250:   PetscCall(PetscPrintf(comm, "Submatrix reuse tests passed\n"));
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

254: int main(int argc, char **argv)
255: {
256:   Mat         mat, submat, submat1;
257:   Mat        *submatrices;
258:   PetscInt    m = 10, n = 10, i = 4, tmp, rstart, rend;
259:   IS          irow, icol;
260:   PetscScalar value = 1.0;
261:   PetscViewer sviewer;
262:   PetscBool   allA = PETSC_FALSE, test_reuse = PETSC_FALSE;

264:   PetscFunctionBeginUser;
265:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
266:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_reuse", &test_reuse, NULL));
267:   if (test_reuse) PetscCall(TestSingleISReuse(PETSC_COMM_WORLD));
268:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
269:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_COMMON));

271:   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
272:   PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
273:   PetscCall(MatSetFromOptions(mat));
274:   PetscCall(MatSetUp(mat));
275:   PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
276:   for (i = rstart; i < rend; i++) {
277:     value = (PetscReal)i + 1;
278:     tmp   = i % 5;
279:     PetscCall(MatSetValues(mat, 1, &tmp, 1, &i, &value, INSERT_VALUES));
280:   }
281:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
282:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
283:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original matrix\n"));
284:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));

286:   /* Test MatCreateSubMatrix_XXX_All(), i.e., submatrix = A */
287:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_all", &allA, NULL));
288:   if (allA) {
289:     PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &irow));
290:     PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &icol));
291:     PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
292:     PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
293:     submat = *submatrices;

295:     /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
296:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with all\n"));
297:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "--------------------\n"));
298:     PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
299:     PetscCall(MatView(submat, sviewer));
300:     PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));

302:     PetscCall(ISDestroy(&irow));
303:     PetscCall(ISDestroy(&icol));

305:     /* test getting a reference on a submat */
306:     PetscCall(PetscObjectReference((PetscObject)submat));
307:     PetscCall(MatDestroySubMatrices(1, &submatrices));
308:     PetscCall(MatDestroy(&submat));
309:   }

311:   /* Form submatrix with rows 2-4 and columns 4-8 */
312:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 2, 1, &irow));
313:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 5, 4, 1, &icol));
314:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
315:   submat = *submatrices;

317:   /* Test reuse submatrices */
318:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));

320:   /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
321:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices\n"));
322:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
323:   PetscCall(MatView(submat, sviewer));
324:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
325:   PetscCall(PetscObjectReference((PetscObject)submat));
326:   PetscCall(MatDestroySubMatrices(1, &submatrices));
327:   PetscCall(MatDestroy(&submat));

329:   /* Form submatrix with rows 2-4 and all columns */
330:   PetscCall(ISDestroy(&icol));
331:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 10, 0, 1, &icol));
332:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
333:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
334:   submat = *submatrices;

336:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with allcolumns\n"));
337:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
338:   PetscCall(MatView(submat, sviewer));
339:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));

341:   /* Test MatDuplicate */
342:   PetscCall(MatDuplicate(submat, MAT_COPY_VALUES, &submat1));
343:   PetscCall(MatDestroy(&submat1));

345:   /* Zero the original matrix */
346:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original zeroed matrix\n"));
347:   PetscCall(MatZeroEntries(mat));
348:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));

350:   PetscCall(ISDestroy(&irow));
351:   PetscCall(ISDestroy(&icol));
352:   PetscCall(PetscObjectReference((PetscObject)submat));
353:   PetscCall(MatDestroySubMatrices(1, &submatrices));
354:   PetscCall(MatDestroy(&submat));
355:   PetscCall(MatDestroy(&mat));
356:   PetscCall(PetscFinalize());
357:   return 0;
358: }

360: /*TEST

362:    test:
363:       args: -mat_type aij

365:    test:
366:       suffix: 2
367:       args: -mat_type dense

369:    test:
370:       suffix: 3
371:       nsize: 3
372:       args: -mat_type aij

374:    test:
375:       suffix: 4
376:       nsize: 3
377:       args: -mat_type dense

379:    test:
380:       suffix: 5
381:       nsize: 3
382:       args: -mat_type aij -test_all

384:    test:
385:       suffix: singleis_reuse
386:       nsize: {{1 2 3}}
387:       args: -mat_type aij -test_reuse -test_empty_rank {{0 1}}
388:       filter: grep "^Submatrix reuse"
389:       output_file: output/ex4_singleis_reuse.out

391:    test:
392:       suffix: singleis_reuse_rectangular
393:       nsize: 3
394:       args: -mat_type aij -test_reuse -test_rectangular -test_empty_rank {{0 1}}
395:       filter: grep "^Submatrix reuse"
396:       output_file: output/ex4_singleis_reuse.out

398:    test:
399:       suffix: singleis_reuse_cuda
400:       requires: cuda
401:       nsize: 3
402:       args: -mat_type aij -test_reuse -test_empty_rank {{0 1}} -reuse_mat_type mpiaijcusparse
403:       filter: grep "^Submatrix reuse"
404:       output_file: output/ex4_singleis_reuse.out

406:    test:
407:       suffix: singleis_reuse_kokkos
408:       requires: kokkos_kernels
409:       nsize: 3
410:       args: -mat_type aij -test_reuse -test_empty_rank {{0 1}} -reuse_mat_type mpiaijkokkos
411:       filter: grep "^Submatrix reuse"
412:       output_file: output/ex4_singleis_reuse.out

414: TEST*/