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 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 25 26 /* Check out if MatLoad() works */ 27 PetscCall(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg)); 28 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Input file not specified"); 29 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 30 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 31 PetscCall(MatSetType(A,MATBAIJ)); 32 PetscCall(MatLoad(A,fd)); 33 34 PetscCall(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B)); 35 PetscCall(PetscViewerDestroy(&fd)); 36 37 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 38 PetscCall(PetscRandomSetFromOptions(rand)); 39 PetscCall(MatGetLocalSize(A,&m,&n)); 40 PetscCall(VecCreate(PETSC_COMM_WORLD,&xx)); 41 PetscCall(VecSetSizes(xx,m,PETSC_DECIDE)); 42 PetscCall(VecSetFromOptions(xx)); 43 PetscCall(VecDuplicate(xx,&s1)); 44 PetscCall(VecDuplicate(xx,&s2)); 45 PetscCall(VecDuplicate(xx,&yy)); 46 PetscCall(MatGetBlockSize(A,&bs)); 47 48 /* Test MatNorm() */ 49 PetscCall(MatNorm(A,NORM_FROBENIUS,&s1norm)); 50 PetscCall(MatNorm(B,NORM_FROBENIUS,&s2norm)); 51 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 52 if (rnorm>tol) { 53 PetscCall(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 PetscCall(MatNorm(A,NORM_INFINITY,&s1norm)); 56 PetscCall(MatNorm(B,NORM_INFINITY,&s2norm)); 57 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 58 if (rnorm>tol) { 59 PetscCall(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 PetscCall(MatNorm(A,NORM_1,&s1norm)); 62 PetscCall(MatNorm(B,NORM_1,&s2norm)); 63 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 64 if (rnorm>tol) { 65 PetscCall(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 PetscCall(VecSetRandom(xx,rand)); 71 PetscCall(MatMult(A,xx,s1)); 72 PetscCall(MatMult(B,xx,s2)); 73 PetscCall(VecAXPY(s2,-1.0,s1)); 74 PetscCall(VecNorm(s2,NORM_2,&rnorm)); 75 if (rnorm>tol) { 76 PetscCall(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 PetscCall(VecSetRandom(xx,rand)); 83 PetscCall(VecSetRandom(yy,rand)); 84 PetscCall(MatMultAdd(A,xx,yy,s1)); 85 PetscCall(MatMultAdd(B,xx,yy,s2)); 86 PetscCall(VecAXPY(s2,-1.0,s1)); 87 PetscCall(VecNorm(s2,NORM_2,&rnorm)); 88 if (rnorm>tol) { 89 PetscCall(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 PetscCall(VecSetRandom(xx,rand)); 96 PetscCall(MatMultTranspose(A,xx,s1)); 97 PetscCall(MatMultTranspose(B,xx,s2)); 98 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 99 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 100 rnorm = s2norm-s1norm; 101 if (rnorm<-tol || rnorm>tol) { 102 PetscCall(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 PetscCall(VecSetRandom(xx,rand)); 108 PetscCall(VecSetRandom(yy,rand)); 109 PetscCall(MatMultTransposeAdd(A,xx,yy,s1)); 110 PetscCall(MatMultTransposeAdd(B,xx,yy,s2)); 111 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 112 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 113 rnorm = s2norm-s1norm; 114 if (rnorm<-tol || rnorm>tol) { 115 PetscCall(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 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 121 PetscCall(MatGetSize(A,&M,&N)); 122 123 for (i=0; i<IMAX; i++) { 124 /* Create random row numbers ad col numbers */ 125 PetscCall(PetscRandomGetValue(rand,&v)); 126 cols[0] = (int)(PetscRealPart(v)*N); 127 PetscCall(PetscRandomGetValue(rand,&v)); 128 cols[1] = (int)(PetscRealPart(v)*N); 129 PetscCall(PetscRandomGetValue(rand,&v)); 130 rows[0] = rstart + (int)(PetscRealPart(v)*m); 131 PetscCall(PetscRandomGetValue(rand,&v)); 132 rows[1] = rstart + (int)(PetscRealPart(v)*m); 133 134 PetscCall(MatGetValues(A,2,rows,2,cols,vals1)); 135 PetscCall(MatGetValues(B,2,rows,2,cols,vals2)); 136 137 for (j=0; j<4; j++) { 138 if (vals1[j] != vals2[j]) { 139 PetscCall(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 PetscCall(PetscRandomGetValue(rand,&v)); 147 row = rstart + (PetscInt)(PetscRealPart(v)*m); 148 PetscCall(MatGetRow(A,row,&ncols1,&cols1,&v1)); 149 PetscCall(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 PetscCheck(v1[i] == v2[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect."); 154 } 155 PetscCheck(j>=ncols2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect"); 156 157 PetscCall(MatRestoreRow(A,row,&ncols1,&cols1,&v1)); 158 PetscCall(MatRestoreRow(B,row,&ncols2,&cols2,&v2)); 159 } 160 161 /* Test MatConvert() */ 162 PetscCall(MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C)); 163 164 /* See if MatMult Says both are same */ 165 for (i=0; i<IMAX; i++) { 166 PetscCall(VecSetRandom(xx,rand)); 167 PetscCall(MatMult(A,xx,s1)); 168 PetscCall(MatMult(C,xx,s2)); 169 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 170 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 171 rnorm = s2norm-s1norm; 172 if (rnorm<-tol || rnorm>tol) { 173 PetscCall(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 PetscCall(MatDestroy(&C)); 177 178 /* Test MatTranspose() */ 179 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&At)); 180 PetscCall(MatTranspose(B,MAT_INITIAL_MATRIX,&Bt)); 181 for (i=0; i<IMAX; i++) { 182 PetscCall(VecSetRandom(xx,rand)); 183 PetscCall(MatMult(At,xx,s1)); 184 PetscCall(MatMult(Bt,xx,s2)); 185 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 186 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 187 rnorm = s2norm-s1norm; 188 if (rnorm<-tol || rnorm>tol) { 189 PetscCall(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 PetscCall(MatDestroy(&At)); 193 PetscCall(MatDestroy(&Bt)); 194 195 PetscCall(MatDestroy(&A)); 196 PetscCall(MatDestroy(&B)); 197 PetscCall(VecDestroy(&xx)); 198 PetscCall(VecDestroy(&yy)); 199 PetscCall(VecDestroy(&s1)); 200 PetscCall(VecDestroy(&s2)); 201 PetscCall(PetscRandomDestroy(&rand)); 202 PetscCall(PetscFinalize()); 203 return 0; 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