Actual source code: ex91.c

  1: /*$Id: ex91.c,v 1.22 2001/08/07 21:30:08 bsmith Exp $*/

  3: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";

 5:  #include petscmat.h

  9: int main(int argc,char **args)
 10: {
 11:   Mat         A,Atrans,sA,*submatA,*submatsA;
 12:   int         bs=1,m=43,ov=1,i,j,k,*rows,*cols,ierr,M,nd=5,*idx,size,mm,nn;
 13:   PetscScalar *vals,rval,one=1.0;
 14:   IS          *is1,*is2;
 15:   PetscRandom rand;
 16:   Vec         xx,s1,s2;
 17:   PetscReal   s1norm,s2norm,rnorm,tol = 1.e-10;
 18:   PetscTruth  flg;

 20:   PetscInitialize(&argc,&args,(char *)0,help);
 21: 

 23:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
 24:   PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
 25:   PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
 26:   PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);

 28:   /* create a SeqBAIJ matrix A */
 29:   M    = m*bs;
 30:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
 31:   PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);

 33:   PetscMalloc(bs*sizeof(int),&rows);
 34:   PetscMalloc(bs*sizeof(int),&cols);
 35:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 36:   PetscMalloc(M*sizeof(PetscScalar),&idx);
 37: 
 38:   /* Now set blocks of random values */
 39:   /* first, set diagonal blocks as zero */
 40:   for (j=0; j<bs*bs; j++) vals[j] = 0.0;
 41:   for (i=0; i<m; i++){
 42:     cols[0] = i*bs; rows[0] = i*bs;
 43:     for (j=1; j<bs; j++) {
 44:       rows[j] = rows[j-1]+1;
 45:       cols[j] = cols[j-1]+1;
 46:     }
 47:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 48:   }
 49:   /* second, add random blocks */
 50:   for (i=0; i<20*bs; i++) {
 51:       PetscRandomGetValue(rand,&rval);
 52:       cols[0] = bs*(int)(PetscRealPart(rval)*m);
 53:       PetscRandomGetValue(rand,&rval);
 54:       rows[0] = bs*(int)(PetscRealPart(rval)*m);
 55:       for (j=1; j<bs; j++) {
 56:         rows[j] = rows[j-1]+1;
 57:         cols[j] = cols[j-1]+1;
 58:       }

 60:       for (j=0; j<bs*bs; j++) {
 61:         PetscRandomGetValue(rand,&rval);
 62:         vals[j] = rval;
 63:       }
 64:       MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 65:   }

 67:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 70:   /* make A a symmetric matrix: A <- A^T + A */
 71:   MatTranspose(A, &Atrans);
 72:   MatAXPY(&one,Atrans,A,DIFFERENT_NONZERO_PATTERN);
 73:   MatDestroy(Atrans);
 74:   MatTranspose(A, &Atrans);
 75:   MatEqual(A, Atrans, &flg);
 76:   if (!flg) {
 77:     SETERRQ(1,"A+A^T is non-symmetric");
 78:   }
 79:   MatDestroy(Atrans);

 81:   /* create a SeqSBAIJ matrix sA (= A) */
 82:   MatConvert(A,MATSEQSBAIJ,&sA);
 83: 
 84:   /* Test sA==A through MatMult() */
 85:   for (i=0; i<nd; i++) {
 86:     MatGetSize(A,&mm,&nn);
 87:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
 88:     VecDuplicate(xx,&s1);
 89:     VecDuplicate(xx,&s2);
 90:     for (j=0; j<3; j++) {
 91:       VecSetRandom(rand,xx);
 92:       MatMult(A,xx,s1);
 93:       MatMult(sA,xx,s2);
 94:       VecNorm(s1,NORM_2,&s1norm);
 95:       VecNorm(s2,NORM_2,&s2norm);
 96:       rnorm = s2norm-s1norm;
 97:       if (rnorm<-tol || rnorm>tol) {
 98:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
 99:       }
100:     }
101:     VecDestroy(xx);
102:     VecDestroy(s1);
103:     VecDestroy(s2);
104:   }

