Actual source code: ex246.cxx

  1: static char help[] = "Tests MATHTOOL with a derived htool::IMatrix<PetscScalar> class\n\n";

  3: #include <petscmat.h>
  4: #include <htool/misc/petsc.hpp>

  6: static PetscErrorCode GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr,void *ctx)
  7: {
  8:   PetscInt  d,j,k;
  9:   PetscReal diff = 0.0,*coords = (PetscReal*)(ctx);

 12:   for (j = 0; j < M; j++) {
 13:     for (k = 0; k < N; k++) {
 14:       diff = 0.0;
 15:       for (d = 0; d < sdim; d++) diff += (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]) * (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]);
 16:       ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
 17:     }
 18:   }
 19:   return 0;
 20: }
 21: class MyIMatrix : public htool::VirtualGenerator<PetscScalar> {
 22:   private:
 23:   PetscReal *coords;
 24:   PetscInt  sdim;
 25:   public:
 26:   MyIMatrix(PetscInt M,PetscInt N,PetscInt spacedim,PetscReal* gcoords) : htool::VirtualGenerator<PetscScalar>(M,N),coords(gcoords),sdim(spacedim) { }

 28:   void copy_submatrix(PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr) const override
 29:   {
 30:     PetscReal diff = 0.0;

 33:     for (PetscInt j = 0; j < M; j++) /* could be optimized by the user how they see fit, e.g., vectorization */
 34:       for (PetscInt k = 0; k < N; k++) {
 35:         diff = 0.0;
 36:         for (PetscInt d = 0; d < sdim; d++) diff += (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]) * (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]);
 37:         ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
 38:       }
 39:     return;
 40:   }
 41: };

 43: int main(int argc,char **argv)
 44: {
 45:   Mat            A,B,P,R;
 46:   PetscInt       m = 100,dim = 3,M,begin = 0;
 47:   PetscMPIInt    size;
 48:   PetscReal      *coords,*gcoords,norm,epsilon,relative;
 49:   PetscBool      sym = PETSC_FALSE;
 50:   PetscRandom    rdm;
 51:   MatHtoolKernel kernel = GenEntries;
 52:   MyIMatrix      *imatrix;

 54:   PetscInitialize(&argc,&argv,(char*)NULL,help);
 55:   PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);
 56:   PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
 57:   PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);
 58:   PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);
 59:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 60:   M = size*m;
 61:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 62:   PetscMalloc1(m*dim,&coords);
 63:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 64:   PetscRandomGetValuesReal(rdm,m*dim,coords);
 65:   PetscCalloc1(M*dim,&gcoords);
 66:   MPI_Exscan(&m,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
 67:   PetscArraycpy(gcoords+begin*dim,coords,m*dim);
 68:   MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
 69:   imatrix = new MyIMatrix(M,M,dim,gcoords);
 70:   MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,NULL,imatrix,&A); /* block-wise assembly using htool::IMatrix<PetscScalar>::copy_submatrix() */
 71:   MatSetOption(A,MAT_SYMMETRIC,sym);
 72:   MatSetFromOptions(A);
 73:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 74:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 75:   MatViewFromOptions(A,NULL,"-A_view");
 76:   MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B); /* entry-wise assembly using GenEntries() */
 77:   MatSetOption(B,MAT_SYMMETRIC,sym);
 78:   MatSetFromOptions(B);
 79:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 81:   MatViewFromOptions(B,NULL,"-B_view");
 82:   MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);
 83:   MatNorm(P,NORM_FROBENIUS,&relative);
 84:   MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);
 85:   MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);
 86:   MatNorm(R,NORM_INFINITY,&norm);
 88:   MatDestroy(&B);
 89:   MatDestroy(&R);
 90:   MatDestroy(&P);
 91:   PetscRandomDestroy(&rdm);
 92:   MatDestroy(&A);
 93:   PetscFree(gcoords);
 94:   PetscFree(coords);
 95:   delete imatrix;
 96:   PetscFinalize();
 97:   return 0;
 98: }

100: /*TEST

102:    build:
103:       requires: htool
104:    test:
105:       requires: htool
106:       suffix: 2
107:       nsize: 4
108:       args: -m_local 120 -mat_htool_epsilon 1.0e-2 -symmetric {{false true}shared output}
109:       output_file: output/ex101.out

111: TEST*/