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=1.e-10; 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 40 #if defined(PETSC_USE_REAL_SINGLE) 41 tol = 1.e-5; 42 #endif 43 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 44 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 45 ierr = VecCreate(PetscObjectComm((PetscObject)A),&x);CHKERRQ(ierr); 46 ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr); 47 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 48 49 ierr = VecCreate(PetscObjectComm((PetscObject)A),&s1);CHKERRQ(ierr); 50 ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr); 51 ierr = VecSetFromOptions(s1);CHKERRQ(ierr); 52 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 53 54 *flg = PETSC_TRUE; 55 for (k=0; k<n; k++) { 56 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 57 ierr = MatMult(A,x,s1);CHKERRQ(ierr); 58 ierr = MatMult(B,x,s2);CHKERRQ(ierr); 59 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 60 if (r2 < tol) { 61 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 62 } else { 63 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 64 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 65 r1 /= r2; 66 } 67 if (r1 > tol) { 68 *flg = PETSC_FALSE; 69 ierr = PetscInfo2(A,"Error: %D-th MatMult() %G\n",k,r1);CHKERRQ(ierr); 70 break; 71 } 72 } 73 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 74 ierr = VecDestroy(&x);CHKERRQ(ierr); 75 ierr = VecDestroy(&s1);CHKERRQ(ierr); 76 ierr = VecDestroy(&s2);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatMultAddEqual" 82 /*@ 83 MatMultAddEqual - Compares matrix-vector products of two matrices. 84 85 Collective on Mat 86 87 Input Parameters: 88 + A - the first matrix 89 - B - the second matrix 90 - n - number of random vectors to be tested 91 92 Output Parameter: 93 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 94 95 Level: intermediate 96 97 Concepts: matrices^equality between 98 @*/ 99 PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 100 { 101 PetscErrorCode ierr; 102 Vec x,y,s1,s2; 103 PetscRandom rctx; 104 PetscReal r1,r2,tol=1.e-10; 105 PetscInt am,an,bm,bn,k; 106 PetscScalar none = -1.0; 107 108 PetscFunctionBegin; 109 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 110 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 111 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); 112 PetscCheckSameComm(A,1,B,2); 113 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 114 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 115 ierr = VecCreate(PetscObjectComm((PetscObject)A),&x);CHKERRQ(ierr); 116 ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr); 117 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 118 119 ierr = VecCreate(PetscObjectComm((PetscObject)A),&s1);CHKERRQ(ierr); 120 ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr); 121 ierr = VecSetFromOptions(s1);CHKERRQ(ierr); 122 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 123 ierr = VecDuplicate(s1,&y);CHKERRQ(ierr); 124 125 *flg = PETSC_TRUE; 126 for (k=0; k<n; k++) { 127 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 128 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr); 129 ierr = MatMultAdd(A,x,y,s1);CHKERRQ(ierr); 130 ierr = MatMultAdd(B,x,y,s2);CHKERRQ(ierr); 131 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 132 if (r2 < tol) { 133 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 134 } else { 135 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 136 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 137 r1 /= r2; 138 } 139 if (r1 > tol) { 140 *flg = PETSC_FALSE; 141 ierr = PetscInfo2(A,"Error: %d-th MatMultAdd() %G\n",k,r1);CHKERRQ(ierr); 142 break; 143 } 144 } 145 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 146 ierr = VecDestroy(&x);CHKERRQ(ierr); 147 ierr = VecDestroy(&y);CHKERRQ(ierr); 148 ierr = VecDestroy(&s1);CHKERRQ(ierr); 149 ierr = VecDestroy(&s2);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "MatMultTransposeEqual" 155 /*@ 156 MatMultTransposeEqual - Compares matrix-vector products of two matrices. 157 158 Collective on Mat 159 160 Input Parameters: 161 + A - the first matrix 162 - B - the second matrix 163 - n - number of random vectors to be tested 164 165 Output Parameter: 166 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 167 168 Level: intermediate 169 170 Concepts: matrices^equality between 171 @*/ 172 PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 173 { 174 PetscErrorCode ierr; 175 Vec x,s1,s2; 176 PetscRandom rctx; 177 PetscReal r1,r2,tol=1.e-10; 178 PetscInt am,an,bm,bn,k; 179 PetscScalar none = -1.0; 180 181 PetscFunctionBegin; 182 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 183 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 184 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); 185 PetscCheckSameComm(A,1,B,2); 186 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 187 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 188 ierr = VecCreate(PetscObjectComm((PetscObject)A),&x);CHKERRQ(ierr); 189 ierr = VecSetSizes(x,am,PETSC_DECIDE);CHKERRQ(ierr); 190 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 191 192 ierr = VecCreate(PetscObjectComm((PetscObject)A),&s1);CHKERRQ(ierr); 193 ierr = VecSetSizes(s1,an,PETSC_DECIDE);CHKERRQ(ierr); 194 ierr = VecSetFromOptions(s1);CHKERRQ(ierr); 195 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 196 197 *flg = PETSC_TRUE; 198 for (k=0; k<n; k++) { 199 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 200 ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr); 201 ierr = MatMultTranspose(B,x,s2);CHKERRQ(ierr); 202 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 203 if (r2 < tol) { 204 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 205 } else { 206 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 207 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 208 r1 /= r2; 209 } 210 if (r1 > tol) { 211 *flg = PETSC_FALSE; 212 ierr = PetscInfo2(A,"Error: %d-th MatMultTranspose() %G\n",k,r1);CHKERRQ(ierr); 213 break; 214 } 215 } 216 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 217 ierr = VecDestroy(&x);CHKERRQ(ierr); 218 ierr = VecDestroy(&s1);CHKERRQ(ierr); 219 ierr = VecDestroy(&s2);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "MatMultTransposeAddEqual" 225 /*@ 226 MatMultTransposeAddEqual - Compares matrix-vector products of two matrices. 227 228 Collective on Mat 229 230 Input Parameters: 231 + A - the first matrix 232 - B - the second matrix 233 - n - number of random vectors to be tested 234 235 Output Parameter: 236 . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise. 237 238 Level: intermediate 239 240 Concepts: matrices^equality between 241 @*/ 242 PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg) 243 { 244 PetscErrorCode ierr; 245 Vec x,y,s1,s2; 246 PetscRandom rctx; 247 PetscReal r1,r2,tol=1.e-10; 248 PetscInt am,an,bm,bn,k; 249 PetscScalar none = -1.0; 250 251 PetscFunctionBegin; 252 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 253 ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr); 254 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); 255 PetscCheckSameComm(A,1,B,2); 256 ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr); 257 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 258 ierr = VecCreate(PetscObjectComm((PetscObject)A),&x);CHKERRQ(ierr); 259 ierr = VecSetSizes(x,am,PETSC_DECIDE);CHKERRQ(ierr); 260 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 261 262 ierr = VecCreate(PetscObjectComm((PetscObject)A),&s1);CHKERRQ(ierr); 263 ierr = VecSetSizes(s1,an,PETSC_DECIDE);CHKERRQ(ierr); 264 ierr = VecSetFromOptions(s1);CHKERRQ(ierr); 265 ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr); 266 ierr = VecDuplicate(s1,&y);CHKERRQ(ierr); 267 268 *flg = PETSC_TRUE; 269 for (k=0; k<n; k++) { 270 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 271 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr); 272 ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr); 273 ierr = MatMultTransposeAdd(B,x,y,s2);CHKERRQ(ierr); 274 ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr); 275 if (r2 < tol) { 276 ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr); 277 } else { 278 ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr); 279 ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr); 280 r1 /= r2; 281 } 282 if (r1 > tol) { 283 *flg = PETSC_FALSE; 284 ierr = PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %G\n",k,r1);CHKERRQ(ierr); 285 break; 286 } 287 } 288 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 289 ierr = VecDestroy(&x);CHKERRQ(ierr); 290 ierr = VecDestroy(&y);CHKERRQ(ierr); 291 ierr = VecDestroy(&s1);CHKERRQ(ierr); 292 ierr = VecDestroy(&s2);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295