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