1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 static PetscErrorCode MatMultEqual_Private(Mat A,Mat B,PetscInt n,PetscBool *flg,PetscBool t,PetscBool add) 5 { 6 PetscErrorCode ierr; 7 Vec Ax = NULL,Bx = NULL,s1 = NULL,s2 = NULL,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 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 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 @*/ 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 @*/ 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 @*/ 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);CHKERRMPI(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 == 0) { 517 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 518 } 519 ierr = MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);CHKERRMPI(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