106:   /* Test MatIncreaseOverlap() */
107:   PetscMalloc(nd*sizeof(IS **),&is1);
108:   PetscMalloc(nd*sizeof(IS **),&is2);

110: 
111:   for (i=0; i<nd; i++) {
112:     PetscRandomGetValue(rand,&rval);
113:     size = (int)(PetscRealPart(rval)*m);
114:     for (j=0; j<size; j++) {
115:       PetscRandomGetValue(rand,&rval);
116:       idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
117:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
118:     }
119:     ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,is1+i);
120:     ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,is2+i);
121:   }
122:   /* for debugging */
123:   /*
124:   MatView(A,PETSC_VIEWER_STDOUT_SELF); 
125:   MatView(sA,PETSC_VIEWER_STDOUT_SELF); 
126:   */

128:   MatIncreaseOverlap(A,nd,is1,ov);
129:   MatIncreaseOverlap(sA,nd,is2,ov);

131:   for (i=0; i<nd; ++i) {
132:     ISSort(is1[i]);
133:     ISSort(is2[i]);
134:   }

136:   for (i=0; i<nd; ++i) {
137:     ISEqual(is1[i],is2[i],&flg);
138:     if (!flg){
139:       /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
140:          ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
141:       SETERRQ1(1,"i=%d, is1 != is2",i);
142:     }
143:   }
144: 
145:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
146:   MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);

148:   /* Test MatMult() */
149:   for (i=0; i<nd; i++) {
150:     MatGetSize(submatA[i],&mm,&nn);
151:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
152:     VecDuplicate(xx,&s1);
153:     VecDuplicate(xx,&s2);
154:     for (j=0; j<3; j++) {
155:       VecSetRandom(rand,xx);
156:       MatMult(submatA[i],xx,s1);
157:       MatMult(submatsA[i],xx,s2);
158:       VecNorm(s1,NORM_2,&s1norm);
159:       VecNorm(s2,NORM_2,&s2norm);
160:       rnorm = s2norm-s1norm;
161:       if (rnorm<-tol || rnorm>tol) {
162:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
163:       }
164:     }
165:     VecDestroy(xx);
166:     VecDestroy(s1);
167:     VecDestroy(s2);
168:   }

170:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
171:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
172:   MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
173: 
174:   /* Test MatMult() */
175:   for (i=0; i<nd; i++) {
176:     MatGetSize(submatA[i],&mm,&nn);
177:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
178:     VecDuplicate(xx,&s1);
179:     VecDuplicate(xx,&s2);
180:     for (j=0; j<3; j++) {
181:       VecSetRandom(rand,xx);
182:       MatMult(submatA[i],xx,s1);
183:       MatMult(submatsA[i],xx,s2);
184:       VecNorm(s1,NORM_2,&s1norm);
185:       VecNorm(s2,NORM_2,&s2norm);
186:       rnorm = s2norm-s1norm;
187:       if (rnorm<-tol || rnorm>tol) {
188:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
189:       }
190:     }
191:     VecDestroy(xx);
192:     VecDestroy(s1);
193:     VecDestroy(s2);
194:   }
195: 
196:   /* Free allocated memory */
197:   for (i=0; i<nd; ++i) {
198:     ISDestroy(is1[i]);
199:     ISDestroy(is2[i]);
200: 
201:     MatDestroy(submatA[i]);
202:     MatDestroy(submatsA[i]);
203: 
204:   }
205: 
206:   PetscFree(submatA);
207:   PetscFree(submatsA);
208: 
209:   PetscFree(is1);
210:   PetscFree(is2);
211:   PetscFree(idx);
212:   PetscFree(rows);
213:   PetscFree(cols);
214:   PetscFree(vals);
215:   MatDestroy(A);
216:   MatDestroy(sA);
217:   PetscRandomDestroy(rand);
218:   PetscFinalize();
219:   return 0;
220: }