Actual source code: ex92.c

  2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatSBAIJ format.\n";

 4:  #include petscmat.h

  8: int main(int argc,char **args)
  9: {
 10:   Mat            A,Atrans,sA,*submatA,*submatsA;
 12:   PetscMPIInt    size,rank;
 13:   PetscInt       bs=1,mbs=11,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,mm,M,N,Mbs;
 14:   PetscScalar    *vals,rval,one=1.0;
 15:   IS             *is1,*is2;
 16:   PetscRandom    rand;
 17:   Vec            xx,s1,s2;
 18:   PetscReal      s1norm,s2norm,rnorm,tol = 1.e-10;
 19:   PetscTruth     flg;
 20:   int            stages[2];

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

 26:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
 27:   PetscOptionsGetInt(PETSC_NULL,"-mat_size",&mbs,PETSC_NULL);
 28:   PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
 29:   PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);

 31:   MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE,
 32:                           PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&A);
 33:   PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);

 35:   MatGetOwnershipRange(A,&rstart,&rend);
 36:   MatGetSize(A,&M,&N);
 37:   Mbs  = M/bs;

 39:   PetscMalloc(bs*sizeof(PetscInt),&rows);
 40:   PetscMalloc(bs*sizeof(PetscInt),&cols);
 41:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 42:   PetscMalloc(M*sizeof(PetscScalar),&idx);

 44:   /* Now set blocks of values */
 45:   for (j=0; j<bs*bs; j++) vals[j] = 0.0;
 46:   for (i=0; i<Mbs; i++){
 47:     cols[0] = i*bs; rows[0] = i*bs;
 48:     for (j=1; j<bs; j++) {
 49:       rows[j] = rows[j-1]+1;
 50:       cols[j] = cols[j-1]+1;
 51:     }
 52:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 53:   }
 54:   /* second, add random blocks */
 55:   for (i=0; i<20*bs; i++) {
 56:       PetscRandomGetValue(rand,&rval);
 57:       cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
 58:       PetscRandomGetValue(rand,&rval);
 59:       rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
 60:       for (j=1; j<bs; j++) {
 61:         rows[j] = rows[j-1]+1;
 62:         cols[j] = cols[j-1]+1;
 63:       }

 65:       for (j=0; j<bs*bs; j++) {
 66:         PetscRandomGetValue(rand,&rval);
 67:         vals[j] = rval;
 68:       }
 69:       MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 70:   }

 72:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 73:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 74: 
 75:   /* make A a symmetric matrix: A <- A^T + A */
 76:   MatTranspose(A, &Atrans);
 77:   MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
 78:   MatDestroy(Atrans);
 79:   MatTranspose(A, &Atrans);
 80:   MatEqual(A, Atrans, &flg);
 81:   if (!flg) {
 82:     SETERRQ(1,"A+A^T is non-symmetric");
 83:   }
 84:   MatDestroy(Atrans);
 85:   /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */

 87:   /* create a SeqSBAIJ matrix sA (= A) */
 88:   MatConvert(A,MATMPISBAIJ,MAT_INITIAL_MATRIX,&sA);
 89:   /* MatView(sA,PETSC_VIEWER_STDOUT_WORLD); */

 91:   /* Test sA==A through MatMult() */
 92:   for (i=0; i<nd; i++) {
 93:     MatGetLocalSize(A, &mm,PETSC_NULL);
 94:     VecCreate(PETSC_COMM_WORLD,&xx);
 95:     VecSetSizes(xx,mm,PETSC_DECIDE);
 96:     VecSetFromOptions(xx);
 97:     VecDuplicate(xx,&s1);
 98:     VecDuplicate(xx,&s2);
 99:     for (j=0; j<3; j++) {
100:       VecSetRandom(xx,rand);
101:       MatMult(A,xx,s1);
102:       MatMult(sA,xx,s2);
103:       VecNorm(s1,NORM_2,&s1norm);
104:       VecNorm(s2,NORM_2,&s2norm);
105:       rnorm = s2norm-s1norm;
106:       if (rnorm<-tol || rnorm>tol) {
107:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
108:       }
109:     }
110:     VecDestroy(xx);
111:     VecDestroy(s1);
112:     VecDestroy(s2);
113:   }
114: 
115:   /* Test MatIncreaseOverlap() */
116:   PetscMalloc(nd*sizeof(IS **),&is1);
117:   PetscMalloc(nd*sizeof(IS **),&is2);

