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 Concepts: matrices^equality between 20 @*/ 21 PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 22 { 23 PetscErrorCode ierr; 24 Vec x,s1,s2; 25 PetscRandom rctx; 26 PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON; 27 PetscInt am,an,bm,bn,k; 28 PetscScalar none = -1.0; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 32 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 33 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 34 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 35 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); 36 PetscCheckSameComm(A,1,B,2); 37 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 38 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 39 ierr = MatCreateVecs(A,&x,&s1);CHKERRQ(ierr); 40 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 41 42 *flg = PETSC_TRUE; 43 for (k=0; k<n; k++) { 44 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 45 ierr = MatMult(A,x,s1);CHKERRQ(ierr); 46 ierr = MatMult(B,x,s2);CHKERRQ(ierr); 47 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 48 if (r2 < tol) { 49 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 50 } else { 51 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 52 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 53 r1 /= r2; 54 } 55 if (r1 > tol) { 56 *flg = PETSC_FALSE; 57 ierr = PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);CHKERRQ(ierr); 58 break; 59 } 60 } 61 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 62 ierr = VecDestroy(&x);CHKERRQ(ierr); 63 ierr = VecDestroy(&s1);CHKERRQ(ierr); 64 ierr = VecDestroy(&s2);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 /*@ 69 MatMultAddEqual - Compares matrix-vector products of two matrices. 70 71 Collective on Mat 72 73 Input Parameters: 74 + A - the first matrix 75 - B - the second matrix 76 - n - number of random vectors to be tested 77 78 Output Parameter: 79 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 80 81 Level: intermediate 82 83 Concepts: matrices^equality between 84 @*/ 85 PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 86 { 87 PetscErrorCode ierr; 88 Vec x,y,s1,s2; 89 PetscRandom rctx; 90 PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON; 91 PetscInt am,an,bm,bn,k; 92 PetscScalar none = -1.0; 93 94 PetscFunctionBegin; 95 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 96 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 97 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); 98 PetscCheckSameComm(A,1,B,2); 99 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 100 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 101 ierr = MatCreateVecs(A,&x,&s1);CHKERRQ(ierr); 102 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 103 ierr = VecDuplicate(s1,&y);CHKERRQ(ierr); 104 105 *flg = PETSC_TRUE; 106 for (k=0; k<n; k++) { 107 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 108 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr); 109 ierr = MatMultAdd(A,x,y,s1);CHKERRQ(ierr); 110 ierr = MatMultAdd(B,x,y,s2);CHKERRQ(ierr); 111 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 112 if (r2 < tol) { 113 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 114 } else { 115 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 116 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 117 r1 /= r2; 118 } 119 if (r1 > tol) { 120 *flg = PETSC_FALSE; 121 ierr = PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);CHKERRQ(ierr); 122 break; 123 } 124 } 125 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 126 ierr = VecDestroy(&x);CHKERRQ(ierr); 127 ierr = VecDestroy(&y);CHKERRQ(ierr); 128 ierr = VecDestroy(&s1);CHKERRQ(ierr); 129 ierr = VecDestroy(&s2);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 /*@ 134 MatMultTransposeEqual - Compares matrix-vector products of two matrices. 135 136 Collective on Mat 137 138 Input Parameters: 139 + A - the first matrix 140 - B - the second matrix 141 - n - number of random vectors to be tested 142 143 Output Parameter: 144 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 145 146 Level: intermediate 147 148 Concepts: matrices^equality between 149 @*/ 150 PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 151 { 152 PetscErrorCode ierr; 153 Vec x,s1,s2; 154 PetscRandom rctx; 155 PetscReal r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON; 156 PetscInt am,an,bm,bn,k; 157 PetscScalar none = -1.0; 158 159 PetscFunctionBegin; 160 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 161 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 162 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); 163 PetscCheckSameComm(A,1,B,2); 164 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 165 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 166 ierr = MatCreateVecs(A,&s1,&x);CHKERRQ(ierr); 167 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 168 169 *flg = PETSC_TRUE; 170 for (k=0; k<n; k++) { 171 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 172 ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr); 173 ierr = MatMultTranspose(B,x,s2);CHKERRQ(ierr); 174 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 175 if (r2 < tol) { 176 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 177 } else { 178 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 179 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 180 r1 /= r2; 181 } 182 if (r1 > tol) { 183 *flg = PETSC_FALSE; 184 ierr = PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);CHKERRQ(ierr); 185 break; 186 } 187 } 188 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 189 ierr = VecDestroy(&x);CHKERRQ(ierr); 190 ierr = VecDestroy(&s1);CHKERRQ(ierr); 191 ierr = VecDestroy(&s2);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 /*@ 196 MatMultTransposeAddEqual - Compares matrix-vector products of two matrices. 197 198 Collective on Mat 199 200 Input Parameters: 201 + A - the first matrix 202 - B - the second matrix 203 - n - number of random vectors to be tested 204 205 Output Parameter: 206 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 207 208 Level: intermediate 209 210 Concepts: matrices^equality between 211 @*/ 212 PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 213 { 214 PetscErrorCode ierr; 215 Vec x,y,s1,s2; 216 PetscRandom rctx; 217 PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON; 218 PetscInt am,an,bm,bn,k; 219 PetscScalar none = -1.0; 220 221 PetscFunctionBegin; 222 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 223 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 224 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); 225 PetscCheckSameComm(A,1,B,2); 226 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 227 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 228 ierr = MatCreateVecs(A,&s1,&x);CHKERRQ(ierr); 229 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 230 ierr = VecDuplicate(s1,&y);CHKERRQ(ierr); 231 232 *flg = PETSC_TRUE; 233 for (k=0; k<n; k++) { 234 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 235 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr); 236 ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr); 237 ierr = MatMultTransposeAdd(B,x,y,s2);CHKERRQ(ierr); 238 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 239 if (r2 < tol) { 240 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 241 } else { 242 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 243 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 244 r1 /= r2; 245 } 246 if (r1 > tol) { 247 *flg = PETSC_FALSE; 248 ierr = PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);CHKERRQ(ierr); 249 break; 250 } 251 } 252 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 253 ierr = VecDestroy(&x);CHKERRQ(ierr); 254 ierr = VecDestroy(&y);CHKERRQ(ierr); 255 ierr = VecDestroy(&s1);CHKERRQ(ierr); 256 ierr = VecDestroy(&s2);CHKERRQ(ierr); 257 PetscFunctionReturn(0); 258 } 259 260 /*@ 261 MatMatMultEqual - Test A*B*x = C*x for n random vector x 262 263 Collective on Mat 264 265 Input Parameters: 266 + A - the first matrix 267 - B - the second matrix 268 - C - the third matrix 269 - n - number of random vectors to be tested 270 271 Output Parameter: 272 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 273 274 Level: intermediate 275 276 Concepts: matrices^equality between 277 @*/ 278 PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg) 279 { 280 PetscErrorCode ierr; 281 Vec x,s1,s2,s3; 282 PetscRandom rctx; 283 PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON; 284 PetscInt am,an,bm,bn,cm,cn,k; 285 PetscScalar none = -1.0; 286 287 PetscFunctionBegin; 288 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 289 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 290 ierr = MatGetLocalSize(C,&cm,&cn);CHKERRQ(ierr); 291 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); 292 293 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);CHKERRQ(ierr); 294 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 295 ierr = MatCreateVecs(B,&x,&s1);CHKERRQ(ierr); 296 ierr = MatCreateVecs(C,NULL,&s2);CHKERRQ(ierr); 297 ierr = VecDuplicate(s2,&s3);CHKERRQ(ierr); 298 299 *flg = PETSC_TRUE; 300 for (k=0; k<n; k++) { 301 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 302 ierr = MatMult(B,x,s1);CHKERRQ(ierr); 303 ierr = MatMult(A,s1,s2);CHKERRQ(ierr); 304 ierr = MatMult(C,x,s3);CHKERRQ(ierr); 305 306 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 307 if (r2 < tol) { 308 ierr = VecNorm(s3,NORM_INFINITY,&r1);CHKERRQ(ierr); 309 } else { 310 ierr = VecAXPY(s2,none,s3);CHKERRQ(ierr); 311 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 312 r1 /= r2; 313 } 314 if (r1 > tol) { 315 *flg = PETSC_FALSE; 316 ierr = PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);CHKERRQ(ierr); 317 break; 318 } 319 } 320 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 321 ierr = VecDestroy(&x);CHKERRQ(ierr); 322 ierr = VecDestroy(&s1);CHKERRQ(ierr); 323 ierr = VecDestroy(&s2);CHKERRQ(ierr); 324 ierr = VecDestroy(&s3);CHKERRQ(ierr); 325 PetscFunctionReturn(0); 326 } 327 328 /*@ 329 MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x 330 331 Collective on Mat 332 333 Input Parameters: 334 + A - the first matrix 335 - B - the second matrix 336 - C - the third matrix 337 - n - number of random vectors to be tested 338 339 Output Parameter: 340 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 341 342 Level: intermediate 343 344 Concepts: matrices^equality between 345 @*/ 346 PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg) 347 { 348 PetscErrorCode ierr; 349 Vec x,s1,s2,s3; 350 PetscRandom rctx; 351 PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON; 352 PetscInt am,an,bm,bn,cm,cn,k; 353 PetscScalar none = -1.0; 354 355 PetscFunctionBegin; 356 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 357 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 358 ierr = MatGetLocalSize(C,&cm,&cn);CHKERRQ(ierr); 359 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); 360 361 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);CHKERRQ(ierr); 362 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 363 ierr = MatCreateVecs(B,&x,&s1);CHKERRQ(ierr); 364 ierr = MatCreateVecs(C,NULL,&s2);CHKERRQ(ierr); 365 ierr = VecDuplicate(s2,&s3);CHKERRQ(ierr); 366 367 *flg = PETSC_TRUE; 368 for (k=0; k<n; k++) { 369 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 370 ierr = MatMult(B,x,s1);CHKERRQ(ierr); 371 ierr = MatMultTranspose(A,s1,s2);CHKERRQ(ierr); 372 ierr = MatMult(C,x,s3);CHKERRQ(ierr); 373 374 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 375 if (r2 < tol) { 376 ierr = VecNorm(s3,NORM_INFINITY,&r1);CHKERRQ(ierr); 377 } else { 378 ierr = VecAXPY(s2,none,s3);CHKERRQ(ierr); 379 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 380 r1 /= r2; 381 } 382 if (r1 > tol) { 383 *flg = PETSC_FALSE; 384 ierr = PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);CHKERRQ(ierr); 385 break; 386 } 387 } 388 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 389 ierr = VecDestroy(&x);CHKERRQ(ierr); 390 ierr = VecDestroy(&s1);CHKERRQ(ierr); 391 ierr = VecDestroy(&s2);CHKERRQ(ierr); 392 ierr = VecDestroy(&s3);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395