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