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;

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

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

 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:     return;
 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      sym = PETSC_FALSE;
 51:   PetscRandom    rdm;
 52:   MatHtoolKernel kernel = GenEntries;
 53:   MyIMatrix     *imatrix;

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

102: /*TEST

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

113: TEST*/