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);
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 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;
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 sym = PETSC_FALSE;
51: PetscRandom rdm;
52: MatHtoolKernel kernel = GenEntries;
53: MyIMatrix *imatrix;
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &argv, (char *)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: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
71: imatrix = new MyIMatrix(M, M, 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(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P));
85: PetscCall(MatNorm(P, NORM_FROBENIUS, &relative));
86: PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R));
87: PetscCall(MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN));
88: PetscCall(MatNorm(R, NORM_INFINITY, &norm));
89: PetscCheck(PetscAbsReal(norm / relative) <= epsilon, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A(!symmetric)-A(symmetric)|| = %g (> %g)", (double)PetscAbsReal(norm / relative), (double)epsilon);
90: PetscCall(MatDestroy(&B));
91: PetscCall(MatDestroy(&R));
92: PetscCall(MatDestroy(&P));
93: PetscCall(PetscRandomDestroy(&rdm));
94: PetscCall(MatDestroy(&A));
95: PetscCall(PetscFree(gcoords));
96: PetscCall(PetscFree(coords));
97: delete imatrix;
98: PetscCall(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*/