1 2 static char help[] = "Tests various routines in MatMPIBAIJ format.\n"; 3 4 #include <petscmat.h> 5 #define IMAX 15 6 int main(int argc,char **args) 7 { 8 Mat A,B,C,At,Bt; 9 PetscViewer fd; 10 char file[PETSC_MAX_PATH_LEN]; 11 PetscRandom rand; 12 Vec xx,yy,s1,s2; 13 PetscReal s1norm,s2norm,rnorm,tol=1.e-10; 14 PetscInt rstart,rend,rows[2],cols[2],m,n,i,j,M,N,ct,row,ncols1,ncols2,bs; 15 PetscMPIInt rank,size; 16 const PetscInt *cols1,*cols2; 17 PetscScalar vals1[4],vals2[4],v; 18 const PetscScalar *v1,*v2; 19 PetscBool flg; 20 21 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 22 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 23 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 25 /* Check out if MatLoad() works */ 26 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg)); 27 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Input file not specified"); 28 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 29 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 30 CHKERRQ(MatSetType(A,MATBAIJ)); 31 CHKERRQ(MatLoad(A,fd)); 32 33 CHKERRQ(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B)); 34 CHKERRQ(PetscViewerDestroy(&fd)); 35 36 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 37 CHKERRQ(PetscRandomSetFromOptions(rand)); 38 CHKERRQ(MatGetLocalSize(A,&m,&n)); 39 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&xx)); 40 CHKERRQ(VecSetSizes(xx,m,PETSC_DECIDE)); 41 CHKERRQ(VecSetFromOptions(xx)); 42 CHKERRQ(VecDuplicate(xx,&s1)); 43 CHKERRQ(VecDuplicate(xx,&s2)); 44 CHKERRQ(VecDuplicate(xx,&yy)); 45 CHKERRQ(MatGetBlockSize(A,&bs)); 46 47 /* Test MatNorm() */ 48 CHKERRQ(MatNorm(A,NORM_FROBENIUS,&s1norm)); 49 CHKERRQ(MatNorm(B,NORM_FROBENIUS,&s2norm)); 50 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 51 if (rnorm>tol) { 52 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 53 } 54 CHKERRQ(MatNorm(A,NORM_INFINITY,&s1norm)); 55 CHKERRQ(MatNorm(B,NORM_INFINITY,&s2norm)); 56 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 57 if (rnorm>tol) { 58 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 59 } 60 CHKERRQ(MatNorm(A,NORM_1,&s1norm)); 61 CHKERRQ(MatNorm(B,NORM_1,&s2norm)); 62 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 63 if (rnorm>tol) { 64 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 65 } 66 67 /* Test MatMult() */ 68 for (i=0; i<IMAX; i++) { 69 CHKERRQ(VecSetRandom(xx,rand)); 70 CHKERRQ(MatMult(A,xx,s1)); 71 CHKERRQ(MatMult(B,xx,s2)); 72 CHKERRQ(VecAXPY(s2,-1.0,s1)); 73 CHKERRQ(VecNorm(s2,NORM_2,&rnorm)); 74 if (rnorm>tol) { 75 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)rnorm,bs)); 76 } 77 } 78 79 /* test MatMultAdd() */ 80 for (i=0; i<IMAX; i++) { 81 CHKERRQ(VecSetRandom(xx,rand)); 82 CHKERRQ(VecSetRandom(yy,rand)); 83 CHKERRQ(MatMultAdd(A,xx,yy,s1)); 84 CHKERRQ(MatMultAdd(B,xx,yy,s2)); 85 CHKERRQ(VecAXPY(s2,-1.0,s1)); 86 CHKERRQ(VecNorm(s2,NORM_2,&rnorm)); 87 if (rnorm>tol) { 88 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)rnorm,bs)); 89 } 90 } 91 92 /* Test MatMultTranspose() */ 93 for (i=0; i<IMAX; i++) { 94 CHKERRQ(VecSetRandom(xx,rand)); 95 CHKERRQ(MatMultTranspose(A,xx,s1)); 96 CHKERRQ(MatMultTranspose(B,xx,s2)); 97 CHKERRQ(VecNorm(s1,NORM_2,&s1norm)); 98 CHKERRQ(VecNorm(s2,NORM_2,&s2norm)); 99 rnorm = s2norm-s1norm; 100 if (rnorm<-tol || rnorm>tol) { 101 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 102 } 103 } 104 /* Test MatMultTransposeAdd() */ 105 for (i=0; i<IMAX; i++) { 106 CHKERRQ(VecSetRandom(xx,rand)); 107 CHKERRQ(VecSetRandom(yy,rand)); 108 CHKERRQ(MatMultTransposeAdd(A,xx,yy,s1)); 109 CHKERRQ(MatMultTransposeAdd(B,xx,yy,s2)); 110 CHKERRQ(VecNorm(s1,NORM_2,&s1norm)); 111 CHKERRQ(VecNorm(s2,NORM_2,&s2norm)); 112 rnorm = s2norm-s1norm; 113 if (rnorm<-tol || rnorm>tol) { 114 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 115 } 116 } 117 118 /* Check MatGetValues() */ 119 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 120 CHKERRQ(MatGetSize(A,&M,&N)); 121 122 for (i=0; i<IMAX; i++) { 123 /* Create random row numbers ad col numbers */ 124 CHKERRQ(PetscRandomGetValue(rand,&v)); 125 cols[0] = (int)(PetscRealPart(v)*N); 126 CHKERRQ(PetscRandomGetValue(rand,&v)); 127 cols[1] = (int)(PetscRealPart(v)*N); 128 CHKERRQ(PetscRandomGetValue(rand,&v)); 129 rows[0] = rstart + (int)(PetscRealPart(v)*m); 130 CHKERRQ(PetscRandomGetValue(rand,&v)); 131 rows[1] = rstart + (int)(PetscRealPart(v)*m); 132 133 CHKERRQ(MatGetValues(A,2,rows,2,cols,vals1)); 134 CHKERRQ(MatGetValues(B,2,rows,2,cols,vals2)); 135 136 for (j=0; j<4; j++) { 137 if (vals1[j] != vals2[j]) { 138 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d]: Error: MatGetValues rstart = %2" PetscInt_FMT " row = %2" PetscInt_FMT " col = %2" PetscInt_FMT " val1 = %e val2 = %e bs = %" PetscInt_FMT "\n",rank,rstart,rows[j/2],cols[j%2],(double)PetscRealPart(vals1[j]),(double)PetscRealPart(vals2[j]),bs)); 139 } 140 } 141 } 142 143 /* Test MatGetRow()/ MatRestoreRow() */ 144 for (ct=0; ct<100; ct++) { 145 CHKERRQ(PetscRandomGetValue(rand,&v)); 146 row = rstart + (PetscInt)(PetscRealPart(v)*m); 147 CHKERRQ(MatGetRow(A,row,&ncols1,&cols1,&v1)); 148 CHKERRQ(MatGetRow(B,row,&ncols2,&cols2,&v2)); 149 150 for (i=0,j=0; i<ncols1 && j<ncols2; j++) { 151 while (cols2[j] != cols1[i]) i++; 152 PetscCheckFalse(v1[i] != v2[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect."); 153 } 154 PetscCheckFalse(j<ncols2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect"); 155 156 CHKERRQ(MatRestoreRow(A,row,&ncols1,&cols1,&v1)); 157 CHKERRQ(MatRestoreRow(B,row,&ncols2,&cols2,&v2)); 158 } 159 160 /* Test MatConvert() */ 161 CHKERRQ(MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C)); 162 163 /* See if MatMult Says both are same */ 164 for (i=0; i<IMAX; i++) { 165 CHKERRQ(VecSetRandom(xx,rand)); 166 CHKERRQ(MatMult(A,xx,s1)); 167 CHKERRQ(MatMult(C,xx,s2)); 168 CHKERRQ(VecNorm(s1,NORM_2,&s1norm)); 169 CHKERRQ(VecNorm(s2,NORM_2,&s2norm)); 170 rnorm = s2norm-s1norm; 171 if (rnorm<-tol || rnorm>tol) { 172 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 173 } 174 } 175 CHKERRQ(MatDestroy(&C)); 176 177 /* Test MatTranspose() */ 178 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&At)); 179 CHKERRQ(MatTranspose(B,MAT_INITIAL_MATRIX,&Bt)); 180 for (i=0; i<IMAX; i++) { 181 CHKERRQ(VecSetRandom(xx,rand)); 182 CHKERRQ(MatMult(At,xx,s1)); 183 CHKERRQ(MatMult(Bt,xx,s2)); 184 CHKERRQ(VecNorm(s1,NORM_2,&s1norm)); 185 CHKERRQ(VecNorm(s2,NORM_2,&s2norm)); 186 rnorm = s2norm-s1norm; 187 if (rnorm<-tol || rnorm>tol) { 188 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",rank,(double)s1norm,(double)s2norm,bs)); 189 } 190 } 191 CHKERRQ(MatDestroy(&At)); 192 CHKERRQ(MatDestroy(&Bt)); 193 194 CHKERRQ(MatDestroy(&A)); 195 CHKERRQ(MatDestroy(&B)); 196 CHKERRQ(VecDestroy(&xx)); 197 CHKERRQ(VecDestroy(&yy)); 198 CHKERRQ(VecDestroy(&s1)); 199 CHKERRQ(VecDestroy(&s2)); 200 CHKERRQ(PetscRandomDestroy(&rand)); 201 CHKERRQ(PetscFinalize()); 202 return 0; 203 } 204 205 /*TEST 206 207 build: 208 requires: !complex 209 210 test: 211 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 212 nsize: 3 213 args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small 214 215 test: 216 suffix: 2 217 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 218 nsize: 3 219 args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small 220 output_file: output/ex53_1.out 221 222 test: 223 suffix: 3 224 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 225 nsize: 3 226 args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small 227 output_file: output/ex53_1.out 228 229 test: 230 TODO: Matrix row/column sizes are not compatible with block size 231 suffix: 4 232 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 233 nsize: 3 234 args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small 235 output_file: output/ex53_1.out 236 237 test: 238 suffix: 5 239 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 240 nsize: 3 241 args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small 242 output_file: output/ex53_1.out 243 244 test: 245 TODO: Matrix row/column sizes are not compatible with block size 246 suffix: 6 247 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 248 nsize: 3 249 args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small 250 output_file: output/ex53_1.out 251 252 test: 253 TODO: Matrix row/column sizes are not compatible with block size 254 suffix: 7 255 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 256 nsize: 3 257 args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small 258 output_file: output/ex53_1.out 259 260 test: 261 suffix: 8 262 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 263 nsize: 4 264 args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small 265 output_file: output/ex53_1.out 266 267 TEST*/ 268