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 PetscInt d, j, k; 8 PetscReal diff = 0.0, *coords = (PetscReal *)(ctx); 9 10 PetscFunctionBeginUser; 11 for (j = 0; j < M; j++) { 12 for (k = 0; k < N; k++) { 13 diff = 0.0; 14 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]); 15 ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff)); 16 } 17 } 18 PetscFunctionReturn(0); 19 } 20 class MyIMatrix : public htool::VirtualGenerator<PetscScalar> { 21 private: 22 PetscReal *coords; 23 PetscInt sdim; 24 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 PetscReal diff = 0.0; 30 31 PetscFunctionBeginUser; 32 for (PetscInt j = 0; j < M; j++) /* could be optimized by the user how they see fit, e.g., vectorization */ 33 for (PetscInt k = 0; k < N; k++) { 34 diff = 0.0; 35 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]); 36 ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff)); 37 } 38 PetscFunctionReturnVoid(); 39 } 40 }; 41 42 int main(int argc, char **argv) { 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 = GenEntries; 50 MyIMatrix *imatrix; 51 52 PetscFunctionBeginUser; 53 PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help)); 54 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL)); 55 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 56 PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL)); 57 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL)); 58 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 59 M = size * m; 60 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 61 PetscCall(PetscMalloc1(m * dim, &coords)); 62 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 63 PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords)); 64 PetscCall(PetscCalloc1(M * dim, &gcoords)); 65 PetscCallMPI(MPI_Exscan(&m, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD)); 66 PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim)); 67 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD)); 68 imatrix = new MyIMatrix(M, M, dim, gcoords); 69 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, NULL, imatrix, &A)); /* block-wise assembly using htool::IMatrix<PetscScalar>::copy_submatrix() */ 70 PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym)); 71 PetscCall(MatSetFromOptions(A)); 72 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 73 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 74 PetscCall(MatViewFromOptions(A, NULL, "-A_view")); 75 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &B)); /* entry-wise assembly using GenEntries() */ 76 PetscCall(MatSetOption(B, MAT_SYMMETRIC, sym)); 77 PetscCall(MatSetFromOptions(B)); 78 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 79 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatViewFromOptions(B, NULL, "-B_view")); 81 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P)); 82 PetscCall(MatNorm(P, NORM_FROBENIUS, &relative)); 83 PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R)); 84 PetscCall(MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN)); 85 PetscCall(MatNorm(R, NORM_INFINITY, &norm)); 86 PetscCheck(PetscAbsReal(norm / relative) <= epsilon, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A(!symmetric)-A(symmetric)|| = %g (> %g)", (double)PetscAbsReal(norm / relative), (double)epsilon); 87 PetscCall(MatDestroy(&B)); 88 PetscCall(MatDestroy(&R)); 89 PetscCall(MatDestroy(&P)); 90 PetscCall(PetscRandomDestroy(&rdm)); 91 PetscCall(MatDestroy(&A)); 92 PetscCall(PetscFree(gcoords)); 93 PetscCall(PetscFree(coords)); 94 delete imatrix; 95 PetscCall(PetscFinalize()); 96 return 0; 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