Actual source code: ex54.c


  2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            A,B,*submatA,*submatB;
  9:   PetscInt       bs=1,m=11,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,mm,nn,M,N,Mbs;
 10:   PetscMPIInt    size,rank;
 11:   PetscScalar    *vals,rval;
 12:   IS             *is1,*is2;
 13:   PetscRandom    rdm;
 14:   Vec            xx,s1,s2;
 15:   PetscReal      s1norm,s2norm,rnorm,tol = 100*PETSC_SMALL;
 16:   PetscBool      flg,test_nd0=PETSC_FALSE;

 18:   PetscInitialize(&argc,&args,(char*)0,help);
 19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 20:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 22:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 23:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 24:   PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
 26:   PetscOptionsGetBool(NULL,NULL,"-test_nd0",&test_nd0,NULL);

 28:   /* Create a AIJ matrix A */
 29:   MatCreate(PETSC_COMM_WORLD,&A);
 30:   MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
 31:   MatSetType(A,MATAIJ);
 32:   MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,NULL);
 33:   MatMPIAIJSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
 34:   MatSetFromOptions(A);
 35:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 37:   /* Create a BAIJ matrix B */
 38:   MatCreate(PETSC_COMM_WORLD,&B);
 39:   MatSetSizes(B,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
 40:   MatSetType(B,MATBAIJ);
 41:   MatSeqBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL);
 42:   MatMPIBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
 43:   MatSetFromOptions(B);
 44:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 46:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 47:   PetscRandomSetFromOptions(rdm);

 49:   MatGetOwnershipRange(A,&rstart,&rend);
 50:   MatGetSize(A,&M,&N);
 51:   Mbs  = M/bs;

 53:   PetscMalloc1(bs,&rows);
 54:   PetscMalloc1(bs,&cols);
 55:   PetscMalloc1(bs*bs,&vals);
 56:   PetscMalloc1(M,&idx);

 58:   /* Now set blocks of values */
 59:   for (i=0; i<40*bs; i++) {
 60:     PetscRandomGetValue(rdm,&rval);
 61:     cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
 62:     PetscRandomGetValue(rdm,&rval);
 63:     rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
 64:     for (j=1; j<bs; j++) {
 65:       rows[j] = rows[j-1]+1;
 66:       cols[j] = cols[j-1]+1;
 67:     }

 69:     for (j=0; j<bs*bs; j++) {
 70:       PetscRandomGetValue(rdm,&rval);
 71:       vals[j] = rval;
 72:     }
 73:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 74:     MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
 75:   }

 77:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 78:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 82:   /* Test MatIncreaseOverlap() */
 83:   PetscMalloc1(nd,&is1);
 84:   PetscMalloc1(nd,&is2);

 86:   if (rank == 0 && test_nd0) nd = 0; /* test case */

 88:   for (i=0; i<nd; i++) {
 89:     PetscRandomGetValue(rdm,&rval);
 90:     sz   = (int)(PetscRealPart(rval)*m);
 91:     for (j=0; j<sz; j++) {
 92:       PetscRandomGetValue(rdm,&rval);
 93:       idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
 94:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
 95:     }
 96:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
 97:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
 98:   }
 99:   MatIncreaseOverlap(A,nd,is1,ov);
100:   MatIncreaseOverlap(B,nd,is2,ov);

102:   for (i=0; i<nd; ++i) {
103:     ISEqual(is1[i],is2[i],&flg);

105:     if (!flg) {
106:       PetscPrintf(PETSC_COMM_SELF,"i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n",i,flg,bs,m,ov,nd,size);
107:     }
108:   }

110:   for (i=0; i<nd; ++i) {
111:     ISSort(is1[i]);
112:     ISSort(is2[i]);
113:   }

115:   MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
116:   MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);

118:   /* Test MatMult() */
119:   for (i=0; i<nd; i++) {
120:     MatGetSize(submatA[i],&mm,&nn);
121:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
122:     VecDuplicate(xx,&s1);
123:     VecDuplicate(xx,&s2);
124:     for (j=0; j<3; j++) {
125:       VecSetRandom(xx,rdm);
126:       MatMult(submatA[i],xx,s1);
127:       MatMult(submatB[i],xx,s2);
128:       VecNorm(s1,NORM_2,&s1norm);
129:       VecNorm(s2,NORM_2,&s2norm);
130:       rnorm = s2norm-s1norm;
131:       if (rnorm<-tol || rnorm>tol) {
132:         PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,(double)s1norm,(double)s2norm);
133:       }
134:     }
135:     VecDestroy(&xx);
136:     VecDestroy(&s1);
137:     VecDestroy(&s2);
138:   }

140:   /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
141:   MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
142:   MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);

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

166:   /* Free allocated memory */
167:   for (i=0; i<nd; ++i) {
168:     ISDestroy(&is1[i]);
169:     ISDestroy(&is2[i]);
170:   }
171:   MatDestroySubMatrices(nd,&submatA);
172:   MatDestroySubMatrices(nd,&submatB);

174:   PetscFree(is1);
175:   PetscFree(is2);
176:   PetscFree(idx);
177:   PetscFree(rows);
178:   PetscFree(cols);
179:   PetscFree(vals);
180:   MatDestroy(&A);
181:   MatDestroy(&B);
182:   PetscRandomDestroy(&rdm);
183:   PetscFinalize();
184:   return 0;
185: }

187: /*TEST

189:    test:
190:       nsize: {{1 3}}
191:       args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}} ; done
192:       output_file: output/ex54.out

194:    test:
195:       suffix: 2
196:       args: -nd 2 -test_nd0
197:       output_file: output/ex54.out

199:    test:
200:       suffix: 3
201:       nsize: 3
202:       args: -nd 2 -test_nd0
203:       output_file: output/ex54.out

205: TEST*/