Actual source code: ex221.c

  1: static char help[] = "Tests various routines for MATSHELL\n\n";

  3: #include <petscmat.h>

  5: typedef struct _n_User *User;
  6: struct _n_User {
  7:   Mat B;
  8: };

 10: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
 11: {
 12:   User user;

 14:   PetscFunctionBegin;
 15:   PetscCall(MatShellGetContext(A, &user));
 16:   PetscCall(MatGetDiagonal(user->B, X));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
 21: {
 22:   User user;

 24:   PetscFunctionBegin;
 25:   PetscCall(MatShellGetContext(A, &user));
 26:   PetscCall(MatMult(user->B, X, Y));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
 31: {
 32:   User user;

 34:   PetscFunctionBegin;
 35:   PetscCall(MatShellGetContext(A, &user));
 36:   PetscCall(MatMultTranspose(user->B, X, Y));
 37:   PetscFunctionReturn(PETSC_SUCCESS);
 38: }

 40: static PetscErrorCode MatCopy_User(Mat A, Mat X, MatStructure str)
 41: {
 42:   User user, userX;

 44:   PetscFunctionBegin;
 45:   PetscCall(MatShellGetContext(A, &user));
 46:   PetscCall(MatShellGetContext(X, &userX));
 47:   PetscCheck(user == userX, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
 48:   PetscCall(PetscObjectReference((PetscObject)user->B));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: static PetscErrorCode MatDestroy_User(Mat A)
 53: {
 54:   User user;

 56:   PetscFunctionBegin;
 57:   PetscCall(MatShellGetContext(A, &user));
 58:   PetscCall(PetscObjectDereference((PetscObject)user->B));
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: int main(int argc, char **args)
 63: {
 64:   User         user;
 65:   Mat          A, S;
 66:   PetscScalar *data, diag = 1.3;
 67:   PetscReal    tol = PETSC_SMALL;
 68:   PetscInt     i, j, m = PETSC_DECIDE, n = PETSC_DECIDE, M = 17, N = 15, s1, s2;
 69:   PetscInt     test, ntest = 2;
 70:   PetscMPIInt  rank, size;
 71:   PetscBool    nc        = PETSC_FALSE, cong, flg;
 72:   PetscBool    ronl      = PETSC_TRUE;
 73:   PetscBool    randomize = PETSC_FALSE, submat = PETSC_FALSE;
 74:   PetscBool    keep         = PETSC_FALSE;
 75:   PetscBool    testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
 76:   PetscBool    testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
 77:   PetscBool    testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;

 79:   PetscFunctionBeginUser;
 80:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 81:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 82:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 83:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 84:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 85:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ml", &m, NULL));
 86:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &n, NULL));
 87:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-square_nc", &nc, NULL));
 88:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-rows_only", &ronl, NULL));
 89:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize", &randomize, NULL));
 90:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-submat", &submat, NULL));
 91:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_zerorows", &testzerorows, NULL));
 92:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_diagscale", &testdiagscale, NULL));
 93:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_getdiag", &testgetdiag, NULL));
 94:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_shift", &testshift, NULL));
 95:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_scale", &testscale, NULL));
 96:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_dup", &testdup, NULL));
 97:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_reset", &testreset, NULL));
 98:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_submat", &testsubmat, NULL));
 99:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy", &testaxpy, NULL));
100:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy_different", &testaxpyd, NULL));
101:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy_error", &testaxpyerr, NULL));
102:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-loop", &ntest, NULL));
103:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
104:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-diag", &diag, NULL));
105:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-keep", &keep, NULL));
106:   /* This tests square matrices with different row/col layout */
107:   if (nc && size > 1) {
108:     M = PetscMax(PetscMax(N, M), 1);
109:     N = M;
110:     m = n = 0;
111:     if (rank == 0) {
112:       m = M - 1;
113:       n = 1;
114:     } else if (rank == 1) {
115:       m = 1;
116:       n = N - 1;
117:     }
118:   }
119:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, NULL, &A));
120:   PetscCall(MatGetLocalSize(A, &m, &n));
121:   PetscCall(MatGetSize(A, &M, &N));
122:   PetscCall(MatGetOwnershipRange(A, &s1, NULL));
123:   s2 = 1;
124:   while (s2 < M) s2 *= 10;
125:   PetscCall(MatDenseGetArray(A, &data));
126:   for (j = 0; j < N; j++) {
127: #if defined(PETSC_USE_COMPLEX)
128:     for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1 + PETSC_i * (s1 - 1);
129: #else
130:     for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1;
131: #endif
132:   }
133:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
134:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

136:   if (submat) {
137:     Mat      A2;
138:     IS       r, c;
139:     PetscInt rst, ren, cst, cen;

141:     PetscCall(MatGetOwnershipRange(A, &rst, &ren));
142:     PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
143:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r));
144:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c));
145:     PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &A2));
146:     PetscCall(ISDestroy(&r));
147:     PetscCall(ISDestroy(&c));
148:     PetscCall(MatDestroy(&A));
149:     A = A2;
150:   }

