1 static char help[] = "Tests MATHTOOL with a derived htool::IMatrix<PetscScalar> class\n\n"; 2 3 #include <petscmat.h> 4 #include <htool/misc/petsc.hpp> 5 6 static PetscScalar GenEntry(PetscInt sdim,PetscInt i,PetscInt j,void *ctx) 7 { 8 PetscInt d; 9 PetscReal diff = 0.0,*coords = (PetscReal*)(ctx); 10 11 PetscFunctionBeginUser; 12 for (d = 0; d < sdim; d++) { diff += (coords[i*sdim+d] - coords[j*sdim+d]) * (coords[i*sdim+d] - coords[j*sdim+d]); } 13 PetscFunctionReturn(1.0/(1.0e-2 + PetscSqrtReal(diff))); 14 } 15 16 class MyIMatrix : public htool::IMatrix<PetscScalar> { 17 private: 18 PetscReal *coords; 19 PetscInt sdim; 20 public: 21 MyIMatrix(PetscInt M,PetscInt N,PetscInt spacedim,PetscReal* gcoords) : htool::IMatrix<PetscScalar>(M,N),coords(gcoords),sdim(spacedim) { } 22 23 PetscScalar get_coef(const PetscInt &i, const PetscInt &j) const override 24 { 25 PetscReal diff = 0.0; 26 27 PetscFunctionBeginUser; 28 for (PetscInt d = 0; d < sdim; d++) { diff += (coords[i*sdim+d] - coords[j*sdim+d]) * (coords[i*sdim+d] - coords[j*sdim+d]); } 29 PetscFunctionReturn(1.0/(1.0e-2 + PetscSqrtReal(diff))); 30 } 31 32 void copy_submatrix(PetscInt M, PetscInt N, const PetscInt *const rows, const PetscInt *const cols, PetscScalar *ptr) const override 33 { 34 PetscFunctionBeginUser; 35 for (PetscInt j = 0; j < M; j++) /* could be optimized by the user how they see fit, e.g., vectorization */ 36 for (PetscInt k = 0; k < N; k++) ptr[j+k*M] = this->get_coef(rows[j],cols[k]); 37 PetscFunctionReturnVoid(); 38 } 39 }; 40 41 int main(int argc,char **argv) 42 { 43 Mat A,B,P,R; 44 PetscInt m = 100,dim = 3,M,begin = 0; 45 PetscMPIInt size; 46 PetscReal *coords,*gcoords,norm,epsilon,relative; 47 PetscBool sym = PETSC_FALSE; 48 PetscRandom rdm; 49 MatHtoolKernel kernel = GenEntry; 50 MyIMatrix *imatrix; 51 PetscErrorCode ierr; 52 53 ierr = PetscInitialize(&argc,&argv,(char*)NULL,help);if (ierr) return ierr; 54 ierr = PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);CHKERRQ(ierr); 55 ierr = PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);CHKERRQ(ierr); 56 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);CHKERRQ(ierr); 57 ierr = PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);CHKERRQ(ierr); 58 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 59 M = size*m; 60 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 61 ierr = PetscMalloc1(m*dim,&coords);CHKERRQ(ierr); 62 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr); 63 ierr = PetscRandomGetValuesReal(rdm,m*dim,coords);CHKERRQ(ierr); 64 ierr = PetscCalloc1(M*dim,&gcoords);CHKERRQ(ierr); 65 ierr = MPI_Exscan(&m,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 66 ierr = PetscArraycpy(gcoords+begin*dim,coords,m*dim);CHKERRQ(ierr); 67 ierr = MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 68 imatrix = new MyIMatrix(M,M,dim,gcoords); 69 ierr = MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,NULL,imatrix,&A);CHKERRQ(ierr); /* block-wise assembly using htool::IMatrix<PetscScalar>::copy_submatrix() */ 70 ierr = MatSetOption(A,MAT_SYMMETRIC,sym);CHKERRQ(ierr); 71 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 72 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 73 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 74 ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr); 75 ierr = MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B);CHKERRQ(ierr); /* entry-wise assembly using GenEntry() */ 76 ierr = MatSetOption(B,MAT_SYMMETRIC,sym);CHKERRQ(ierr); 77 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 78 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 79 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 80 ierr = MatViewFromOptions(B,NULL,"-B_view");CHKERRQ(ierr); 81 ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);CHKERRQ(ierr); 82 ierr = MatNorm(P,NORM_FROBENIUS,&relative);CHKERRQ(ierr); 83 ierr = MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr); 84 ierr = MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 85 ierr = MatNorm(R,NORM_INFINITY,&norm);CHKERRQ(ierr); 86 if (PetscAbsReal(norm/relative) > epsilon) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A(!symmetric)-A(symmetric)|| = %g (> %g)",(double)PetscAbsReal(norm/relative),(double)epsilon); 87 ierr = MatDestroy(&B);CHKERRQ(ierr); 88 ierr = MatDestroy(&R);CHKERRQ(ierr); 89 ierr = MatDestroy(&P);CHKERRQ(ierr); 90 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 91 ierr = MatDestroy(&A);CHKERRQ(ierr); 92 ierr = PetscFree(gcoords);CHKERRQ(ierr); 93 ierr = PetscFree(coords);CHKERRQ(ierr); 94 delete imatrix; 95 ierr = PetscFinalize(); 96 return ierr; 97 } 98 99 /*TEST 100 101 build: 102 requires: htool 103 test: 104 requires: htool 105 suffix: 2 106 nsize: 4 107 args: -m_local 120 -mat_htool_epsilon 1.0e-2 -symmetric {{false true}shared output} 108 output_file: output/ex101.out 109 110 TEST*/ 111