Actual source code: ex53.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Tests the vatious routines in MatMPIBAIJ format.\n";


  5: #include <petscmat.h>
  6: #define IMAX 15
  9: int main(int argc,char **args)
 10: {
 11:   Mat               A,B,C,At,Bt;
 12:   PetscViewer       fd;
 13:   char              file[PETSC_MAX_PATH_LEN];
 14:   PetscRandom       rand;
 15:   Vec               xx,yy,s1,s2;
 16:   PetscReal         s1norm,s2norm,rnorm,tol=1.e-10;
 17:   PetscInt          rstart,rend,rows[2],cols[2],m,n,i,j,M,N,ct,row,ncols1,ncols2,bs;
 18:   PetscMPIInt       rank,size;
 19:   PetscErrorCode    ierr;
 20:   const PetscInt    *cols1,*cols2;
 21:   PetscScalar       vals1[4],vals2[4],v;
 22:   const PetscScalar *v1,*v2;
 23:   PetscBool         flg;

 25:   PetscInitialize(&argc,&args,(char*)0,help);
 26:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 29: #if defined(PETSC_USE_COMPLEX)
 30:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
 31: #else

 33:   /* Check out if MatLoad() works */
 34:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 35:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Input file not specified");
 36:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 37:   MatCreate(PETSC_COMM_WORLD,&A);
 38:   MatSetType(A,MATBAIJ);
 39:   MatLoad(A,fd);

 41:   MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
 42:   PetscViewerDestroy(&fd);

 44:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 45:   PetscRandomSetFromOptions(rand);
 46:   MatGetLocalSize(A,&m,&n);
 47:   VecCreate(PETSC_COMM_WORLD,&xx);
 48:   VecSetSizes(xx,m,PETSC_DECIDE);
 49:   VecSetFromOptions(xx);
 50:   VecDuplicate(xx,&s1);
 51:   VecDuplicate(xx,&s2);
 52:   VecDuplicate(xx,&yy);
 53:   MatGetBlockSize(A,&bs);

 55:   /* Test MatNorm() */
 56:   MatNorm(A,NORM_FROBENIUS,&s1norm);
 57:   MatNorm(B,NORM_FROBENIUS,&s2norm);
 58:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 59:   if (rnorm>tol) {
 60:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 61:   }
 62:   MatNorm(A,NORM_INFINITY,&s1norm);
 63:   MatNorm(B,NORM_INFINITY,&s2norm);
 64:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 65:   if (rnorm>tol) {
 66:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 67:   }
 68:   MatNorm(A,NORM_1,&s1norm);
 69:   MatNorm(B,NORM_1,&s2norm);
 70:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 71:   if (rnorm>tol) {
 72:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 73:   }

 75:   /* Test MatMult() */
 76:   for (i=0; i<IMAX; i++) {
 77:     VecSetRandom(xx,rand);
 78:     MatMult(A,xx,s1);
 79:     MatMult(B,xx,s2);
 80:     VecAXPY(s2,-1.0,s1);
 81:     VecNorm(s2,NORM_2,&rnorm);
 82:     if (rnorm>tol) {
 83:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
 84:     }
 85:   }

 87:   /* test MatMultAdd() */
 88:   for (i=0; i<IMAX; i++) {
 89:     VecSetRandom(xx,rand);
 90:     VecSetRandom(yy,rand);
 91:     MatMultAdd(A,xx,yy,s1);
 92:     MatMultAdd(B,xx,yy,s2);
 93:     VecAXPY(s2,-1.0,s1);
 94:     VecNorm(s2,NORM_2,&rnorm);
 95:     if (rnorm>tol) {
 96:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
 97:     }
 98:   }

100:   /* Test MatMultTranspose() */
101:   for (i=0; i<IMAX; i++) {
102:     VecSetRandom(xx,rand);
103:     MatMultTranspose(A,xx,s1);
104:     MatMultTranspose(B,xx,s2);
105:     VecNorm(s1,NORM_2,&s1norm);
106:     VecNorm(s2,NORM_2,&s2norm);
107:     rnorm = s2norm-s1norm;
108:     if (rnorm<-tol || rnorm>tol) {
109:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
110:     }
111:   }
112:   /* Test MatMultTransposeAdd() */
113:   for (i=0; i<IMAX; i++) {
114:     VecSetRandom(xx,rand);
115:     VecSetRandom(yy,rand);
116:     MatMultTransposeAdd(A,xx,yy,s1);
117:     MatMultTransposeAdd(B,xx,yy,s2);
118:     VecNorm(s1,NORM_2,&s1norm);
119:     VecNorm(s2,NORM_2,&s2norm);
120:     rnorm = s2norm-s1norm;
121:     if (rnorm<-tol || rnorm>tol) {
122:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
123:     }
124:   }

126:   /* Check MatGetValues() */
127:   MatGetOwnershipRange(A,&rstart,&rend);
128:   MatGetSize(A,&M,&N);


131:   for (i=0; i<IMAX; i++) {
132:     /* Create random row numbers ad col numbers */
133:     PetscRandomGetValue(rand,&v);
134:     cols[0] = (int)(PetscRealPart(v)*N);
135:     PetscRandomGetValue(rand,&v);
136:     cols[1] = (int)(PetscRealPart(v)*N);
137:     PetscRandomGetValue(rand,&v);
138:     rows[0] = rstart + (int)(PetscRealPart(v)*m);
139:     PetscRandomGetValue(rand,&v);
140:     rows[1] = rstart + (int)(PetscRealPart(v)*m);

142:     MatGetValues(A,2,rows,2,cols,vals1);
143:     MatGetValues(B,2,rows,2,cols,vals2);


146:     for (j=0; j<4; j++) {
147:       if (vals1[j] != vals2[j]) {
148:         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);
149:       }
150:     }
151:   }