152:   PetscCall(MatGetSize(A, &M, &N));
153:   PetscCall(MatGetLocalSize(A, &m, &n));
154:   PetscCall(MatHasCongruentLayouts(A, &cong));

156:   PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
157:   PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, keep));
158:   PetscCall(PetscObjectSetName((PetscObject)A, "initial"));
159:   PetscCall(MatViewFromOptions(A, NULL, "-view_mat"));

161:   PetscCall(PetscNew(&user));
162:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, user, &S));
163:   PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User));
164:   PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User));
165:   if (cong) PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User));
166:   PetscCall(MatShellSetOperation(S, MATOP_COPY, (void (*)(void))MatCopy_User));
167:   PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_User));
168:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &user->B));

170:   /* Square and rows only scaling */
171:   ronl = cong ? ronl : PETSC_TRUE;

173:   for (test = 0; test < ntest; test++) {
174:     PetscReal err;

176:     PetscCall(MatMultAddEqual(A, S, 10, &flg));
177:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add\n", test));
178:     PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg));
179:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult transpose add\n", test));
180:     PetscCall(MatMultHermitianTransposeAddEqual(A, S, 10, &flg));
181:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult hermitian transpose add\n", test));
182:     if (testzerorows) {
183:       Mat       ST, B, C, BT, BTT;
184:       IS        zr;
185:       Vec       x = NULL, b1 = NULL, b2 = NULL;
186:       PetscInt *idxs = NULL, nr = 0;

188:       if (rank == (test % size)) {
189:         nr = 1;
190:         PetscCall(PetscMalloc1(nr, &idxs));
191:         if (test % 2) {
192:           idxs[0] = (2 * M - 1 - test / 2) % M;
193:         } else {
194:           idxs[0] = (test / 2) % M;
195:         }
196:         idxs[0] = PetscMax(idxs[0], 0);
197:       }
198:       PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr));
199:       PetscCall(PetscObjectSetName((PetscObject)zr, "ZR"));
200:       PetscCall(ISViewFromOptions(zr, NULL, "-view_is"));
201:       PetscCall(MatCreateVecs(A, &x, &b1));
202:       if (randomize) {
203:         PetscCall(VecSetRandom(x, NULL));
204:         PetscCall(VecSetRandom(b1, NULL));
205:       } else {
206:         PetscCall(VecSet(x, 11.4));
207:         PetscCall(VecSet(b1, -14.2));
208:       }
209:       PetscCall(VecDuplicate(b1, &b2));
210:       PetscCall(VecCopy(b1, b2));
211:       PetscCall(PetscObjectSetName((PetscObject)b1, "A_B1"));
212:       PetscCall(PetscObjectSetName((PetscObject)b2, "A_B2"));
213:       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
214:         PetscCall(VecDestroy(&b1));
215:       }
216:       if (ronl) {
217:         PetscCall(MatZeroRowsIS(A, zr, diag, x, b1));
218:         PetscCall(MatZeroRowsIS(S, zr, diag, x, b2));
219:       } else {
220:         PetscCall(MatZeroRowsColumnsIS(A, zr, diag, x, b1));
221:         PetscCall(MatZeroRowsColumnsIS(S, zr, diag, x, b2));
222:         PetscCall(ISDestroy(&zr));
223:         /* Mix zerorows and zerorowscols */
224:         nr   = 0;
225:         idxs = NULL;
226:         if (rank == 0) {
227:           nr = 1;
228:           PetscCall(PetscMalloc1(nr, &idxs));
229:           if (test % 2) {
230:             idxs[0] = (3 * M - 2 - test / 2) % M;
231:           } else {
232:             idxs[0] = (test / 2 + 1) % M;
233:           }
234:           idxs[0] = PetscMax(idxs[0], 0);
235:         }
236:         PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr));
237:         PetscCall(PetscObjectSetName((PetscObject)zr, "ZR2"));
238:         PetscCall(ISViewFromOptions(zr, NULL, "-view_is"));
239:         PetscCall(MatZeroRowsIS(A, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL));
240:         PetscCall(MatZeroRowsIS(S, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL));
241:       }
242:       PetscCall(ISDestroy(&zr));

