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 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); 10 11 PetscFunctionBeginUser; 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 PetscFunctionReturn(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) { } 27 28 void copy_submatrix(PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr) const override 29 { 30 PetscReal diff = 0.0; 31 32 PetscFunctionBeginUser; 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 PetscFunctionReturnVoid(); 40 } 41 }; 42 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; 53 54 PetscFunctionBeginUser; 55 PetscCall(PetscInitialize(&argc,&argv,(char*)NULL,help)); 56 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL)); 57 PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 58 PetscCall(PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL)); 59 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL)); 60 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 61 M = size*m; 62 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 63 PetscCall(PetscMalloc1(m*dim,&coords)); 64 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); 65 PetscCall(PetscRandomGetValuesReal(rdm,m*dim,coords)); 66 PetscCall(PetscCalloc1(M*dim,&gcoords)); 67 PetscCallMPI(MPI_Exscan(&m,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD)); 68 PetscCall(PetscArraycpy(gcoords+begin*dim,coords,m*dim)); 69 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD)); 70 imatrix = new MyIMatrix(M,M,dim,gcoords); 71 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,NULL,imatrix,&A)); /* block-wise assembly using htool::IMatrix<PetscScalar>::copy_submatrix() */ 72 PetscCall(MatSetOption(A,MAT_SYMMETRIC,sym)); 73 PetscCall(MatSetFromOptions(A)); 74 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 75 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 76 PetscCall(MatViewFromOptions(A,NULL,"-A_view")); 77 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B)); /* entry-wise assembly using GenEntries() */ 78 PetscCall(MatSetOption(B,MAT_SYMMETRIC,sym)); 79 PetscCall(MatSetFromOptions(B)); 80 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 81 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 82 PetscCall(MatViewFromOptions(B,NULL,"-B_view")); 83 PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P)); 84 PetscCall(MatNorm(P,NORM_FROBENIUS,&relative)); 85 PetscCall(MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R)); 86 PetscCall(MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN)); 87 PetscCall(MatNorm(R,NORM_INFINITY,&norm)); 88 PetscCheck(PetscAbsReal(norm/relative) <= epsilon,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A(!symmetric)-A(symmetric)|| = %g (> %g)",(double)PetscAbsReal(norm/relative),(double)epsilon); 89 PetscCall(MatDestroy(&B)); 90 PetscCall(MatDestroy(&R)); 91 PetscCall(MatDestroy(&P)); 92 PetscCall(PetscRandomDestroy(&rdm)); 93 PetscCall(MatDestroy(&A)); 94 PetscCall(PetscFree(gcoords)); 95 PetscCall(PetscFree(coords)); 96 delete imatrix; 97 PetscCall(PetscFinalize()); 98 return 0; 99 } 100 101 /*TEST 102 103 build: 104 requires: htool 105 test: 106 requires: htool 107 suffix: 2 108 nsize: 4 109 args: -m_local 120 -mat_htool_epsilon 1.0e-2 -symmetric {{false true}shared output} 110 output_file: output/ex101.out 111 112 TEST*/ 113