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