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