244:       if (b1) {
245:         Vec b;

247:         PetscCall(VecViewFromOptions(b1, NULL, "-view_b"));
248:         PetscCall(VecViewFromOptions(b2, NULL, "-view_b"));
249:         PetscCall(VecDuplicate(b1, &b));
250:         PetscCall(VecCopy(b1, b));
251:         PetscCall(VecAXPY(b, -1.0, b2));
252:         PetscCall(VecNorm(b, NORM_INFINITY, &err));
253:         if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error b %g\n", test, (double)err));
254:         PetscCall(VecDestroy(&b));
255:       }
256:       PetscCall(VecDestroy(&b1));
257:       PetscCall(VecDestroy(&b2));
258:       PetscCall(VecDestroy(&x));
259:       PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &B));

261:       PetscCall(MatCreateTranspose(S, &ST));
262:       PetscCall(MatComputeOperator(ST, MATDENSE, &BT));
263:       PetscCall(MatTranspose(BT, MAT_INITIAL_MATRIX, &BTT));
264:       PetscCall(PetscObjectSetName((PetscObject)B, "S"));
265:       PetscCall(PetscObjectSetName((PetscObject)BTT, "STT"));
266:       PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &C));
267:       PetscCall(PetscObjectSetName((PetscObject)C, "A"));

269:       PetscCall(MatViewFromOptions(C, NULL, "-view_mat"));
270:       PetscCall(MatViewFromOptions(B, NULL, "-view_mat"));
271:       PetscCall(MatViewFromOptions(BTT, NULL, "-view_mat"));

273:       PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN));
274:       PetscCall(MatNorm(C, NORM_FROBENIUS, &err));
275:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err));

277:       PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &C));
278:       PetscCall(MatAXPY(C, -1.0, BTT, SAME_NONZERO_PATTERN));
279:       PetscCall(MatNorm(C, NORM_FROBENIUS, &err));
280:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err));

282:       PetscCall(MatDestroy(&ST));
283:       PetscCall(MatDestroy(&BTT));
284:       PetscCall(MatDestroy(&BT));
285:       PetscCall(MatDestroy(&B));
286:       PetscCall(MatDestroy(&C));
287:     }
288:     if (testdiagscale) { /* MatDiagonalScale() */
289:       Vec vr, vl;

291:       PetscCall(MatCreateVecs(A, &vr, &vl));
292:       if (randomize) {
293:         PetscCall(VecSetRandom(vr, NULL));
294:         PetscCall(VecSetRandom(vl, NULL));
295:       } else {
296:         PetscCall(VecSet(vr, test % 2 ? 0.15 : 1.0 / 0.15));
297:         PetscCall(VecSet(vl, test % 2 ? -1.2 : 1.0 / -1.2));
298:       }
299:       PetscCall(MatDiagonalScale(A, vl, vr));
300:       PetscCall(MatDiagonalScale(S, vl, vr));
301:       PetscCall(VecDestroy(&vr));
302:       PetscCall(VecDestroy(&vl));
303:     }

305:     if (testscale) { /* MatScale() */
306:       PetscCall(MatScale(A, test % 2 ? 1.4 : 1.0 / 1.4));
307:       PetscCall(MatScale(S, test % 2 ? 1.4 : 1.0 / 1.4));
308:     }

310:     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
311:       PetscCall(MatShift(A, test % 2 ? -77.5 : 77.5));
312:       PetscCall(MatShift(S, test % 2 ? -77.5 : 77.5));
313:     }

315:     if (testgetdiag && cong) { /* MatGetDiagonal() */
316:       Vec dA, dS;

318:       PetscCall(MatCreateVecs(A, &dA, NULL));
319:       PetscCall(MatCreateVecs(S, &dS, NULL));
320:       PetscCall(MatGetDiagonal(A, dA));
321:       PetscCall(MatGetDiagonal(S, dS));
322:       PetscCall(VecAXPY(dA, -1.0, dS));
323:       PetscCall(VecNorm(dA, NORM_INFINITY, &err));
324:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error diag %g\n", test, (double)err));
325:       PetscCall(VecDestroy(&dA));
326:       PetscCall(VecDestroy(&dS));
327:     }

