2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for 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=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize;
 13:   PetscScalar    *vals,rval;
 14:   IS             *is1,*is2;
 15:   PetscRandom    rdm;
 16:   Vec            xx,s1,s2;
 17:   PetscReal      s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
 18:   PetscBool      flg;

 20:   PetscInitialize(&argc,&args,(char*)0,help);


 23:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 24:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
 26:   PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
 27:   M    = m*bs;

 29:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
 30:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 31:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
 32:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 33:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 34:   PetscRandomSetFromOptions(rdm);

 36:   PetscMalloc1(bs,&rows);
 37:   PetscMalloc1(bs,&cols);
 38:   PetscMalloc1(bs*bs,&vals);
 39:   PetscMalloc1(M,&idx);

 41:   /* Now set blocks of values */
 42:   for (i=0; i<20*bs; i++) {
 43:     PetscRandomGetValue(rdm,&rval);
 44:     cols[0] = bs*(int)(PetscRealPart(rval)*m);
 45:     PetscRandomGetValue(rdm,&rval);
 46:     rows[0] = bs*(int)(PetscRealPart(rval)*m);
 47:     for (j=1; j<bs; j++) {
 48:       rows[j] = rows[j-1]+1;
 49:       cols[j] = cols[j-1]+1;
 50:     }

 52:     for (j=0; j<bs*bs; j++) {
 53:       PetscRandomGetValue(rdm,&rval);
 54:       vals[j] = rval;
 55:     }
 56:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 57:     MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
 58:   }

 60:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 61:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 62:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 63:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 65:   /* Test MatIncreaseOverlap() */
 66:   PetscMalloc1(nd,&is1);
 67:   PetscMalloc1(nd,&is2);


 70:   for (i=0; i<nd; i++) {
 71:     PetscRandomGetValue(rdm,&rval);
 72:     lsize = (int)(PetscRealPart(rval)*m);
 73:     for (j=0; j<lsize; j++) {
 74:       PetscRandomGetValue(rdm,&rval);
 75:       idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
 76:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
 77:     }
 78:     ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is1+i);
 79:     ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is2+i);
 80:   }
 81:   MatIncreaseOverlap(A,nd,is1,ov);
 82:   MatIncreaseOverlap(B,nd,is2,ov);

 84:   for (i=0; i<nd; ++i) {
 85:     ISEqual(is1[i],is2[i],&flg);
 86:     PetscPrintf(PETSC_COMM_SELF,"i=%D, flg =%d\n",i,(int)flg);
 87:   }

 89:   for (i=0; i<nd; ++i) {
 90:     ISSort(is1[i]);
 91:     ISSort(is2[i]);
 92:   }

 94:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
 95:   MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);

 97:   /* Test MatMult() */
 98:   for (i=0; i<nd; i++) {
 99:     MatGetSize(submatA[i],&mm,&nn);
100:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
101:     VecDuplicate(xx,&s1);
102:     VecDuplicate(xx,&s2);
103:     for (j=0; j<3; j++) {
104:       VecSetRandom(xx,rdm);
105:       MatMult(submatA[i],xx,s1);
106:       MatMult(submatB[i],xx,s2);
107:       VecNorm(s1,NORM_2,&s1norm);
108:       VecNorm(s2,NORM_2,&s2norm);
109:       rnorm = s2norm-s1norm;
110:       if (rnorm<-tol || rnorm>tol) {
111:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
112:       }
113:     }
114:     VecDestroy(&xx);
115:     VecDestroy(&s1);
116:     VecDestroy(&s2);
117:   }
118:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
119:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
120:   MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);

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

144:   /* Free allocated memory */
145:   for (i=0; i<nd; ++i) {
146:     ISDestroy(&is1[i]);
147:     ISDestroy(&is2[i]);
148:     MatDestroy(&submatA[i]);
149:     MatDestroy(&submatB[i]);
150:   }
151:   PetscFree(is1);
152:   PetscFree(is2);
153:   PetscFree(idx);
154:   PetscFree(rows);
155:   PetscFree(cols);
156:   PetscFree(vals);
157:   MatDestroy(&A);
158:   MatDestroy(&B);
159:   PetscFree(submatA);
160:   PetscFree(submatB);
161:   PetscRandomDestroy(&rdm);
162:   PetscFinalize();
163:   return 0;
164: }