1 static char help[] = "Tests various routines for MATSHELL\n\n";
2
3 #include <petscmat.h>
4
5 typedef struct _n_User *User;
6 struct _n_User {
7 Mat B;
8 };
9
MatGetDiagonal_User(Mat A,Vec X)10 static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
11 {
12 User user;
13
14 PetscFunctionBegin;
15 PetscCall(MatShellGetContext(A, &user));
16 PetscCall(MatGetDiagonal(user->B, X));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
MatMult_User(Mat A,Vec X,Vec Y)20 static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
21 {
22 User user;
23
24 PetscFunctionBegin;
25 PetscCall(MatShellGetContext(A, &user));
26 PetscCall(MatMult(user->B, X, Y));
27 PetscFunctionReturn(PETSC_SUCCESS);
28 }
29
MatMultTranspose_User(Mat A,Vec X,Vec Y)30 static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
31 {
32 User user;
33
34 PetscFunctionBegin;
35 PetscCall(MatShellGetContext(A, &user));
36 PetscCall(MatMultTranspose(user->B, X, Y));
37 PetscFunctionReturn(PETSC_SUCCESS);
38 }
39
MatCopy_User(Mat A,Mat X,MatStructure str)40 static PetscErrorCode MatCopy_User(Mat A, Mat X, MatStructure str)
41 {
42 User user, userX;
43
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 }
51
MatDestroy_User(Mat A)52 static PetscErrorCode MatDestroy_User(Mat A)
53 {
54 User user;
55
56 PetscFunctionBegin;
57 PetscCall(MatShellGetContext(A, &user));
58 PetscCall(PetscObjectDereference((PetscObject)user->B));
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
main(int argc,char ** args)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;
78
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));
135
136 if (submat) {
137 Mat A2;
138 IS r, c;
139 PetscInt rst, ren, cst, cen;
140
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 }
151
152 PetscCall(MatGetSize(A, &M, &N));
153 PetscCall(MatGetLocalSize(A, &m, &n));
154 PetscCall(MatHasCongruentLayouts(A, &cong));
155
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"));
160
161 PetscCall(PetscNew(&user));
162 PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, user, &S));
163 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User));
164 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_User));
165 if (cong) PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_User));
166 PetscCall(MatShellSetOperation(S, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_User));
167 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_User));
168 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &user->B));
169
170 /* Square and rows only scaling */
171 ronl = cong ? ronl : PETSC_TRUE;
172
173 for (test = 0; test < ntest; test++) {
174 PetscReal err;
175
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;
187
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));
243
244 if (b1) {
245 Vec b;
246
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));
260
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"));
268
269 PetscCall(MatViewFromOptions(C, NULL, "-view_mat"));
270 PetscCall(MatViewFromOptions(B, NULL, "-view_mat"));
271 PetscCall(MatViewFromOptions(BTT, NULL, "-view_mat"));
272
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));
276
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));
281
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;
290
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 }
304
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 }
309
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 }
314
315 if (testgetdiag && cong) { /* MatGetDiagonal() */
316 Vec dA, dS;
317
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 }
328
329 if (testdup && !test) {
330 Mat A2, S2;
331
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 }
339
340 if (testsubmat) {
341 Mat sA, sS, dA, dS, At, St;
342 IS r, c;
343 PetscInt rst, ren, cst, cen;
344
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 }
386
387 if (testaxpy) {
388 Mat tA, tS, dA, dS;
389 MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN};
390
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 }
417
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 }
426
427 PetscCall(MatDestroy(&A));
428 PetscCall(MatDestroy(&S));
429 PetscCall(PetscFree(user));
430 PetscCall(PetscFinalize());
431 return 0;
432 }
433
434 /*TEST
435
436 testset:
437 suffix: rect
438 requires: !single
439 output_file: output/empty.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}}
442
443 testset:
444 suffix: square
445 requires: !single
446 output_file: output/empty.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*/
450