119:   for (i=0; i<nd; i++) {
120:     PetscRandomGetValue(rand,&rval);
121:     sz = (PetscInt)((0.5 + 0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
122:     sz /= (size*nd*10);
123:     /*
124:     if (!rank){
125:       PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
126:     } 
127:     */
128:     for (j=0; j<sz; j++) {
129:       PetscRandomGetValue(rand,&rval);
130:       idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
131:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
132:     }
133:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is1+i);
134:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is2+i);
135:   }

137:   PetscLogStageRegister(&stages[0],"MatOv_SBAIJ");
138:   PetscLogStageRegister(&stages[1],"MatOv_ BAIJ");

140:   PetscLogStagePush(stages[0]);
141:   MatIncreaseOverlap(sA,nd,is2,ov);
142:   PetscLogStagePop();

144:   PetscLogStagePush(stages[1]);
145:   MatIncreaseOverlap(A,nd,is1,ov);
146:   PetscLogStagePop();

148:   for (i=0; i<nd; ++i) {
149:     ISEqual(is1[i],is2[i],&flg);
150:     if (!flg ){
151:       if (rank == size){
152:         ISSort(is1[i]);
153:         ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
154:         ISSort(is2[i]);
155:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
156:       }
157:       SETERRQ1(1,"i=%D, is1 != is2",i);
158:     }
159:   }
160: 
161:   for (i=0; i<nd; ++i) {
162:     ISSort(is1[i]);
163:     ISSort(is2[i]);
164:   }
165:   MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
166:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);

168:   /* Test MatMult() */
169:   for (i=0; i<nd; i++) {
170:     MatGetSize(submatA[i],&mm,PETSC_NULL);
171:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
172:     VecDuplicate(xx,&s1);
173:     VecDuplicate(xx,&s2);
174:     for (j=0; j<3; j++) {
175:       VecSetRandom(xx,rand);
176:       MatMult(submatA[i],xx,s1);
177:       MatMult(submatsA[i],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:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
183:       }
184:     }
185:     VecDestroy(xx);
186:     VecDestroy(s1);
187:     VecDestroy(s2);
188:   }

190:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
191:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
192:   MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);

194:   /* Test MatMult() */
195:   for (i=0; i<nd; i++) {
196:     MatGetSize(submatA[i],&mm,PETSC_NULL);
197:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
198:     VecDuplicate(xx,&s1);
199:     VecDuplicate(xx,&s2);
200:     for (j=0; j<3; j++) {
201:       VecSetRandom(xx,rand);
202:       MatMult(submatA[i],xx,s1);
203:       MatMult(submatsA[i],xx,s2);
204:       VecNorm(s1,NORM_2,&s1norm);
205:       VecNorm(s2,NORM_2,&s2norm);
206:       rnorm = s2norm-s1norm;
207:       if (rnorm<-tol || rnorm>tol) {
208:         PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
209:       }
210:     }
211:     VecDestroy(xx);
212:     VecDestroy(s1);
213:     VecDestroy(s2);
214:   }

216:   /* Free allocated memory */
217:   for (i=0; i<nd; ++i) {
218:     ISDestroy(is1[i]);
219:     ISDestroy(is2[i]);
220:     MatDestroy(submatA[i]);
221:     MatDestroy(submatsA[i]);
222:   }
223:   PetscFree(submatA);
224:   PetscFree(submatsA);
225:   PetscFree(is1);
226:   PetscFree(is2);
227:   PetscFree(idx);
228:   PetscFree(rows);
229:   PetscFree(cols);
230:   PetscFree(vals);
231:   MatDestroy(A);
232:   MatDestroy(sA);
233:   PetscRandomDestroy(rand);
234:   PetscFinalize();
235:   return 0;
236: }