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