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 = 1.e-10;
 18:   PetscTruth     flg;

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

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

 29:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
 30:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_NULL,&B);
 31:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 32:   PetscRandomSetFromOptions(rdm);

 34:   PetscMalloc(bs*sizeof(PetscInt),&rows);
 35:   PetscMalloc(bs*sizeof(PetscInt),&cols);
 36:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 37:   PetscMalloc(M*sizeof(PetscScalar),&idx);
 38: 
 39:   /* Now set blocks of values */
 40:   for (i=0; i<20*bs; i++) {
 41:       PetscRandomGetValue(rdm,&rval);
 42:       cols[0] = bs*(int)(PetscRealPart(rval)*m);
 43:       PetscRandomGetValue(rdm,&rval);
 44:       rows[0] = bs*(int)(PetscRealPart(rval)*m);
 45:       for (j=1; j<bs; j++) {
 46:         rows[j] = rows[j-1]+1;
 47:         cols[j] = cols[j-1]+1;
 48:       }

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

 58:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 60:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 61:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 63:     /* Test MatIncreaseOverlap() */
 64:   PetscMalloc(nd*sizeof(IS **),&is1);
 65:   PetscMalloc(nd*sizeof(IS **),&is2);

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

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

 87:   for (i=0; i<nd; ++i) {
 88:     ISSort(is1[i]);
 89:     ISSort(is2[i]);
 90:   }
 91: 
 92:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
 93:   MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);

 95:   /* Test MatMult() */
 96:   for (i=0; i<nd; i++) {
 97:     MatGetSize(submatA[i],&mm,&nn);
 98:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
 99:     VecDuplicate(xx,&s1);
100:     VecDuplicate(xx,&s2);
101:     for (j=0; j<3; j++) {
102:       VecSetRandom(xx,rdm);
103:       MatMult(submatA[i],xx,s1);
104:       MatMult(submatB[i],xx,s2);
105:       VecNorm(s1,NORM_2,&s1norm);
106:       VecNorm(s2,NORM_2,&s2norm);
107:       rnorm = s2norm-s1norm;
108:       if (rnorm<-tol || rnorm>tol) {
109:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
110:       }
111:     }
112:     VecDestroy(xx);
113:     VecDestroy(s1);
114:     VecDestroy(s2);
115:   }
116:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
117:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
118:   MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
119: 
120:   /* Test MatMult() */
121:   for (i=0; i<nd; i++) {
122:     MatGetSize(submatA[i],&mm,&nn);
123:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
124:     VecDuplicate(xx,&s1);
125:     VecDuplicate(xx,&s2);
126:     for (j=0; j<3; j++) {
127:       VecSetRandom(xx,rdm);
128:       MatMult(submatA[i],xx,s1);
129:       MatMult(submatB[i],xx,s2);
130:       VecNorm(s1,NORM_2,&s1norm);
131:       VecNorm(s2,NORM_2,&s2norm);
132:       rnorm = s2norm-s1norm;
133:       if (rnorm<-tol || rnorm>tol) {
134:         PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
135:       }
136:     }
137:     VecDestroy(xx);
138:     VecDestroy(s1);
139:     VecDestroy(s2);
140:   }
141: 
142:   /* Free allocated memory */
143:   for (i=0; i<nd; ++i) {
144:     ISDestroy(is1[i]);
145:     ISDestroy(is2[i]);
146:     MatDestroy(submatA[i]);
147:     MatDestroy(submatB[i]);
148:  }
149:   PetscFree(is1);
150:   PetscFree(is2);
151:   PetscFree(idx);
152:   PetscFree(rows);
153:   PetscFree(cols);
154:   PetscFree(vals);
155:   MatDestroy(A);
156:   MatDestroy(B);
157:   PetscFree(submatA);
158:   PetscFree(submatB);
159:   PetscRandomDestroy(rdm);
160:   PetscFinalize();
161:   return 0;
162: }