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/hmatrix/interfaces/virtual_generator.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);

 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(PETSC_SUCCESS);
 20: }
 21: class MyIMatrix : public htool::VirtualGenerator<PetscScalar> {
 22: private:
 23:   PetscReal *coords;
 24:   PetscInt   sdim;

 26: public:
 27:   MyIMatrix(PetscInt spacedim, PetscReal *gcoords) : htool::VirtualGenerator<PetscScalar>(), coords(gcoords), sdim(spacedim) { }

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

 33:     PetscFunctionBeginUser;
 34:     for (PetscInt j = 0; j < M; j++) /* could be optimized by the user how they see fit, e.g., vectorization */
 35:       for (PetscInt k = 0; k < N; k++) {
 36:         diff = 0.0;
 37:         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]);
 38:         ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff));
 39:       }
 40:     PetscFunctionReturnVoid();
 41:   }
 42: };

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

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

104: /*TEST

106:    build:
107:       requires: htool
108:    testset:
109:       requires: htool
110:       nsize: 4
111:       args: -m_local 120 -mat_htool_epsilon 1.0e-2
112:       test:
113:         args: -symmetric
114:       test:
115:         args: -mat_htool_block_tree_consistency false

117: TEST*/