153:   /* Test MatGetRow()/ MatRestoreRow() */
154:   for (ct=0; ct<100; ct++) {
155:     PetscRandomGetValue(rand,&v);
156:     row  = rstart + (int)(PetscRealPart(v)*m);
157:     MatGetRow(A,row,&ncols1,&cols1,&v1);
158:     MatGetRow(B,row,&ncols2,&cols2,&v2);

160:     for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
161:       while (cols2[j] != cols1[i]) i++;
162:       if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - vals incorrect.");
163:     }
164:     if (j<ncols2) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - cols incorrect");

166:     MatRestoreRow(A,row,&ncols1,&cols1,&v1);
167:     MatRestoreRow(B,row,&ncols2,&cols2,&v2);
168:   }

170:   /* Test MatConvert() */
171:   MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);

173:   /* See if MatMult Says both are same */
174:   for (i=0; i<IMAX; i++) {
175:     VecSetRandom(xx,rand);
176:     MatMult(A,xx,s1);
177:     MatMult(C,xx,s2);
178:     VecNorm(s1,NORM_2,&s1norm);
179:     VecNorm(s2,NORM_2,&s2norm);
180:     rnorm = s2norm-s1norm;
181:     if (rnorm<-tol || rnorm>tol) {
182:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
183:     }
184:   }
185:   MatDestroy(&C);

187:   /* Test MatTranspose() */
188:   MatTranspose(A,MAT_INITIAL_MATRIX,&At);
189:   MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);
190:   for (i=0; i<IMAX; i++) {
191:     VecSetRandom(xx,rand);
192:     MatMult(At,xx,s1);
193:     MatMult(Bt,xx,s2);
194:     VecNorm(s1,NORM_2,&s1norm);
195:     VecNorm(s2,NORM_2,&s2norm);
196:     rnorm = s2norm-s1norm;
197:     if (rnorm<-tol || rnorm>tol) {
198:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
199:     }
200:   }
201:   MatDestroy(&At);
202:   MatDestroy(&Bt);

204:   MatDestroy(&A);
205:   MatDestroy(&B);
206:   VecDestroy(&xx);
207:   VecDestroy(&yy);
208:   VecDestroy(&s1);
209:   VecDestroy(&s2);
210:   PetscRandomDestroy(&rand);
211:   PetscFinalize();
212: #endif
213:   return 0;
214: }