Actual source code: ex54.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatBAIJ format.\n";

  4: #include <petscmat.h>

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

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

 25:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 26:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);

 30:   /* MatCreateBAIJ(PETSC_COMM_WORLD,bs,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL,&A); */
 31:   MatCreate(PETSC_COMM_WORLD,&A);
 32:   MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
 33:   MatSetType(A,MATBAIJ);
 34:   MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
 35:   MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
 36:   MatSetFromOptions(A);
 37:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 39:   MatCreateAIJ(PETSC_COMM_WORLD,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL,&B);
 40:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 41:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 42:   PetscRandomSetFromOptions(rdm);

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

 48:   PetscMalloc1(bs,&rows);
 49:   PetscMalloc1(bs,&cols);
 50:   PetscMalloc1(bs*bs,&vals);
 51:   PetscMalloc1(M,&idx);

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

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

 72:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 73:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 74:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 75:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 77:   /* Test MatIncreaseOverlap() */
 78:   PetscMalloc1(nd,&is1);
 79:   PetscMalloc1(nd,&is2);


 82:   for (i=0; i<nd; i++) {
 83:     PetscRandomGetValue(rdm,&rval);
 84:     sz   = (int)(PetscRealPart(rval)*m);
 85:     for (j=0; j<sz; j++) {
 86:       PetscRandomGetValue(rdm,&rval);
 87:       idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
 88:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
 89:     }
 90:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
 91:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
 92:   }
 93:   MatIncreaseOverlap(A,nd,is1,ov);
 94:   MatIncreaseOverlap(B,nd,is2,ov);

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

 99:     if (!flg) {
100:       PetscPrintf(PETSC_COMM_SELF,"i=%D, flg=%d :bs=%D m=%D ov=%D nd=%D np=%D\n",i,flg,bs,m,ov,nd,size);
101:     }
102:   }

104:   for (i=0; i<nd; ++i) {
105:     ISSort(is1[i]);
106:     ISSort(is2[i]);
107:   }

109:   MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
110:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);


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

135:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */

137:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
138:   MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);

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

162:   /* Free allocated memory */
163:   for (i=0; i<nd; ++i) {
164:     ISDestroy(&is1[i]);
165:     ISDestroy(&is2[i]);
166:     MatDestroy(&submatA[i]);
167:     MatDestroy(&submatB[i]);
168:   }
169:   PetscFree(is1);
170:   PetscFree(is2);
171:   PetscFree(idx);
172:   PetscFree(rows);
173:   PetscFree(cols);
174:   PetscFree(vals);
175:   MatDestroy(&A);
176:   MatDestroy(&B);
177:   PetscFree(submatA);
178:   PetscFree(submatB);
179:   PetscRandomDestroy(&rdm);
180:   PetscFinalize();
181:   return 0;
182: }