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