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