2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            A,B,*submatA,*submatB;
  9:   PetscInt       bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize;
 11:   PetscScalar    *vals,rval;
 12:   IS             *is1,*is2;
 13:   PetscRandom    rdm;
 14:   Vec            xx,s1,s2;
 15:   PetscReal      s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
 16:   PetscBool      flg;

 18:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 19:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 20:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 21:   PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
 22:   PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
 23:   M    = m*bs;

 25:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
 26:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 27:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
 28:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 29:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 30:   PetscRandomSetFromOptions(rdm);

 32:   PetscMalloc1(bs,&rows);
 33:   PetscMalloc1(bs,&cols);
 34:   PetscMalloc1(bs*bs,&vals);
 35:   PetscMalloc1(M,&idx);

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

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

 56:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 57:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 58:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 61:   /* Test MatIncreaseOverlap() */
 62:   PetscMalloc1(nd,&is1);
 63:   PetscMalloc1(nd,&is2);


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

 80:   for (i=0; i<nd; ++i) {
 81:     ISEqual(is1[i],is2[i],&flg);
 82:     if (!flg) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%D, flg =%d\n",i,(int)flg);
 83:   }

 85:   for (i=0; i<nd; ++i) {
 86:     ISSort(is1[i]);
 87:     ISSort(is2[i]);
 88:   }

 90:   MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
 91:   MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);

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

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,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
133:       }
134:     }
135:     VecDestroy(&xx);
136:     VecDestroy(&s1);
137:     VecDestroy(&s2);
138:   }

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


161: /*TEST

163:    test:
164:       args: -mat_block_size {{1 2  5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}

166: TEST*/