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, (char *)0, 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:     for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1;
128:   }
129:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
130:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

132:   if (submat) {
133:     Mat      A2;
134:     IS       r, c;
135:     PetscInt rst, ren, cst, cen;

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

148:   PetscCall(MatGetSize(A, &M, &N));
149:   PetscCall(MatGetLocalSize(A, &m, &n));
150:   PetscCall(MatHasCongruentLayouts(A, &cong));

152:   PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
153:   PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, keep));
154:   PetscCall(PetscObjectSetName((PetscObject)A, "initial"));
155:   PetscCall(MatViewFromOptions(A, NULL, "-view_mat"));

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

166:   /* Square and rows only scaling */
167:   ronl = cong ? ronl : PETSC_TRUE;

169:   for (test = 0; test < ntest; test++) {
170:     PetscReal err;

172:     PetscCall(MatMultAddEqual(A, S, 10, &flg));
173:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add\n", test));
174:     PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg));
175:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add (T)\n", test));
176:     if (testzerorows) {
177:       Mat       ST, B, C, BT, BTT;
178:       IS        zr;
179:       Vec       x = NULL, b1 = NULL, b2 = NULL;
180:       PetscInt *idxs = NULL, nr = 0;

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

238:       if (b1) {
239:         Vec b;

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

255:       PetscCall(MatCreateTranspose(S, &ST));
256:       PetscCall(MatComputeOperator(ST, MATDENSE, &BT));
257:       PetscCall(MatTranspose(BT, MAT_INITIAL_MATRIX, &BTT));
258:       PetscCall(PetscObjectSetName((PetscObject)B, "S"));
259:       PetscCall(PetscObjectSetName((PetscObject)BTT, "STT"));
260:       PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &C));
261:       PetscCall(PetscObjectSetName((PetscObject)C, "A"));

263:       PetscCall(MatViewFromOptions(C, NULL, "-view_mat"));
264:       PetscCall(MatViewFromOptions(B, NULL, "-view_mat"));
265:       PetscCall(MatViewFromOptions(BTT, NULL, "-view_mat"));

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

271:       PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &C));
272:       PetscCall(MatAXPY(C, -1.0, BTT, SAME_NONZERO_PATTERN));
273:       PetscCall(MatNorm(C, NORM_FROBENIUS, &err));
274:       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));

276:       PetscCall(MatDestroy(&ST));
277:       PetscCall(MatDestroy(&BTT));
278:       PetscCall(MatDestroy(&BT));
279:       PetscCall(MatDestroy(&B));
280:       PetscCall(MatDestroy(&C));
281:     }
282:     if (testdiagscale) { /* MatDiagonalScale() */
283:       Vec vr, vl;

285:       PetscCall(MatCreateVecs(A, &vr, &vl));
286:       if (randomize) {
287:         PetscCall(VecSetRandom(vr, NULL));
288:         PetscCall(VecSetRandom(vl, NULL));
289:       } else {
290:         PetscCall(VecSet(vr, test % 2 ? 0.15 : 1.0 / 0.15));
291:         PetscCall(VecSet(vl, test % 2 ? -1.2 : 1.0 / -1.2));
292:       }
293:       PetscCall(MatDiagonalScale(A, vl, vr));
294:       PetscCall(MatDiagonalScale(S, vl, vr));
295:       PetscCall(VecDestroy(&vr));
296:       PetscCall(VecDestroy(&vl));
297:     }

299:     if (testscale) { /* MatScale() */
300:       PetscCall(MatScale(A, test % 2 ? 1.4 : 1.0 / 1.4));
301:       PetscCall(MatScale(S, test % 2 ? 1.4 : 1.0 / 1.4));
302:     }

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

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

312:       PetscCall(MatCreateVecs(A, &dA, NULL));
313:       PetscCall(MatCreateVecs(S, &dS, NULL));
314:       PetscCall(MatGetDiagonal(A, dA));
315:       PetscCall(MatGetDiagonal(S, dS));
316:       PetscCall(VecAXPY(dA, -1.0, dS));
317:       PetscCall(VecNorm(dA, NORM_INFINITY, &err));
318:       if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error diag %g\n", test, (double)err));
319:       PetscCall(VecDestroy(&dA));
320:       PetscCall(VecDestroy(&dS));
321:     }

323:     if (testdup && !test) {
324:       Mat A2, S2;

326:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
327:       PetscCall(MatDuplicate(S, MAT_COPY_VALUES, &S2));
328:       PetscCall(MatDestroy(&A));
329:       PetscCall(MatDestroy(&S));
330:       A = A2;
331:       S = S2;
332:     }

334:     if (testsubmat) {
335:       Mat      sA, sS, dA, dS, At, St;
336:       IS       r, c;
337:       PetscInt rst, ren, cst, cen;

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

381:     if (testaxpy) {
382:       Mat          tA, tS, dA, dS;
383:       MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN};

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

412:     if (testreset && (ntest == 1 || test == ntest - 2)) {
413:       /* reset MATSHELL */
414:       PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
415:       PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
416:       /* reset A */
417:       PetscCall(MatCopy(user->B, A, DIFFERENT_NONZERO_PATTERN));
418:     }
419:   }

421:   PetscCall(MatDestroy(&A));
422:   PetscCall(MatDestroy(&S));
423:   PetscCall(PetscFree(user));
424:   PetscCall(PetscFinalize());
425:   return 0;
426: }

428: /*TEST

430:    testset:
431:      suffix: rect
432:      requires: !single
433:      output_file: output/ex221_1.out
434:      nsize: {{1 3}}
435:      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}

437:    testset:
438:      suffix: square
439:      requires: !single
440:      output_file: output/ex221_1.out
441:      nsize: {{1 3}}
442:      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
443: TEST*/