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