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