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