Actual source code: ex92.c

petsc-3.3-p6 2013-02-11
  2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatSBAIJ format.\n";
  3: /* Example of usage:
  4:       mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat 
  5: */
  6: #include <petscmat.h>

 10: int main(int argc,char **args)
 11: {
 12:   Mat            A,Atrans,sA,*submatA,*submatsA;
 14:   PetscMPIInt    size,rank;
 15:   PetscInt       bs=1,mbs=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs;
 16:   PetscScalar    *vals,rval,one=1.0;
 17:   IS             *is1,*is2;
 18:   PetscRandom    rand;
 19:   PetscBool      flg,TestOverlap,TestSubMat,TestAllcols;
 20:   PetscLogStage  stages[2];
 21:   PetscInt       vid = -1;

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

 27:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
 28:   PetscOptionsGetInt(PETSC_NULL,"-mat_mbs",&mbs,PETSC_NULL);
 29:   PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
 30:   PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
 31:   PetscOptionsGetInt(PETSC_NULL,"-view_id",&vid,PETSC_NULL);
 32:   PetscOptionsHasName(PETSC_NULL, "-test_overlap", &TestOverlap);
 33:   PetscOptionsHasName(PETSC_NULL, "-test_submat", &TestSubMat);
 34:   PetscOptionsHasName(PETSC_NULL, "-test_allcols", &TestAllcols);
 35: 
 36:   MatCreate(PETSC_COMM_WORLD,&A);
 37:   MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
 38:   MatSetType(A,MATBAIJ);
 39:   MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL);
 40:   MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL);

 42:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 43:   PetscRandomSetFromOptions(rand);

 45:   MatGetOwnershipRange(A,&rstart,&rend);
 46:   MatGetSize(A,&M,&N);
 47:   Mbs  = M/bs;

 49:   PetscMalloc(bs*sizeof(PetscInt),&rows);
 50:   PetscMalloc(bs*sizeof(PetscInt),&cols);
 51:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 52:   PetscMalloc(M*sizeof(PetscScalar),&idx);

 54:   /* Now set blocks of values */
 55:   for (j=0; j<bs*bs; j++) vals[j] = 0.0;
 56:   for (i=0; i<Mbs; i++){
 57:     cols[0] = i*bs; rows[0] = i*bs;
 58:     for (j=1; j<bs; j++) {
 59:       rows[j] = rows[j-1]+1;
 60:       cols[j] = cols[j-1]+1;
 61:     }
 62:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 63:   }
 64:   /* second, add random blocks */
 65:   for (i=0; i<20*bs; i++) {
 66:       PetscRandomGetValue(rand,&rval);
 67:       cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
 68:       PetscRandomGetValue(rand,&rval);
 69:       rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
 70:       for (j=1; j<bs; j++) {
 71:         rows[j] = rows[j-1]+1;
 72:         cols[j] = cols[j-1]+1;
 73:       }

 75:       for (j=0; j<bs*bs; j++) {
 76:         PetscRandomGetValue(rand,&rval);
 77:         vals[j] = rval;
 78:       }
 79:       MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 80:   }

 82:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 84: 
 85:   /* make A a symmetric matrix: A <- A^T + A */
 86:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
 87:   MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
 88:   MatDestroy(&Atrans);
 89:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
 90:   MatEqual(A, Atrans, &flg);
 91:   if (flg) {
 92:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 93:   } else {
 94:     SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
 95:   }
 96:   MatDestroy(&Atrans);

 98:   /* create a SeqSBAIJ matrix sA (= A) */
 99:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
100:   if (vid >= 0 && vid < size){
101:     if (!rank) printf("A: \n");
102:     MatView(A,PETSC_VIEWER_STDOUT_WORLD);
103:     if (!rank) printf("sA: \n");
104:     MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
105:   }

107:   /* Test sA==A through MatMult() */
108:   MatMultEqual(A,sA,10,&flg);
109:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG ,"Error in MatConvert(): A != sA");
110: 
111:   /* Test MatIncreaseOverlap() */
112:   PetscMalloc(nd*sizeof(IS **),&is1);
113:   PetscMalloc(nd*sizeof(IS **),&is2);

115:   for (i=0; i<nd; i++) {
116:     if (!TestAllcols){
117:       PetscRandomGetValue(rand,&rval);
118:       sz = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
119: 
120:       for (j=0; j<sz; j++) {
121:         PetscRandomGetValue(rand,&rval);
122:         idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
123:         for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
124:       }
125:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
126:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
127:       if (rank == vid){
128:         PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
129:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
130:       }
131:     } else { /* Test all rows and colums */
132:       sz = M;
133:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);
134:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);

136:       if (rank == vid){
137:         PetscBool      colflag;
138:         ISIdentity(is2[i],&colflag);
139:         printf("[%d] is2[%d], colflag %d\n",rank,i,colflag);
140:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
141:       }
142:     }
143:   }

145:   PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
146:   PetscLogStageRegister("MatOv_BAIJ",&stages[1]);

148:   /* Test MatIncreaseOverlap */
149:   if (TestOverlap){
150:     PetscLogStagePush(stages[0]);
151:     MatIncreaseOverlap(sA,nd,is2,ov);
152:     PetscLogStagePop();

154:     PetscLogStagePush(stages[1]);
155:     MatIncreaseOverlap(A,nd,is1,ov);
156:     PetscLogStagePop();

158:     if (rank == vid){
159:       printf("\n[%d] IS from BAIJ:\n",rank);
160:       ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
161:       printf("\n[%d] IS from SBAIJ:\n",rank);
162:       ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
163:     }

165:     for (i=0; i<nd; ++i) {
166:       ISEqual(is1[i],is2[i],&flg);
167:       if (!flg ){
168:         if (rank == 0){
169:           ISSort(is1[i]);
170:           /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF); */
171:           ISSort(is2[i]);
172:           /* ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
173:         }
174:         SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
175:       }
176:     }
177:   }

179:   /* Test MatGetSubmatrices */
180:   if (TestSubMat){
181:     for(i = 0; i < nd; ++i) {
182:       ISSort(is1[i]);
183:       ISSort(is2[i]);
184:     }
185:     MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
186:     MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);

188:     MatMultEqual(A,sA,10,&flg);
189:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");

191:     /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
192:     MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
193:     MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
194:     MatMultEqual(A,sA,10,&flg);
195:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatGetSubmatrices(): A != sA");
196: 
197:     for (i=0; i<nd; ++i) {
198:       MatDestroy(&submatA[i]);
199:       MatDestroy(&submatsA[i]);
200:     }
201:     PetscFree(submatA);
202:     PetscFree(submatsA);
203:   }

205:   /* Free allocated memory */
206:   for (i=0; i<nd; ++i) {
207:     ISDestroy(&is1[i]);
208:     ISDestroy(&is2[i]);
209:   }
210:   PetscFree(is1);
211:   PetscFree(is2);
212:   PetscFree(idx);
213:   PetscFree(rows);
214:   PetscFree(cols);
215:   PetscFree(vals);
216:   MatDestroy(&A);
217:   MatDestroy(&sA);
218:   PetscRandomDestroy(&rand);
219:   PetscFinalize();
220:   return 0;
221: }