Actual source code: ex92.c

  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
  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=1,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,M,N,Mbs;
 16:   PetscScalar    *vals,rval,one=1.0;
 17:   IS             *is1,*is2;
 18:   PetscRandom    rand;
 19:   PetscBool      flg;
 20:   PETSC_UNUSED 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_size",&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: 
 33:   MatCreate(PETSC_COMM_WORLD,&A);
 34:   MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
 35:   MatSetType(A,MATBAIJ);
 36:   MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL);
 37:   MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL);

 39:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 40:   PetscRandomSetFromOptions(rand);

 42:   MatGetOwnershipRange(A,&rstart,&rend);
 43:   MatGetSize(A,&M,&N);
 44:   Mbs  = M/bs;

 46:   PetscMalloc(bs*sizeof(PetscInt),&rows);
 47:   PetscMalloc(bs*sizeof(PetscInt),&cols);
 48:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 49:   PetscMalloc(M*sizeof(PetscScalar),&idx);

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

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

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

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

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

112:   for (i=0; i<nd; i++) {
113:     PetscRandomGetValue(rand,&rval);
114:     sz = (PetscInt)((0.5 + 0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
115:     //sz /= (size*nd*10);
116: 
117:     if (rank == vid){
118:       PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
119:     }
120: 
121:     for (j=0; j<sz; j++) {
122:       PetscRandomGetValue(rand,&rval);
123:       idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
124:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
125:     }

127:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
128:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);

130:     if (rank == vid){
131:       printf("[%d] is2[%d]\n",rank,i);
132:       ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
133:     }
134:   }

136:   PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
137:   PetscLogStageRegister("MatOv_BAIJ",&stages[1]);

139:   /* Test MatIncreaseOverlap */
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:   if (rank == vid){
149:     printf("\n[%d] IS from BAIJ:\n",rank);
150:     ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
151:     printf("\n[%d] IS from SBAIJ:\n",rank);
152:     ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
153:   }

155:   for (i=0; i<nd; ++i) {
156:     ISEqual(is1[i],is2[i],&flg);
157:     if (!flg ){
158:       if (rank == 0){
159:         ISSort(is1[i]);
160:         ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
161:         ISSort(is2[i]);
162:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
163:       }
164:       SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
165:     }
166:   }
167: 
168:   /* Test MatGetSubmatrices */
169:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
170:   MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);

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

175:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
176:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
177:   MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
178:   MatMultEqual(A,sA,10,&flg);
179:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");

181:   /* Free allocated memory */
182:   for (i=0; i<nd; ++i) {
183:     ISDestroy(&is1[i]);
184:     ISDestroy(&is2[i]);
185:     MatDestroy(&submatA[i]);
186:     MatDestroy(&submatsA[i]);
187:   }
188:   PetscFree(submatA);
189:   PetscFree(submatsA);
190:   PetscFree(is1);
191:   PetscFree(is2);
192:   PetscFree(idx);
193:   PetscFree(rows);
194:   PetscFree(cols);
195:   PetscFree(vals);
196:   MatDestroy(&A);
197:   MatDestroy(&sA);
198:   PetscRandomDestroy(&rand);
199:   PetscFinalize();
200:   return 0;
201: }