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