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