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 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 24 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 25 26 /* Check out if MatLoad() works */ 27 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);CHKERRQ(ierr); 28 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Input file not specified"); 29 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr); 30 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 31 ierr = MatSetType(A,MATBAIJ);CHKERRQ(ierr); 32 ierr = MatLoad(A,fd);CHKERRQ(ierr); 33 34 ierr = MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 35 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 36 37 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 38 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 39 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 40 ierr = VecCreate(PETSC_COMM_WORLD,&xx);CHKERRQ(ierr); 41 ierr = VecSetSizes(xx,m,PETSC_DECIDE);CHKERRQ(ierr); 42 ierr = VecSetFromOptions(xx);CHKERRQ(ierr); 43 ierr = VecDuplicate(xx,&s1);CHKERRQ(ierr); 44 ierr = VecDuplicate(xx,&s2);CHKERRQ(ierr); 45 ierr = VecDuplicate(xx,&yy);CHKERRQ(ierr); 46 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 47 48 /* Test MatNorm() */ 49 ierr = MatNorm(A,NORM_FROBENIUS,&s1norm);CHKERRQ(ierr); 50 ierr = MatNorm(B,NORM_FROBENIUS,&s2norm);CHKERRQ(ierr); 51 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 52 if (rnorm>tol) { 53 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);CHKERRQ(ierr); 54 } 55 ierr = MatNorm(A,NORM_INFINITY,&s1norm);CHKERRQ(ierr); 56 ierr = MatNorm(B,NORM_INFINITY,&s2norm);CHKERRQ(ierr); 57 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 58 if (rnorm>tol) { 59 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);CHKERRQ(ierr); 60 } 61 ierr = MatNorm(A,NORM_1,&s1norm);CHKERRQ(ierr); 62 ierr = MatNorm(B,NORM_1,&s2norm);CHKERRQ(ierr); 63 rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm; 64 if (rnorm>tol) { 65 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); 66 } 67 68 /* Test MatMult() */ 69 for (i=0; i<IMAX; i++) { 70 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 71 ierr = MatMult(A,xx,s1);CHKERRQ(ierr); 72 ierr = MatMult(B,xx,s2);CHKERRQ(ierr); 73 ierr = VecAXPY(s2,-1.0,s1);CHKERRQ(ierr); 74 ierr = VecNorm(s2,NORM_2,&rnorm);CHKERRQ(ierr); 75 if (rnorm>tol) { 76 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);CHKERRQ(ierr); 77 } 78 } 79 80 /* test MatMultAdd() */ 81 for (i=0; i<IMAX; i++) { 82 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 83 ierr = VecSetRandom(yy,rand);CHKERRQ(ierr); 84 ierr = MatMultAdd(A,xx,yy,s1);CHKERRQ(ierr); 85 ierr = MatMultAdd(B,xx,yy,s2);CHKERRQ(ierr); 86 ierr = VecAXPY(s2,-1.0,s1);CHKERRQ(ierr); 87 ierr = VecNorm(s2,NORM_2,&rnorm);CHKERRQ(ierr); 88 if (rnorm>tol) { 89 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);CHKERRQ(ierr); 90 } 91 } 92 93 /* Test MatMultTranspose() */ 94 for (i=0; i<IMAX; i++) { 95 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 96 ierr = MatMultTranspose(A,xx,s1);CHKERRQ(ierr); 97 ierr = MatMultTranspose(B,xx,s2);CHKERRQ(ierr); 98 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 99 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 100 rnorm = s2norm-s1norm; 101 if (rnorm<-tol || rnorm>tol) { 102 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);CHKERRQ(ierr); 103 } 104 } 105 /* Test MatMultTransposeAdd() */ 106 for (i=0; i<IMAX; i++) { 107 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 108 ierr = VecSetRandom(yy,rand);CHKERRQ(ierr); 109 ierr = MatMultTransposeAdd(A,xx,yy,s1);CHKERRQ(ierr); 110 ierr = MatMultTransposeAdd(B,xx,yy,s2);CHKERRQ(ierr); 111 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 112 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 113 rnorm = s2norm-s1norm; 114 if (rnorm<-tol || rnorm>tol) { 115 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);CHKERRQ(ierr); 116 } 117 } 118 119 /* Check MatGetValues() */ 120 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 121 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 122 123 for (i=0; i<IMAX; i++) { 124 /* Create random row numbers ad col numbers */ 125 ierr = PetscRandomGetValue(rand,&v);CHKERRQ(ierr); 126 cols[0] = (int)(PetscRealPart(v)*N); 127 ierr = PetscRandomGetValue(rand,&v);CHKERRQ(ierr); 128 cols[1] = (int)(PetscRealPart(v)*N); 129 ierr = PetscRandomGetValue(rand,&v);CHKERRQ(ierr); 130 rows[0] = rstart + (int)(PetscRealPart(v)*m); 131 ierr = PetscRandomGetValue(rand,&v);CHKERRQ(ierr); 132 rows[1] = rstart + (int)(PetscRealPart(v)*m); 133 134 ierr = MatGetValues(A,2,rows,2,cols,vals1);CHKERRQ(ierr); 135 ierr = MatGetValues(B,2,rows,2,cols,vals2);CHKERRQ(ierr); 136 137 for (j=0; j<4; j++) { 138 if (vals1[j] != vals2[j]) { 139 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); 140 } 141 } 142 } 143 144 /* Test MatGetRow()/ MatRestoreRow() */ 145 for (ct=0; ct<100; ct++) { 146 ierr = PetscRandomGetValue(rand,&v);CHKERRQ(ierr); 147 row = rstart + (int)(PetscRealPart(v)*m); 148 ierr = MatGetRow(A,row,&ncols1,&cols1,&v1);CHKERRQ(ierr); 149 ierr = MatGetRow(B,row,&ncols2,&cols2,&v2);CHKERRQ(ierr); 150 151 for (i=0,j=0; i<ncols1 && j<ncols2; j++) { 152 while (cols2[j] != cols1[i]) i++; 153 if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect."); 154 } 155 if (j<ncols2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect"); 156 157 ierr = MatRestoreRow(A,row,&ncols1,&cols1,&v1);CHKERRQ(ierr); 158 ierr = MatRestoreRow(B,row,&ncols2,&cols2,&v2);CHKERRQ(ierr); 159 } 160 161 /* Test MatConvert() */ 162 ierr = MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 163 164 /* See if MatMult Says both are same */ 165 for (i=0; i<IMAX; i++) { 166 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 167 ierr = MatMult(A,xx,s1);CHKERRQ(ierr); 168 ierr = MatMult(C,xx,s2);CHKERRQ(ierr); 169 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 170 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 171 rnorm = s2norm-s1norm; 172 if (rnorm<-tol || rnorm>tol) { 173 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); 174 } 175 } 176 ierr = MatDestroy(&C);CHKERRQ(ierr); 177 178 /* Test MatTranspose() */ 179 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr); 180 ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);CHKERRQ(ierr); 181 for (i=0; i<IMAX; i++) { 182 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 183 ierr = MatMult(At,xx,s1);CHKERRQ(ierr); 184 ierr = MatMult(Bt,xx,s2);CHKERRQ(ierr); 185 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 186 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 187 rnorm = s2norm-s1norm; 188 if (rnorm<-tol || rnorm>tol) { 189 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); 190 } 191 } 192 ierr = MatDestroy(&At);CHKERRQ(ierr); 193 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 194 195 ierr = MatDestroy(&A);CHKERRQ(ierr); 196 ierr = MatDestroy(&B);CHKERRQ(ierr); 197 ierr = VecDestroy(&xx);CHKERRQ(ierr); 198 ierr = VecDestroy(&yy);CHKERRQ(ierr); 199 ierr = VecDestroy(&s1);CHKERRQ(ierr); 200 ierr = VecDestroy(&s2);CHKERRQ(ierr); 201 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 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