Actual source code: ex51.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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*/