1 static char help[] = "Tests MATHTOOL with a derived htool::IMatrix<PetscScalar> class\n\n";
2
3 #include <petscmat.h>
4 #include <htool/hmatrix/interfaces/virtual_generator.hpp>
5
GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt * J,const PetscInt * K,PetscScalar * ptr,PetscCtx ctx)6 static PetscErrorCode GenEntries(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, PetscCtx 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(PETSC_SUCCESS);
20 }
21 class MyIMatrix : public htool::VirtualGenerator<PetscScalar> {
22 private:
23 PetscReal *coords;
24 PetscInt sdim;
25
26 public:
MyIMatrix(PetscInt spacedim,PetscReal * gcoords)27 MyIMatrix(PetscInt spacedim, PetscReal *gcoords) : htool::VirtualGenerator<PetscScalar>(), coords(gcoords), sdim(spacedim) { }
28
copy_submatrix(PetscInt M,PetscInt N,const PetscInt * J,const PetscInt * K,PetscScalar * ptr) const29 virtual void copy_submatrix(PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr) const override
30 {
31 PetscReal diff = 0.0;
32
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 };
43
main(int argc,char ** argv)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;
54
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 }
103
104 /*TEST
105
106 build:
107 requires: htool
108 testset:
109 requires: htool
110 output_file: output/empty.out
111 nsize: 4
112 args: -m_local 120 -mat_htool_epsilon 1.0e-2
113 test:
114 args: -symmetric
115 test:
116 args: -mat_htool_block_tree_consistency false
117
118 TEST*/
119