1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
3 /*
4 n; try the MatMult variant n times
5 flg: return the boolean result, equal or not
6 t: 0 => no transpose; 1 => transpose; 2 => Hermitian transpose
7 add: 0 => no add (e.g., y = Ax); 1 => add third vector (e.g., z = Ax + y); 2 => add update (e.g., y = Ax + y)
8 */
MatMultEqual_Private(Mat A,Mat B,PetscInt n,PetscBool * flg,PetscInt t,PetscInt add)9 static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscInt add)
10 {
11 Vec Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
12 PetscRandom rctx;
13 PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
14 PetscInt am, an, bm, bn, k;
15 PetscScalar none = -1.0;
16 #if defined(PETSC_USE_INFO)
17 const char *sops[] = {"MatMult", "MatMultAdd", "MatMultAdd (update)", "MatMultTranspose", "MatMultTransposeAdd", "MatMultTransposeAdd (update)", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd", "MatMultHermitianTransposeAdd (update)"};
18 const char *sop;
19 #endif
20
21 PetscFunctionBegin;
22 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
23 PetscValidHeaderSpecific(B, MAT_CLASSID, 2);
24 PetscCheckSameComm(A, 1, B, 2);
25 PetscValidLogicalCollectiveInt(A, n, 3);
26 PetscAssertPointer(flg, 4);
27 PetscValidLogicalCollectiveInt(A, t, 5);
28 PetscValidLogicalCollectiveInt(A, add, 6);
29 PetscCall(MatGetLocalSize(A, &am, &an));
30 PetscCall(MatGetLocalSize(B, &bm, &bn));
31 PetscCheck(am == bm && an == bn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A,Mat B: local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, bm, an, bn);
32 #if defined(PETSC_USE_INFO)
33 sop = sops[add + 3 * t];
34 #endif
35 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
36 PetscCall(PetscRandomSetFromOptions(rctx));
37 if (t) {
38 PetscCall(MatCreateVecs(A, &s1, &Ax));
39 PetscCall(MatCreateVecs(B, &s2, &Bx));
40 } else {
41 PetscCall(MatCreateVecs(A, &Ax, &s1));
42 PetscCall(MatCreateVecs(B, &Bx, &s2));
43 }
44 if (add) {
45 PetscCall(VecDuplicate(s1, &Ay));
46 PetscCall(VecDuplicate(s2, &By));
47 }
48
49 *flg = PETSC_TRUE;
50 for (k = 0; k < n; k++) {
51 Vec Aadd = NULL, Badd = NULL;
52
53 PetscCall(VecSetRandom(Ax, rctx));
54 PetscCall(VecCopy(Ax, Bx));
55 if (add) {
56 PetscCall(VecSetRandom(Ay, rctx));
57 PetscCall(VecCopy(Ay, By));
58 Aadd = Ay;
59 Badd = By;
60 if (add == 2) {
61 PetscCall(VecCopy(Ay, s1));
62 PetscCall(VecCopy(By, s2));
63 Aadd = s1;
64 Badd = s2;
65 }
66 }
67 if (t == 1) {
68 if (add) {
69 PetscCall(MatMultTransposeAdd(A, Ax, Aadd, s1));
70 PetscCall(MatMultTransposeAdd(B, Bx, Badd, s2));
71 } else {
72 PetscCall(MatMultTranspose(A, Ax, s1));
73 PetscCall(MatMultTranspose(B, Bx, s2));
74 }
75 } else if (t == 2) {
76 if (add) {
77 PetscCall(MatMultHermitianTransposeAdd(A, Ax, Aadd, s1));
78 PetscCall(MatMultHermitianTransposeAdd(B, Bx, Badd, s2));
79 } else {
80 PetscCall(MatMultHermitianTranspose(A, Ax, s1));
81 PetscCall(MatMultHermitianTranspose(B, Bx, s2));
82 }
83 } else {
84 if (add) {
85 PetscCall(MatMultAdd(A, Ax, Aadd, s1));
86 PetscCall(MatMultAdd(B, Bx, Badd, s2));
87 } else {
88 PetscCall(MatMult(A, Ax, s1));
89 PetscCall(MatMult(B, Bx, s2));
90 }
91 }
92 PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
93 if (r2 < tol) {
94 PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
95 } else {
96 PetscCall(VecAXPY(s2, none, s1));
97 PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
98 r1 /= r2;
99 }
100 if (r1 > tol) {
101 *flg = PETSC_FALSE;
102 PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
103 break;
104 }
105 }
106 PetscCall(PetscRandomDestroy(&rctx));
107 PetscCall(VecDestroy(&Ax));
108 PetscCall(VecDestroy(&Bx));
109 PetscCall(VecDestroy(&Ay));
110 PetscCall(VecDestroy(&By));
111 PetscCall(VecDestroy(&s1));
112 PetscCall(VecDestroy(&s2));
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
MatMatMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg,PetscBool At,PetscBool Bt)116 static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
117 {
118 Vec Ax, Bx, Cx, s1, s2, s3;
119 PetscRandom rctx;
120 PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
121 PetscInt am, an, bm, bn, cm, cn, k;
122 PetscScalar none = -1.0;
123 #if defined(PETSC_USE_INFO)
124 const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
125 const char *sop;
126 #endif
127
128 PetscFunctionBegin;
129 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
130 PetscValidHeaderSpecific(B, MAT_CLASSID, 2);
131 PetscCheckSameComm(A, 1, B, 2);
132 PetscValidHeaderSpecific(C, MAT_CLASSID, 3);
133 PetscCheckSameComm(A, 1, C, 3);
134 PetscValidLogicalCollectiveInt(A, n, 4);
135 PetscAssertPointer(flg, 5);
136 PetscValidLogicalCollectiveBool(A, At, 6);
137 PetscValidLogicalCollectiveBool(B, Bt, 7);
138 PetscCall(MatGetLocalSize(A, &am, &an));
139 PetscCall(MatGetLocalSize(B, &bm, &bn));
140 PetscCall(MatGetLocalSize(C, &cm, &cn));
141 if (At) {
142 PetscInt tt = an;
143 an = am;
144 am = tt;
145 }
146 if (Bt) {
147 PetscInt tt = bn;
148 bn = bm;
149 bm = tt;
150 }
151 PetscCheck(an == bm && am == cm && bn == cn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, an, bm, bn, cm, cn);
152
153 #if defined(PETSC_USE_INFO)
154 sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
155 #endif
156 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
157 PetscCall(PetscRandomSetFromOptions(rctx));
158 if (Bt) {
159 PetscCall(MatCreateVecs(B, &s1, &Bx));
160 } else {
161 PetscCall(MatCreateVecs(B, &Bx, &s1));
162 }
163 if (At) {
164 PetscCall(MatCreateVecs(A, &s2, &Ax));
165 } else {
166 PetscCall(MatCreateVecs(A, &Ax, &s2));
167 }
168 PetscCall(MatCreateVecs(C, &Cx, &s3));
169
170 *flg = PETSC_TRUE;
171 for (k = 0; k < n; k++) {
172 PetscCall(VecSetRandom(Bx, rctx));
173 if (Bt) {
174 PetscCall(MatMultTranspose(B, Bx, s1));
175 } else {
176 PetscCall(MatMult(B, Bx, s1));
177 }
178 PetscCall(VecCopy(s1, Ax));
179 if (At) {
180 PetscCall(MatMultTranspose(A, Ax, s2));
181 } else {
182 PetscCall(MatMult(A, Ax, s2));
183 }
184 PetscCall(VecCopy(Bx, Cx));
185 PetscCall(MatMult(C, Cx, s3));
186
187 PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
188 if (r2 < tol) {
189 PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
190 } else {
191 PetscCall(VecAXPY(s2, none, s3));
192 PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
193 r1 /= r2;
194 }
195 if (r1 > tol) {
196 *flg = PETSC_FALSE;
197 PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
198 break;
199 }
200 }
201 PetscCall(PetscRandomDestroy(&rctx));
202 PetscCall(VecDestroy(&Ax));
203 PetscCall(VecDestroy(&Bx));
204 PetscCall(VecDestroy(&Cx));
205 PetscCall(VecDestroy(&s1));
206 PetscCall(VecDestroy(&s2));
207 PetscCall(VecDestroy(&s3));
208 PetscFunctionReturn(PETSC_SUCCESS);
209 }
210
211 /*@
212 MatMultEqual - Compares matrix-vector products of two matrices using `n` random vectors
213
214 Collective
215
216 Input Parameters:
217 + A - the first matrix
218 . B - the second matrix
219 - n - number of random vectors to be tested
220
221 Output Parameter:
222 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
223
224 Level: intermediate
225
226 Note:
227 The tolerance for equality is a generous `PETSC_SQRT_MACHINE_EPSILON` in the norm of the difference of the two computed vectors to
228 allow for differences in the numerical computations. Hence this routine may indicate equality even if there is a small systematic difference
229 between the two matrices.
230
231 .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`, `MatEqual()`
232 @*/
MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)233 PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
234 {
235 PetscFunctionBegin;
236 PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
237 PetscFunctionReturn(PETSC_SUCCESS);
238 }
239
240 /*@
241 MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.
242
243 Collective
244
245 Input Parameters:
246 + A - the first matrix
247 . B - the second matrix
248 - n - number of random vectors to be tested
249
250 Output Parameter:
251 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
252
253 Level: intermediate
254
255 .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
256 @*/
MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)257 PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
258 {
259 PetscFunctionBegin;
260 PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
261 PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
265 /*@
266 MatMultTransposeEqual - Compares matrix-vector products of two matrices.
267
268 Collective
269
270 Input Parameters:
271 + A - the first matrix
272 . B - the second matrix
273 - n - number of random vectors to be tested
274
275 Output Parameter:
276 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
277
278 Level: intermediate
279
280 .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
281 @*/
MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)282 PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
283 {
284 PetscFunctionBegin;
285 PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 0));
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
289 /*@
290 MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
291
292 Collective
293
294 Input Parameters:
295 + A - the first matrix
296 . B - the second matrix
297 - n - number of random vectors to be tested
298
299 Output Parameter:
300 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
301
302 Level: intermediate
303
304 .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
305 @*/
MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)306 PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
307 {
308 PetscFunctionBegin;
309 PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
310 PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
311 PetscFunctionReturn(PETSC_SUCCESS);
312 }
313
314 /*@
315 MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
316
317 Collective
318
319 Input Parameters:
320 + A - the first matrix
321 . B - the second matrix
322 - n - number of random vectors to be tested
323
324 Output Parameter:
325 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
326
327 Level: intermediate
328
329 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
330 @*/
MatMultHermitianTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)331 PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
332 {
333 PetscFunctionBegin;
334 PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
335 PetscFunctionReturn(PETSC_SUCCESS);
336 }
337
338 /*@
339 MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
340
341 Collective
342
343 Input Parameters:
344 + A - the first matrix
345 . B - the second matrix
346 - n - number of random vectors to be tested
347
348 Output Parameter:
349 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
350
351 Level: intermediate
352
353 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
354 @*/
MatMultHermitianTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool * flg)355 PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
356 {
357 PetscFunctionBegin;
358 PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
359 PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
360 PetscFunctionReturn(PETSC_SUCCESS);
361 }
362
363 /*@
364 MatMatMultEqual - Test A*B*x = C*x for n random vector x
365
366 Collective
367
368 Input Parameters:
369 + A - the first matrix
370 . B - the second matrix
371 . C - the third matrix
372 - n - number of random vectors to be tested
373
374 Output Parameter:
375 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
376
377 Level: intermediate
378
379 .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
380 @*/
MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg)381 PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
382 {
383 PetscFunctionBegin;
384 PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
385 PetscFunctionReturn(PETSC_SUCCESS);
386 }
387
388 /*@
389 MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
390
391 Collective
392
393 Input Parameters:
394 + A - the first matrix
395 . B - the second matrix
396 . C - the third matrix
397 - n - number of random vectors to be tested
398
399 Output Parameter:
400 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
401
402 Level: intermediate
403
404 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
405 @*/
MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg)406 PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
407 {
408 PetscFunctionBegin;
409 PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
410 PetscFunctionReturn(PETSC_SUCCESS);
411 }
412
413 /*@
414 MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
415
416 Collective
417
418 Input Parameters:
419 + A - the first matrix
420 . B - the second matrix
421 . C - the third matrix
422 - n - number of random vectors to be tested
423
424 Output Parameter:
425 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
426
427 Level: intermediate
428
429 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
430 @*/
MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg)431 PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
432 {
433 PetscFunctionBegin;
434 PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
435 PetscFunctionReturn(PETSC_SUCCESS);
436 }
437
MatProjMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool rart,PetscBool * flg)438 static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
439 {
440 Vec x, v1, v2, v3, v4, Cx, Bx;
441 PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
442 PetscInt i, am, an, bm, bn, cm, cn;
443 PetscRandom rdm;
444
445 PetscFunctionBegin;
446 PetscCall(MatGetLocalSize(A, &am, &an));
447 PetscCall(MatGetLocalSize(B, &bm, &bn));
448 if (rart) {
449 PetscInt t = bm;
450 bm = bn;
451 bn = t;
452 }
453 PetscCall(MatGetLocalSize(C, &cm, &cn));
454 PetscCheck(an == bm && bn == cm && bn == cn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, an, bm, bn, cm, cn);
455
456 /* Create left vector of A: v2 */
457 PetscCall(MatCreateVecs(A, &Bx, &v2));
458
459 /* Create right vectors of B: x, v3, v4 */
460 if (rart) {
461 PetscCall(MatCreateVecs(B, &v1, &x));
462 } else {
463 PetscCall(MatCreateVecs(B, &x, &v1));
464 }
465 PetscCall(VecDuplicate(x, &v3));
466
467 PetscCall(MatCreateVecs(C, &Cx, &v4));
468 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
469 PetscCall(PetscRandomSetFromOptions(rdm));
470
471 *flg = PETSC_TRUE;
472 for (i = 0; i < n; i++) {
473 PetscCall(VecSetRandom(x, rdm));
474 PetscCall(VecCopy(x, Cx));
475 PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
476 if (rart) {
477 PetscCall(MatMultTranspose(B, x, v1));
478 } else {
479 PetscCall(MatMult(B, x, v1));
480 }
481 PetscCall(VecCopy(v1, Bx));
482 PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
483 PetscCall(VecCopy(v2, v1));
484 if (rart) {
485 PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
486 } else {
487 PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
488 }
489 PetscCall(VecNorm(v4, NORM_2, &norm_abs));
490 PetscCall(VecAXPY(v4, -1.0, v3));
491 PetscCall(VecNorm(v4, NORM_2, &norm_rel));
492
493 if (norm_abs > tol) norm_rel /= norm_abs;
494 if (norm_rel > tol || PetscIsInfOrNanReal(norm_rel)) {
495 *flg = PETSC_FALSE;
496 PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
497 break;
498 }
499 }
500
501 PetscCall(PetscRandomDestroy(&rdm));
502 PetscCall(VecDestroy(&x));
503 PetscCall(VecDestroy(&Bx));
504 PetscCall(VecDestroy(&Cx));
505 PetscCall(VecDestroy(&v1));
506 PetscCall(VecDestroy(&v2));
507 PetscCall(VecDestroy(&v3));
508 PetscCall(VecDestroy(&v4));
509 PetscFunctionReturn(PETSC_SUCCESS);
510 }
511
512 /*@
513 MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
514
515 Collective
516
517 Input Parameters:
518 + A - the first matrix
519 . B - the second matrix
520 . C - the third matrix
521 - n - number of random vectors to be tested
522
523 Output Parameter:
524 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
525
526 Level: intermediate
527
528 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
529 @*/
MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg)530 PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
531 {
532 PetscFunctionBegin;
533 PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
534 PetscFunctionReturn(PETSC_SUCCESS);
535 }
536
537 /*@
538 MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
539
540 Collective
541
542 Input Parameters:
543 + A - the first matrix
544 . B - the second matrix
545 . C - the third matrix
546 - n - number of random vectors to be tested
547
548 Output Parameter:
549 . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
550
551 Level: intermediate
552
553 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
554 @*/
MatRARtMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool * flg)555 PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
556 {
557 PetscFunctionBegin;
558 PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
559 PetscFunctionReturn(PETSC_SUCCESS);
560 }
561
562 /*@
563 MatIsLinear - Check if a shell matrix `A` is a linear operator.
564
565 Collective
566
567 Input Parameters:
568 + A - the shell matrix
569 - n - number of random vectors to be tested
570
571 Output Parameter:
572 . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
573
574 Level: intermediate
575
576 .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
577 @*/
MatIsLinear(Mat A,PetscInt n,PetscBool * flg)578 PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
579 {
580 Vec x, y, s1, s2;
581 PetscRandom rctx;
582 PetscScalar a;
583 PetscInt k;
584 PetscReal norm, normA;
585 MPI_Comm comm;
586 PetscMPIInt rank;
587
588 PetscFunctionBegin;
589 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
590 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
591 PetscCallMPI(MPI_Comm_rank(comm, &rank));
592
593 PetscCall(PetscRandomCreate(comm, &rctx));
594 PetscCall(PetscRandomSetFromOptions(rctx));
595 PetscCall(MatCreateVecs(A, &x, &s1));
596 PetscCall(VecDuplicate(x, &y));
597 PetscCall(VecDuplicate(s1, &s2));
598
599 *flg = PETSC_TRUE;
600 for (k = 0; k < n; k++) {
601 PetscCall(VecSetRandom(x, rctx));
602 PetscCall(VecSetRandom(y, rctx));
603 if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
604 PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
605
606 /* s2 = a*A*x + A*y */
607 PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
608 PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
609 PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
610
611 /* s1 = A * (a x + y) */
612 PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
613 PetscCall(MatMult(A, y, s1));
614 PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
615
616 PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
617 PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
618 if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
619 *flg = PETSC_FALSE;
620 PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th |A*(ax+y) - (a*A*x+A*y)|/|A(ax+y)| %g > tol %g\n", k, (double)(norm / normA), (double)(100 * PETSC_MACHINE_EPSILON)));
621 break;
622 }
623 }
624 PetscCall(PetscRandomDestroy(&rctx));
625 PetscCall(VecDestroy(&x));
626 PetscCall(VecDestroy(&y));
627 PetscCall(VecDestroy(&s1));
628 PetscCall(VecDestroy(&s2));
629 PetscFunctionReturn(PETSC_SUCCESS);
630 }
631