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*/