1 2 static char help[] = "Tests various routines in MatSeqBAIJ format.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,B,C,D,Fact; 9 Vec xx,s1,s2,yy; 10 PetscInt m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M; 11 PetscScalar rval,vals1[4],vals2[4]; 12 PetscRandom rdm; 13 IS is1,is2; 14 PetscReal s1norm,s2norm,rnorm,tol = 1.e-4; 15 PetscBool flg; 16 MatFactorInfo info; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 20 /* Test MatSetValues() and MatGetValues() */ 21 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL)); 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL)); 23 M = m*bs; 24 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A)); 25 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 26 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B)); 27 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 28 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 29 PetscCall(PetscRandomSetFromOptions(rdm)); 30 PetscCall(VecCreateSeq(PETSC_COMM_SELF,M,&xx)); 31 PetscCall(VecDuplicate(xx,&s1)); 32 PetscCall(VecDuplicate(xx,&s2)); 33 PetscCall(VecDuplicate(xx,&yy)); 34 35 /* For each row add at least 15 elements */ 36 for (row=0; row<M; row++) { 37 for (i=0; i<25*bs; i++) { 38 PetscCall(PetscRandomGetValue(rdm,&rval)); 39 col = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 40 PetscCall(MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES)); 41 PetscCall(MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES)); 42 } 43 } 44 45 /* Now set blocks of values */ 46 for (i=0; i<20*bs; i++) { 47 PetscCall(PetscRandomGetValue(rdm,&rval)); 48 cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 49 vals1[0] = rval; 50 PetscCall(PetscRandomGetValue(rdm,&rval)); 51 cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 52 vals1[1] = rval; 53 PetscCall(PetscRandomGetValue(rdm,&rval)); 54 rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 55 vals1[2] = rval; 56 PetscCall(PetscRandomGetValue(rdm,&rval)); 57 rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 58 vals1[3] = rval; 59 PetscCall(MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES)); 60 PetscCall(MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES)); 61 } 62 63 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 64 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 65 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 66 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 67 68 /* Test MatNorm() */ 69 PetscCall(MatNorm(A,NORM_FROBENIUS,&s1norm)); 70 PetscCall(MatNorm(B,NORM_FROBENIUS,&s2norm)); 71 rnorm = PetscAbsReal(s2norm-s1norm)/s2norm; 72 if (rnorm>tol) { 73 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 74 } 75 PetscCall(MatNorm(A,NORM_INFINITY,&s1norm)); 76 PetscCall(MatNorm(B,NORM_INFINITY,&s2norm)); 77 rnorm = PetscAbsReal(s2norm-s1norm)/s2norm; 78 if (rnorm>tol) { 79 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 80 } 81 PetscCall(MatNorm(A,NORM_1,&s1norm)); 82 PetscCall(MatNorm(B,NORM_1,&s2norm)); 83 rnorm = PetscAbsReal(s2norm-s1norm)/s2norm; 84 if (rnorm>tol) { 85 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 86 } 87 88 /* MatShift() */ 89 rval = 10*s1norm; 90 PetscCall(MatShift(A,rval)); 91 PetscCall(MatShift(B,rval)); 92 93 /* Test MatTranspose() */ 94 PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A)); 95 PetscCall(MatTranspose(B,MAT_INPLACE_MATRIX,&B)); 96 97 /* Now do MatGetValues() */ 98 for (i=0; i<30; i++) { 99 PetscCall(PetscRandomGetValue(rdm,&rval)); 100 cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 101 PetscCall(PetscRandomGetValue(rdm,&rval)); 102 cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 103 PetscCall(PetscRandomGetValue(rdm,&rval)); 104 rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 105 PetscCall(PetscRandomGetValue(rdm,&rval)); 106 rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M)); 107 PetscCall(MatGetValues(A,2,rows,2,cols,vals1)); 108 PetscCall(MatGetValues(B,2,rows,2,cols,vals2)); 109 PetscCall(PetscArraycmp(vals1,vals2,4,&flg)); 110 if (!flg) { 111 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %" PetscInt_FMT "\n",bs)); 112 } 113 } 114 115 /* Test MatMult(), MatMultAdd() */ 116 for (i=0; i<40; i++) { 117 PetscCall(VecSetRandom(xx,rdm)); 118 PetscCall(VecSet(s2,0.0)); 119 PetscCall(MatMult(A,xx,s1)); 120 PetscCall(MatMultAdd(A,xx,s2,s2)); 121 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 122 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 123 rnorm = s2norm-s1norm; 124 if (rnorm<-tol || rnorm>tol) { 125 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 126 } 127 } 128 129 /* Test MatMult() */ 130 PetscCall(MatMultEqual(A,B,10,&flg)); 131 if (!flg) { 132 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n")); 133 } 134 135 /* Test MatMultAdd() */ 136 PetscCall(MatMultAddEqual(A,B,10,&flg)); 137 if (!flg) { 138 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n")); 139 } 140 141 /* Test MatMultTranspose() */ 142 PetscCall(MatMultTransposeEqual(A,B,10,&flg)); 143 if (!flg) { 144 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n")); 145 } 146 147 /* Test MatMultTransposeAdd() */ 148 PetscCall(MatMultTransposeAddEqual(A,B,10,&flg)); 149 if (!flg) { 150 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n")); 151 } 152 153 /* Test MatMatMult() */ 154 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&C)); 155 PetscCall(MatSetRandom(C,rdm)); 156 PetscCall(MatMatMult(A,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D)); 157 PetscCall(MatMatMultEqual(A,C,D,40,&flg)); 158 if (!flg) { 159 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n")); 160 } 161 PetscCall(MatDestroy(&D)); 162 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&D)); 163 PetscCall(MatMatMult(A,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&D)); 164 PetscCall(MatMatMultEqual(A,C,D,40,&flg)); 165 if (!flg) { 166 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n")); 167 } 168 169 /* Do LUFactor() on both the matrices */ 170 PetscCall(PetscMalloc1(M,&idx)); 171 for (i=0; i<M; i++) idx[i] = i; 172 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1)); 173 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2)); 174 PetscCall(PetscFree(idx)); 175 PetscCall(ISSetPermutation(is1)); 176 PetscCall(ISSetPermutation(is2)); 177 178 PetscCall(MatFactorInfoInitialize(&info)); 179 180 info.fill = 2.0; 181 info.dtcol = 0.0; 182 info.zeropivot = 1.e-14; 183 info.pivotinblocks = 1.0; 184 185 if (bs < 4) { 186 PetscCall(MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact)); 187 PetscCall(MatLUFactorSymbolic(Fact,A,is1,is2,&info)); 188 PetscCall(MatLUFactorNumeric(Fact,A,&info)); 189 PetscCall(VecSetRandom(yy,rdm)); 190 PetscCall(MatForwardSolve(Fact,yy,xx)); 191 PetscCall(MatBackwardSolve(Fact,xx,s1)); 192 PetscCall(MatDestroy(&Fact)); 193 PetscCall(VecScale(s1,-1.0)); 194 PetscCall(MatMultAdd(A,s1,yy,yy)); 195 PetscCall(VecNorm(yy,NORM_2,&rnorm)); 196 if (rnorm<-tol || rnorm>tol) { 197 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n",(double)rnorm,bs)); 198 } 199 } 200 201 PetscCall(MatLUFactor(B,is1,is2,&info)); 202 PetscCall(MatLUFactor(A,is1,is2,&info)); 203 204 /* Test MatSolveAdd() */ 205 for (i=0; i<10; i++) { 206 PetscCall(VecSetRandom(xx,rdm)); 207 PetscCall(VecSetRandom(yy,rdm)); 208 PetscCall(MatSolveAdd(B,xx,yy,s2)); 209 PetscCall(MatSolveAdd(A,xx,yy,s1)); 210 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 211 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 212 rnorm = s2norm-s1norm; 213 if (rnorm<-tol || rnorm>tol) { 214 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 215 } 216 } 217 218 /* Test MatSolveAdd() when x = A'b +x */ 219 for (i=0; i<10; i++) { 220 PetscCall(VecSetRandom(xx,rdm)); 221 PetscCall(VecSetRandom(s1,rdm)); 222 PetscCall(VecCopy(s2,s1)); 223 PetscCall(MatSolveAdd(B,xx,s2,s2)); 224 PetscCall(MatSolveAdd(A,xx,s1,s1)); 225 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 226 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 227 rnorm = s2norm-s1norm; 228 if (rnorm<-tol || rnorm>tol) { 229 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 230 } 231 } 232 233 /* Test MatSolve() */ 234 for (i=0; i<10; i++) { 235 PetscCall(VecSetRandom(xx,rdm)); 236 PetscCall(MatSolve(B,xx,s2)); 237 PetscCall(MatSolve(A,xx,s1)); 238 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 239 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 240 rnorm = s2norm-s1norm; 241 if (rnorm<-tol || rnorm>tol) { 242 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 243 } 244 } 245 246 /* Test MatSolveTranspose() */ 247 if (bs < 8) { 248 for (i=0; i<10; i++) { 249 PetscCall(VecSetRandom(xx,rdm)); 250 PetscCall(MatSolveTranspose(B,xx,s2)); 251 PetscCall(MatSolveTranspose(A,xx,s1)); 252 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 253 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 254 rnorm = s2norm-s1norm; 255 if (rnorm<-tol || rnorm>tol) { 256 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs)); 257 } 258 } 259 } 260 261 PetscCall(MatDestroy(&A)); 262 PetscCall(MatDestroy(&B)); 263 PetscCall(MatDestroy(&C)); 264 PetscCall(MatDestroy(&D)); 265 PetscCall(VecDestroy(&xx)); 266 PetscCall(VecDestroy(&s1)); 267 PetscCall(VecDestroy(&s2)); 268 PetscCall(VecDestroy(&yy)); 269 PetscCall(ISDestroy(&is1)); 270 PetscCall(ISDestroy(&is2)); 271 PetscCall(PetscRandomDestroy(&rdm)); 272 PetscCall(PetscFinalize()); 273 return 0; 274 } 275 276 /*TEST 277 278 test: 279 args: -mat_block_size {{1 2 3 4 5 6 7 8}} 280 281 TEST*/ 282