329:     if (testdup && !test) {
330:       Mat A2, S2;

332:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
333:       PetscCall(MatDuplicate(S, MAT_COPY_VALUES, &S2));
334:       PetscCall(MatDestroy(&A));
335:       PetscCall(MatDestroy(&S));
336:       A = A2;
337:       S = S2;
338:     }

340:     if (testsubmat) {
341:       Mat      sA, sS, dA, dS, At, St;
342:       IS       r, c;
343:       PetscInt rst, ren, cst, cen;

345:       PetscCall(MatGetOwnershipRange(A, &rst, &ren));
346:       PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
347:       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r));
348:       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c));
349:       PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &sA));
350:       PetscCall(MatCreateSubMatrix(S, r, c, MAT_INITIAL_MATRIX, &sS));
351:       PetscCall(MatMultAddEqual(sA, sS, 10, &flg));
352:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add\n", test));
353:       PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg));
354:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add (T)\n", test));
355:       PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA));
356:       PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS));
357:       PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN));
358:       PetscCall(MatNorm(dA, NORM_FROBENIUS, &err));
359:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err));
360:       PetscCall(MatDestroy(&sA));
361:       PetscCall(MatDestroy(&sS));
362:       PetscCall(MatDestroy(&dA));
363:       PetscCall(MatDestroy(&dS));
364:       PetscCall(MatCreateTranspose(A, &At));
365:       PetscCall(MatCreateTranspose(S, &St));
366:       PetscCall(MatCreateSubMatrix(At, c, r, MAT_INITIAL_MATRIX, &sA));
367:       PetscCall(MatCreateSubMatrix(St, c, r, MAT_INITIAL_MATRIX, &sS));
368:       PetscCall(MatMultAddEqual(sA, sS, 10, &flg));
369:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add\n", test));
370:       PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg));
371:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n", test));
372:       PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA));
373:       PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS));
374:       PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN));
375:       PetscCall(MatNorm(dA, NORM_FROBENIUS, &err));
376:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n", test, (double)err));
377:       PetscCall(MatDestroy(&sA));
378:       PetscCall(MatDestroy(&sS));
379:       PetscCall(MatDestroy(&dA));
380:       PetscCall(MatDestroy(&dS));
381:       PetscCall(MatDestroy(&At));
382:       PetscCall(MatDestroy(&St));
383:       PetscCall(ISDestroy(&r));
384:       PetscCall(ISDestroy(&c));
385:     }

387:     if (testaxpy) {
388:       Mat          tA, tS, dA, dS;
389:       MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN};

391:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &tA));
392:       if (testaxpyd && !(test % 2)) {
393:         PetscCall(PetscObjectReference((PetscObject)tA));
394:         tS = tA;
395:       } else {
396:         PetscCall(PetscObjectReference((PetscObject)S));
397:         tS = S;
398:       }
399:       PetscCall(MatAXPY(A, 0.5, tA, str[test % 3]));
400:       PetscCall(MatAXPY(S, 0.5, tS, str[test % 3]));
401:       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
402:       if (testaxpyerr) PetscCall(MatScale(tA, 0));
403:       PetscCall(MatDestroy(&tA));
404:       PetscCall(MatDestroy(&tS));
405:       PetscCall(MatMultAddEqual(A, S, 10, &flg));
406:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add\n", test));
407:       PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg));
408:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add (T)\n", test));
409:       PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &dA));
410:       PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &dS));
411:       PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN));
412:       PetscCall(MatNorm(dA, NORM_FROBENIUS, &err));
413:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err));
414:       PetscCall(MatDestroy(&dA));
415:       PetscCall(MatDestroy(&dS));
416:     }

418:     if (testreset && (ntest == 1 || test == ntest - 2)) {
419:       /* reset MATSHELL */
420:       PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
421:       PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
422:       /* reset A */
423:       PetscCall(MatCopy(user->B, A, DIFFERENT_NONZERO_PATTERN));
424:     }
425:   }

427:   PetscCall(MatDestroy(&A));
428:   PetscCall(MatDestroy(&S));
429:   PetscCall(PetscFree(user));
430:   PetscCall(PetscFinalize());
431:   return 0;
432: }

434: /*TEST

436:    testset:
437:      suffix: rect
438:      requires: !single
439:      output_file: output/ex221_1.out
440:      nsize: {{1 3}}
441:      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}

443:    testset:
444:      suffix: square
445:      requires: !single
446:      output_file: output/ex221_1.out
447:      nsize: {{1 3}}
448:      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
449: TEST*/