Actual source code: ex53.c

  1: /*$Id: ex53.c,v 1.25 2001/09/11 16:32:50 bsmith Exp $*/

  3: static char help[] = "Tests the vatious routines in MatMPIBAIJ format.\n";


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

 23:   PetscInitialize(&argc,&args,(char *)0,help);
 24:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 26: #if defined(PETSC_USE_COMPLEX)
 27:   SETERRQ(1,"This example does not work with complex numbers");
 28: #else

 30:  /* Check out if MatLoad() works */
 31:   PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);
 32:   if (!flg) SETERRQ(1,"Input file not specified");
 33:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);
 34:   MatLoad(fd,MATMPIBAIJ,&A);
 35:   PetscViewerDestroy(fd);

 37:   MatConvert(A,MATMPIAIJ,&B);
 38: 
 39:   PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);
 40:   MatGetLocalSize(A,&m,&n);
 41:   VecCreate(PETSC_COMM_WORLD,&xx);
 42:   VecSetSizes(xx,m,PETSC_DECIDE);
 43:   VecSetFromOptions(xx);
 44:   VecDuplicate(xx,&s1);
 45:   VecDuplicate(xx,&s2);
 46:   VecDuplicate(xx,&yy);

 48:   MatGetBlockSize(A,&bs);
 49:   /* Test MatMult() */
 50:   for (i=0; i<IMAX; i++) {
 51:     VecSetRandom(rand,xx);
 52:     MatMult(A,xx,s1);
 53:     MatMult(B,xx,s2);
 54:     VecNorm(s1,NORM_2,&s1norm);
 55:     VecNorm(s2,NORM_2,&s2norm);
 56:     rnorm = s2norm-s1norm;
 57:     if (rnorm<-tol || rnorm>tol) {
 58:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
 59: s1norm,s2norm,bs);
 60:     }
 61:   }
 62:   /* test MatMultAdd() */
 63:   for (i=0; i<IMAX; i++) {
 64:     VecSetRandom(rand,xx);
 65:     VecSetRandom(rand,yy);
 66:     MatMultAdd(A,xx,yy,s1);
 67:     MatMultAdd(B,xx,yy,s2);
 68:     VecNorm(s1,NORM_2,&s1norm);
 69:     VecNorm(s2,NORM_2,&s2norm);
 70:     rnorm = s2norm-s1norm;
 71:     if (rnorm<-tol || rnorm>tol) {
 72:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
 73:     }
 74:   }
 75:     /* Test MatMultTranspose() */
 76:   for (i=0; i<IMAX; i++) {
 77:     VecSetRandom(rand,xx);
 78:     MatMultTranspose(A,xx,s1);
 79:     MatMultTranspose(B,xx,s2);
 80:     VecNorm(s1,NORM_2,&s1norm);
 81:     VecNorm(s2,NORM_2,&s2norm);
 82:     rnorm = s2norm-s1norm;
 83:     if (rnorm<-tol || rnorm>tol) {
 84:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
 85:     }
 86:   }
 87:   /* Test MatMultTransposeAdd() */
 88:   for (i=0; i<IMAX; i++) {
 89:     VecSetRandom(rand,xx);
 90:     VecSetRandom(rand,yy);
 91:     MatMultTransposeAdd(A,xx,yy,s1);
 92:     MatMultTransposeAdd(B,xx,yy,s2);
 93:     VecNorm(s1,NORM_2,&s1norm);
 94:     VecNorm(s2,NORM_2,&s2norm);
 95:     rnorm = s2norm-s1norm;
 96:     if (rnorm<-tol || rnorm>tol) {
 97:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
 98:     }
 99:   }

101:   /* Check MatGetValues() */
102:   MatGetOwnershipRange(A,&rstart,&rend);
103:   MatGetSize(A,&M,&N);


106:   for (i=0; i<IMAX; i++) {
107:     /* Create random row numbers ad col numbers */
108:     PetscRandomGetValue(rand,&v);
109:     cols[0] = (int)(PetscRealPart(v)*N);
110:     PetscRandomGetValue(rand,&v);
111:     cols[1] = (int)(PetscRealPart(v)*N);
112:     PetscRandomGetValue(rand,&v);
113:     rows[0] = rstart + (int)(PetscRealPart(v)*m);
114:     PetscRandomGetValue(rand,&v);
115:     rows[1] = rstart + (int)(PetscRealPart(v)*m);
116: 
117:     MatGetValues(A,2,rows,2,cols,vals1);
118:     MatGetValues(B,2,rows,2,cols,vals2);


121:     for (j=0; j<4; j++) {
122:       if(vals1[j] != vals2[j])
123:         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);
124:     }
125:   }

127:   /* Test MatGetRow()/ MatRestoreRow() */
128:   for (ct=0; ct<100; ct++) {
129:     PetscRandomGetValue(rand,&v);
130:     row  = rstart + (int)(PetscRealPart(v)*m);
131:     MatGetRow(A,row,&ncols1,&cols1,&v1);
132:     MatGetRow(B,row,&ncols2,&cols2,&v2);
133: 
134:     for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
135:       while (cols2[j] != cols1[i]) i++;
136:       if (v1[i] != v2[j]) SETERRQ(1,"MatGetRow() failed - vals incorrect.");
137:     }
138:     if (j<ncols2) SETERRQ(1,"MatGetRow() failed - cols incorrect");
139: 
140:     MatRestoreRow(A,row,&ncols1,&cols1,&v1);
141:     MatRestoreRow(B,row,&ncols2,&cols2,&v2);
142:   }
143: 
144:   /* Test MatConvert() */
145:   MatConvert(A,MATSAME,&C);
146: 
147:   /* See if MatMult Says both are same */
148:   for (i=0; i<IMAX; i++) {
149:     VecSetRandom(rand,xx);
150:     MatMult(A,xx,s1);
151:     MatMult(C,xx,s2);
152:     VecNorm(s1,NORM_2,&s1norm);
153:     VecNorm(s2,NORM_2,&s2norm);
154:     rnorm = s2norm-s1norm;
155:     if (rnorm<-tol || rnorm>tol) {
156:       PetscPrintf(PETSC_COMM_SELF,"Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
157: s1norm,s2norm,bs);
158:     }
159:   }
160:   MatDestroy(C);

162:   /* Test MatTranspose() */
163:   MatTranspose(A,&At);
164:   MatTranspose(B,&Bt);
165:   for (i=0; i<IMAX; i++) {
166:     VecSetRandom(rand,xx);
167:     MatMult(At,xx,s1);
168:     MatMult(Bt,xx,s2);
169:     VecNorm(s1,NORM_2,&s1norm);
170:     VecNorm(s2,NORM_2,&s2norm);
171:     rnorm = s2norm-s1norm;
172:     if (rnorm<-tol || rnorm>tol) {
173:       PetscPrintf(PETSC_COMM_SELF,"Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
174:                   s1norm,s2norm,bs);
175:     }
176:   }
177:   MatDestroy(At);
178:   MatDestroy(Bt);

180:   MatDestroy(A);
181:   MatDestroy(B);
182:   VecDestroy(xx);
183:   VecDestroy(yy);
184:   VecDestroy(s1);
185:   VecDestroy(s2);
186:   PetscRandomDestroy(rand);
187:   PetscFinalize();
188: #endif
189:   return 0;
190: }