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