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