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