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