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