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