Actual source code: htool.cxx
1: #include <../src/mat/impls/htool/htool.hpp>
2: #include <set>
4: const char *const MatHtoolCompressorTypes[] = {"sympartialACA", "fullACA", "SVD"};
5: const char *const MatHtoolClusteringTypes[] = {"PCARegular", "PCAGeometric", "BoundingBox1Regular", "BoundingBox1Geometric"};
6: const char *HtoolCitations[2] = {"@article{marchand2020two,\n"
7: " Author = {Marchand, Pierre and Claeys, Xavier and Jolivet, Pierre and Nataf, Fr\\'ed\\'eric and Tournier, Pierre-Henri},\n"
8: " Title = {Two-level preconditioning for $h$-version boundary element approximation of hypersingular operator with {GenEO}},\n"
9: " Year = {2020},\n"
10: " Publisher = {Elsevier},\n"
11: " Journal = {Numerische Mathematik},\n"
12: " Volume = {146},\n"
13: " Pages = {597--628},\n"
14: " Url = {https://github.com/htool-ddm/htool}\n"
15: "}\n",
16: "@article{Marchand2026,\n"
17: " Author = {Marchand, Pierre and Tournier, Pierre-Henri and Jolivet, Pierre},\n"
18: " Title = {{Htool-DDM}: A {C++} library for parallel solvers and compressed linear systems},\n"
19: " Year = {2026},\n"
20: " Publisher = {The Open Journal},\n"
21: " Journal = {Journal of Open Source Software},\n"
22: " Volume = {11},\n"
23: " Number = {118},\n"
24: " Pages = {9279},\n"
25: " Url = {https://doi.org/10.21105/joss.09279}\n"
26: "}\n"};
27: static PetscBool HtoolCite[2] = {PETSC_FALSE, PETSC_FALSE};
29: static PetscErrorCode MatGetDiagonal_Htool(Mat A, Vec v)
30: {
31: Mat_Htool *a;
32: PetscScalar *x;
33: PetscBool flg;
35: PetscFunctionBegin;
36: PetscCall(MatHasCongruentLayouts(A, &flg));
37: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent layouts supported");
38: PetscCall(MatShellGetContext(A, &a));
39: PetscCall(VecGetArrayWrite(v, &x));
40: PetscCallExternalVoid("copy_diagonal_in_user_numbering", htool::copy_diagonal_in_user_numbering(a->distributed_operator_holder->hmatrix, x));
41: PetscCall(VecRestoreArrayWrite(v, &x));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode MatGetDiagonalBlock_Htool(Mat A, Mat *b)
46: {
47: Mat_Htool *a;
48: Mat B;
49: PetscScalar *ptr, shift, scale;
50: PetscBool flg;
51: PetscMPIInt rank;
52: htool::Cluster<PetscReal> *source_cluster = nullptr;
54: PetscFunctionBegin;
55: PetscCall(MatHasCongruentLayouts(A, &flg));
56: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent layouts supported");
57: PetscCall(MatShellGetContext(A, &a));
58: PetscCall(PetscObjectQuery((PetscObject)A, "DiagonalBlock", (PetscObject *)&B)); /* same logic as in MatGetDiagonalBlock_MPIDense() */
59: if (!B) {
60: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
61: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, A->rmap->n, A->rmap->n, A->rmap->n, nullptr, &B));
62: PetscCall(MatDenseGetArrayWrite(B, &ptr));
63: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
64: source_cluster = a->source_cluster ? a->source_cluster.get() : a->target_cluster.get();
65: PetscCallExternalVoid("copy_to_dense_in_user_numbering", htool::copy_to_dense_in_user_numbering(*a->distributed_operator_holder->hmatrix.get_sub_hmatrix(a->target_cluster->get_cluster_on_partition(rank), source_cluster->get_cluster_on_partition(rank)), ptr));
66: PetscCall(MatDenseRestoreArrayWrite(B, &ptr));
67: PetscCall(MatPropagateSymmetryOptions(A, B));
68: PetscCall(PetscObjectCompose((PetscObject)A, "DiagonalBlock", (PetscObject)B));
69: *b = B;
70: PetscCall(MatDestroy(&B));
71: PetscCall(MatShift(*b, shift));
72: PetscCall(MatScale(*b, scale));
73: } else {
74: PetscCall(MatShellGetScalingShifts(A, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
75: *b = B;
76: }
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode MatMult_Htool(Mat A, Vec x, Vec y)
81: {
82: Mat_Htool *a;
83: const PetscScalar *in;
84: PetscScalar *out;
86: PetscFunctionBegin;
87: PetscCall(MatShellGetContext(A, &a));
88: PetscCall(VecGetArrayRead(x, &in));
89: PetscCall(VecGetArrayWrite(y, &out));
90: if (a->permutation == PETSC_TRUE) htool::add_distributed_operator_vector_product_local_to_local<PetscScalar>('N', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, nullptr);
91: else htool::internal_add_distributed_operator_vector_product_local_to_local<PetscScalar>('N', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, nullptr);
92: PetscCall(VecRestoreArrayRead(x, &in));
93: PetscCall(VecRestoreArrayWrite(y, &out));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode MatMultTranspose_Htool(Mat A, Vec x, Vec y)
98: {
99: Mat_Htool *a;
100: const PetscScalar *in;
101: PetscScalar *out;
103: PetscFunctionBegin;
104: PetscCall(MatShellGetContext(A, &a));
105: PetscCall(VecGetArrayRead(x, &in));
106: PetscCall(VecGetArrayWrite(y, &out));
107: if (a->permutation == PETSC_TRUE) htool::add_distributed_operator_vector_product_local_to_local<PetscScalar>('T', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, nullptr);
108: else htool::internal_add_distributed_operator_vector_product_local_to_local<PetscScalar>('T', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, nullptr);
109: PetscCall(VecRestoreArrayRead(x, &in));
110: PetscCall(VecRestoreArrayWrite(y, &out));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatIncreaseOverlap_Htool(Mat A, PetscInt is_max, IS is[], PetscInt ov)
115: {
116: std::set<PetscInt> set;
117: const PetscInt *idx;
118: PetscInt *oidx, size, bs[2];
119: PetscMPIInt csize;
121: PetscFunctionBegin;
122: PetscCall(MatGetBlockSizes(A, bs, bs + 1));
123: if (bs[0] != bs[1]) bs[0] = 1;
124: for (PetscInt i = 0; i < is_max; ++i) {
125: /* basic implementation that adds indices by shifting an IS by -ov, -ov+1..., -1, 1..., ov-1, ov */
126: /* needed to avoid subdomain matrices to replicate A since it is dense */
127: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is[i]), &csize));
128: PetscCheck(csize == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel IS");
129: PetscCall(ISGetSize(is[i], &size));
130: PetscCall(ISGetIndices(is[i], &idx));
131: for (PetscInt j = 0; j < size; ++j) {
132: set.insert(idx[j]);
133: for (PetscInt k = 1; k <= ov; ++k) { /* for each layer of overlap */
134: if (idx[j] - k >= 0) set.insert(idx[j] - k); /* do not insert negative indices */
135: if (idx[j] + k < A->rmap->N && idx[j] + k < A->cmap->N) set.insert(idx[j] + k); /* do not insert indices greater than the dimension of A */
136: }
137: }
138: PetscCall(ISRestoreIndices(is[i], &idx));
139: PetscCall(ISDestroy(is + i));
140: if (bs[0] > 1) {
141: for (std::set<PetscInt>::iterator it = set.cbegin(); it != set.cend(); it++) {
142: std::vector<PetscInt> block(bs[0]);
143: std::iota(block.begin(), block.end(), (*it / bs[0]) * bs[0]);
144: set.insert(block.cbegin(), block.cend());
145: }
146: }
147: size = set.size(); /* size with overlap */
148: PetscCall(PetscMalloc1(size, &oidx));
149: for (const PetscInt j : set) *oidx++ = j;
150: oidx -= size;
151: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, oidx, PETSC_OWN_POINTER, is + i));
152: }
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode MatCreateSubMatrices_Htool(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
157: {
158: Mat_Htool *a;
159: Mat D, B, BT;
160: const PetscScalar *copy;
161: PetscScalar *ptr, shift, scale;
162: const PetscInt *idxr, *idxc, *it;
163: PetscInt nrow, m, i;
164: PetscBool flg;
166: PetscFunctionBegin;
167: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
168: PetscCall(MatShellGetContext(A, &a));
169: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
170: for (i = 0; i < n; ++i) {
171: PetscCall(ISGetLocalSize(irow[i], &nrow));
172: PetscCall(ISGetLocalSize(icol[i], &m));
173: PetscCall(ISGetIndices(irow[i], &idxr));
174: PetscCall(ISGetIndices(icol[i], &idxc));
175: if (scall != MAT_REUSE_MATRIX) PetscCall(MatCreateDense(PETSC_COMM_SELF, nrow, m, nrow, m, nullptr, (*submat) + i));
176: PetscCall(MatDenseGetArrayWrite((*submat)[i], &ptr));
177: if (irow[i] == icol[i]) { /* same row and column IS? */
178: PetscCall(MatHasCongruentLayouts(A, &flg));
179: if (flg) {
180: PetscCall(ISSorted(irow[i], &flg));
181: if (flg) { /* sorted IS? */
182: it = std::lower_bound(idxr, idxr + nrow, A->rmap->rstart);
183: if (it != idxr + nrow && *it == A->rmap->rstart) { /* rmap->rstart in IS? */
184: if (std::distance(idxr, it) + A->rmap->n <= nrow) { /* long enough IS to store the local diagonal block? */
185: for (PetscInt j = 0; j < A->rmap->n && flg; ++j)
186: if (PetscUnlikely(it[j] != A->rmap->rstart + j)) flg = PETSC_FALSE;
187: if (flg) { /* complete local diagonal block in IS? */
188: /* fast extraction when the local diagonal block is part of the submatrix, e.g., for PCASM or PCHPDDM
189: * [ B C E ]
190: * A = [ B D E ]
191: * [ B F E ]
192: */
193: m = std::distance(idxr, it); /* shift of the coefficient (0,0) of block D from above */
194: PetscCall(MatGetDiagonalBlock(A, &D));
195: PetscCall(MatDenseGetArrayRead(D, ©));
196: for (PetscInt k = 0; k < A->rmap->n; ++k) PetscCall(PetscArraycpy(ptr + (m + k) * nrow + m, copy + k * A->rmap->n, A->rmap->n)); /* block D from above */
197: PetscCall(MatDenseRestoreArrayRead(D, ©));
198: if (m) {
199: a->wrapper->copy_submatrix(nrow, m, idxr, idxc, ptr); /* vertical block B from above */
200: /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
201: if (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE) {
202: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, m, A->rmap->n, m, ptr + m, &B));
203: PetscCall(MatDenseSetLDA(B, nrow));
204: PetscCall(MatCreateDense(PETSC_COMM_SELF, m, A->rmap->n, m, A->rmap->n, ptr + m * nrow, &BT));
205: PetscCall(MatDenseSetLDA(BT, nrow));
206: if (A->hermitian == PETSC_BOOL3_TRUE && PetscDefined(USE_COMPLEX)) {
207: PetscCall(MatHermitianTranspose(B, MAT_REUSE_MATRIX, &BT));
208: } else {
209: PetscCall(MatTransposeSetPrecursor(B, BT));
210: PetscCall(MatTranspose(B, MAT_REUSE_MATRIX, &BT));
211: }
212: PetscCall(MatDestroy(&B));
213: PetscCall(MatDestroy(&BT));
214: } else {
215: for (PetscInt k = 0; k < A->rmap->n; ++k) { /* block C from above */
216: a->wrapper->copy_submatrix(m, 1, idxr, idxc + m + k, ptr + (m + k) * nrow);
217: }
218: }
219: }
220: if (m + A->rmap->n != nrow) {
221: a->wrapper->copy_submatrix(nrow, std::distance(it + A->rmap->n, idxr + nrow), idxr, idxc + m + A->rmap->n, ptr + (m + A->rmap->n) * nrow); /* vertical block E from above */
222: /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
223: if (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE) {
224: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, nrow - (m + A->rmap->n), A->rmap->n, nrow - (m + A->rmap->n), ptr + (m + A->rmap->n) * nrow + m, &B));
225: PetscCall(MatDenseSetLDA(B, nrow));
226: PetscCall(MatCreateDense(PETSC_COMM_SELF, nrow - (m + A->rmap->n), A->rmap->n, nrow - (m + A->rmap->n), A->rmap->n, ptr + m * nrow + m + A->rmap->n, &BT));
227: PetscCall(MatDenseSetLDA(BT, nrow));
228: if (A->hermitian == PETSC_BOOL3_TRUE && PetscDefined(USE_COMPLEX)) {
229: PetscCall(MatHermitianTranspose(B, MAT_REUSE_MATRIX, &BT));
230: } else {
231: PetscCall(MatTransposeSetPrecursor(B, BT));
232: PetscCall(MatTranspose(B, MAT_REUSE_MATRIX, &BT));
233: }
234: PetscCall(MatDestroy(&B));
235: PetscCall(MatDestroy(&BT));
236: } else {
237: for (PetscInt k = 0; k < A->rmap->n; ++k) { /* block F from above */
238: a->wrapper->copy_submatrix(std::distance(it + A->rmap->n, idxr + nrow), 1, it + A->rmap->n, idxc + m + k, ptr + (m + k) * nrow + m + A->rmap->n);
239: }
240: }
241: }
242: } /* complete local diagonal block not in IS */
243: } else flg = PETSC_FALSE; /* IS not long enough to store the local diagonal block */
244: } else flg = PETSC_FALSE; /* rmap->rstart not in IS */
245: } /* unsorted IS */
246: }
247: } else flg = PETSC_FALSE; /* different row and column IS */
248: if (!flg) a->wrapper->copy_submatrix(nrow, m, idxr, idxc, ptr); /* reassemble everything */
249: PetscCall(ISRestoreIndices(irow[i], &idxr));
250: PetscCall(ISRestoreIndices(icol[i], &idxc));
251: PetscCall(MatDenseRestoreArrayWrite((*submat)[i], &ptr));
252: PetscCall(MatShift((*submat)[i], shift));
253: PetscCall(MatScale((*submat)[i], scale));
254: }
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode MatDestroy_Htool(Mat A)
259: {
260: Mat_Htool *a;
261: PetscContainer container;
262: MatHtoolKernelTranspose *kernelt;
264: PetscFunctionBegin;
265: PetscCall(MatShellGetContext(A, &a));
266: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", nullptr));
267: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", nullptr));
268: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", nullptr));
269: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", nullptr));
270: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", nullptr));
271: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", nullptr));
272: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", nullptr));
273: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", nullptr));
274: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", nullptr));
275: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUseRecompression_C", nullptr));
276: PetscCall(PetscObjectQuery((PetscObject)A, "KernelTranspose", (PetscObject *)&container));
277: if (container) { /* created in MatTranspose_Htool() */
278: PetscCall(PetscContainerGetPointer(container, &kernelt));
279: PetscCall(MatDestroy(&kernelt->A));
280: PetscCall(PetscObjectCompose((PetscObject)A, "KernelTranspose", nullptr));
281: }
282: if (a->gcoords_source != a->gcoords_target) PetscCall(PetscFree(a->gcoords_source));
283: PetscCall(PetscFree(a->gcoords_target));
284: delete a->wrapper;
285: a->target_cluster.reset();
286: a->source_cluster.reset();
287: a->distributed_operator_holder.reset();
288: PetscCall(PetscFree(a));
289: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", nullptr)); // needed to avoid a call to MatShellSetContext_Immutable()
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode MatView_Htool(Mat A, PetscViewer pv)
294: {
295: Mat_Htool *a;
296: PetscScalar shift, scale;
297: PetscBool flg;
298: std::map<std::string, std::string> hmatrix_information;
300: PetscFunctionBegin;
301: PetscCall(MatShellGetContext(A, &a));
302: hmatrix_information = htool::get_distributed_hmatrix_information(a->distributed_operator_holder->hmatrix, PetscObjectComm((PetscObject)A));
303: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &flg));
304: if (flg) {
305: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
306: PetscCall(PetscViewerASCIIPrintf(pv, "symmetry: %c\n", a->distributed_operator_holder->block_diagonal_hmatrix->get_symmetry()));
307: if (PetscAbsScalar(scale - 1.0) > PETSC_MACHINE_EPSILON) {
308: #if defined(PETSC_USE_COMPLEX)
309: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g+%gi\n", (double)PetscRealPart(scale), (double)PetscImaginaryPart(scale)));
310: #else
311: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g\n", (double)scale));
312: #endif
313: }
314: if (PetscAbsScalar(shift) > PETSC_MACHINE_EPSILON) {
315: #if defined(PETSC_USE_COMPLEX)
316: PetscCall(PetscViewerASCIIPrintf(pv, "shift: %g+%gi\n", (double)PetscRealPart(shift), (double)PetscImaginaryPart(shift)));
317: #else
318: PetscCall(PetscViewerASCIIPrintf(pv, "shift: %g\n", (double)shift));
319: #endif
320: }
321: PetscCall(PetscViewerASCIIPrintf(pv, "maximal cluster leaf size: %" PetscInt_FMT "\n", a->max_cluster_leaf_size));
322: PetscCall(PetscViewerASCIIPrintf(pv, "epsilon: %g\n", (double)a->epsilon));
323: PetscCall(PetscViewerASCIIPrintf(pv, "eta: %g\n", (double)a->eta));
324: PetscCall(PetscViewerASCIIPrintf(pv, "minimum target depth: %" PetscInt_FMT "\n", a->depth[0]));
325: PetscCall(PetscViewerASCIIPrintf(pv, "minimum source depth: %" PetscInt_FMT "\n", a->depth[1]));
326: PetscCall(PetscViewerASCIIPrintf(pv, "compressor: %s\n", MatHtoolCompressorTypes[a->compressor]));
327: PetscCall(PetscViewerASCIIPrintf(pv, "clustering: %s\n", MatHtoolClusteringTypes[a->clustering]));
328: PetscCall(PetscViewerASCIIPrintf(pv, "compression ratio: %s\n", hmatrix_information["Compression_ratio"].c_str()));
329: PetscCall(PetscViewerASCIIPrintf(pv, "space saving: %s\n", hmatrix_information["Space_saving"].c_str()));
330: PetscCall(PetscViewerASCIIPrintf(pv, "block tree consistency: %s\n", PetscBools[a->distributed_operator_holder->hmatrix.is_block_tree_consistent()]));
331: PetscCall(PetscViewerASCIIPrintf(pv, "recompression: %s\n", PetscBools[a->recompression]));
332: PetscCall(PetscViewerASCIIPrintf(pv, "number of dense (resp. low rank) matrices: %s (resp. %s)\n", hmatrix_information["Number_of_dense_blocks"].c_str(), hmatrix_information["Number_of_low_rank_blocks"].c_str()));
333: PetscCall(
334: PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) dense block sizes: (%s, %s, %s)\n", hmatrix_information["Dense_block_size_min"].c_str(), hmatrix_information["Dense_block_size_mean"].c_str(), hmatrix_information["Dense_block_size_max"].c_str()));
335: PetscCall(PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) low rank block sizes: (%s, %s, %s)\n", hmatrix_information["Low_rank_block_size_min"].c_str(), hmatrix_information["Low_rank_block_size_mean"].c_str(),
336: hmatrix_information["Low_rank_block_size_max"].c_str()));
337: PetscCall(PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) ranks: (%s, %s, %s)\n", hmatrix_information["Rank_min"].c_str(), hmatrix_information["Rank_mean"].c_str(), hmatrix_information["Rank_max"].c_str()));
338: }
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */
343: static PetscErrorCode MatGetRow_Htool(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
344: {
345: Mat_Htool *a;
346: PetscScalar shift, scale;
347: PetscInt *idxc;
348: PetscBLASInt one = 1, bn;
350: PetscFunctionBegin;
351: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
352: PetscCall(MatShellGetContext(A, &a));
353: if (nz) *nz = A->cmap->N;
354: if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */
355: PetscCall(PetscMalloc1(A->cmap->N, &idxc));
356: for (PetscInt i = 0; i < A->cmap->N; ++i) idxc[i] = i;
357: }
358: if (idx) *idx = idxc;
359: if (v) {
360: PetscCall(PetscMalloc1(A->cmap->N, v));
361: if (a->wrapper) a->wrapper->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
362: else reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx)->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
363: PetscCall(PetscBLASIntCast(A->cmap->N, &bn));
364: PetscCallCXX(htool::Blas<PetscScalar>::scal(&bn, &scale, *v, &one));
365: if (row < A->cmap->N) (*v)[row] += shift;
366: }
367: if (!idx) PetscCall(PetscFree(idxc));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode MatRestoreRow_Htool(Mat, PetscInt, PetscInt *, PetscInt **idx, PetscScalar **v)
372: {
373: PetscFunctionBegin;
374: if (idx) PetscCall(PetscFree(*idx));
375: if (v) PetscCall(PetscFree(*v));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: static PetscErrorCode MatSetFromOptions_Htool(Mat A, PetscOptionItems PetscOptionsObject)
380: {
381: Mat_Htool *a;
382: PetscInt n;
383: PetscBool flg;
385: PetscFunctionBegin;
386: PetscCall(MatShellGetContext(A, &a));
387: PetscOptionsHeadBegin(PetscOptionsObject, "Htool options");
388: PetscCall(PetscOptionsBoundedInt("-mat_htool_max_cluster_leaf_size", "Maximal leaf size in cluster tree", nullptr, a->max_cluster_leaf_size, &a->max_cluster_leaf_size, nullptr, 0));
389: PetscCall(PetscOptionsBoundedReal("-mat_htool_epsilon", "Relative error in Frobenius norm when approximating a block", nullptr, a->epsilon, &a->epsilon, nullptr, 0.0));
390: PetscCall(PetscOptionsReal("-mat_htool_eta", "Admissibility condition tolerance", nullptr, a->eta, &a->eta, nullptr));
391: PetscCall(PetscOptionsBoundedInt("-mat_htool_min_target_depth", "Minimal cluster tree depth associated with the rows", nullptr, a->depth[0], a->depth, nullptr, 0));
392: PetscCall(PetscOptionsBoundedInt("-mat_htool_min_source_depth", "Minimal cluster tree depth associated with the columns", nullptr, a->depth[1], a->depth + 1, nullptr, 0));
393: PetscCall(PetscOptionsBool("-mat_htool_block_tree_consistency", "Block tree consistency", nullptr, a->block_tree_consistency, &a->block_tree_consistency, nullptr));
394: PetscCall(PetscOptionsBool("-mat_htool_recompression", "Use recompression", nullptr, a->recompression, &a->recompression, nullptr));
396: n = 0;
397: PetscCall(PetscOptionsEList("-mat_htool_compressor", "Type of compression", "MatHtoolCompressorType", MatHtoolCompressorTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes), MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA], &n, &flg));
398: if (flg) a->compressor = MatHtoolCompressorType(n);
399: n = 0;
400: PetscCall(PetscOptionsEList("-mat_htool_clustering", "Type of clustering", "MatHtoolClusteringType", MatHtoolClusteringTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes), MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR], &n, &flg));
401: if (flg) a->clustering = MatHtoolClusteringType(n);
402: PetscOptionsHeadEnd();
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode MatAssemblyEnd_Htool(Mat A, MatAssemblyType)
407: {
408: Mat_Htool *a;
409: const PetscInt *ranges;
410: PetscInt *offset;
411: PetscMPIInt size, rank;
412: char S = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 'H' : (A->symmetric == PETSC_BOOL3_TRUE ? 'S' : 'N'), uplo = S == 'N' ? 'N' : 'U';
413: htool::VirtualGenerator<PetscScalar> *generator = nullptr;
414: htool::ClusterTreeBuilder<PetscReal> recursive_build_strategy;
415: htool::Cluster<PetscReal> *source_cluster;
416: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> compressor;
418: PetscFunctionBegin;
419: for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(HtoolCite); ++i) PetscCall(PetscCitationsRegister(HtoolCitations[i], HtoolCite + i));
420: PetscCall(MatShellGetContext(A, &a));
421: delete a->wrapper;
422: a->target_cluster.reset();
423: a->source_cluster.reset();
424: a->distributed_operator_holder.reset();
425: // clustering
426: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
427: PetscCall(PetscMalloc1(2 * size, &offset));
428: PetscCall(MatGetOwnershipRanges(A, &ranges));
429: for (PetscInt i = 0; i < size; ++i) {
430: offset[2 * i] = ranges[i];
431: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
432: }
433: switch (a->clustering) {
434: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
435: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
436: break;
437: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
438: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
439: break;
440: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
441: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
442: break;
443: default:
444: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
445: }
446: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
447: a->target_cluster = std::make_unique<htool::Cluster<PetscReal>>(recursive_build_strategy.create_cluster_tree_from_local_partition(A->rmap->N, a->dim, a->gcoords_target, 2, size, offset));
448: if (a->gcoords_target != a->gcoords_source) {
449: PetscCall(MatGetOwnershipRangesColumn(A, &ranges));
450: for (PetscInt i = 0; i < size; ++i) {
451: offset[2 * i] = ranges[i];
452: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
453: }
454: switch (a->clustering) {
455: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
456: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
457: break;
458: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
459: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
460: break;
461: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
462: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
463: break;
464: default:
465: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
466: }
467: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
468: a->source_cluster = std::make_unique<htool::Cluster<PetscReal>>(recursive_build_strategy.create_cluster_tree_from_local_partition(A->cmap->N, a->dim, a->gcoords_source, 2, size, offset));
469: S = uplo = 'N';
470: source_cluster = a->source_cluster.get();
471: } else source_cluster = a->target_cluster.get();
472: PetscCall(PetscFree(offset));
473: // generator
474: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
475: else {
476: a->wrapper = nullptr;
477: generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx);
478: }
479: // compressor
480: switch (a->compressor) {
481: case MAT_HTOOL_COMPRESSOR_FULL_ACA:
482: compressor = std::make_shared<htool::fullACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
483: break;
484: case MAT_HTOOL_COMPRESSOR_SVD:
485: compressor = std::make_shared<htool::SVD<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
486: break;
487: default:
488: compressor = std::make_shared<htool::sympartialACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
489: }
490: // local hierarchical matrix
491: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
492: auto hmatrix_builder = htool::HMatrixTreeBuilder<PetscScalar>(a->epsilon, a->eta, S, uplo);
493: if (a->recompression) {
494: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> RecompressedLowRankGenerator = std::make_shared<htool::RecompressedLowRankGenerator<PetscScalar>>(*compressor, std::function<void(htool::LowRankMatrix<PetscScalar> &)>(htool::SVD_recompression<PetscScalar>));
495: hmatrix_builder.set_low_rank_generator(RecompressedLowRankGenerator);
496: } else {
497: hmatrix_builder.set_low_rank_generator(compressor);
498: }
499: hmatrix_builder.set_minimal_target_depth(a->depth[0]);
500: hmatrix_builder.set_minimal_source_depth(a->depth[1]);
501: PetscCheck(a->block_tree_consistency || (!a->block_tree_consistency && !(A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE)), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot have a MatHtool with inconsistent block tree which is either symmetric or Hermitian");
502: hmatrix_builder.set_block_tree_consistency(a->block_tree_consistency);
503: a->distributed_operator_holder = std::make_unique<htool::DefaultApproximationBuilder<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, *a->target_cluster, *source_cluster, hmatrix_builder, PetscObjectComm((PetscObject)A));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: static PetscErrorCode MatProductNumeric_Htool(Mat C)
508: {
509: Mat_Product *product = C->product;
510: Mat_Htool *a;
511: const PetscScalar *in;
512: PetscScalar *out;
513: PetscInt N, lda;
515: PetscFunctionBegin;
516: MatCheckProduct(C, 1);
517: PetscCall(MatGetSize(C, nullptr, &N));
518: PetscCall(MatDenseGetLDA(C, &lda));
519: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
520: PetscCall(MatDenseGetArrayRead(product->B, &in));
521: PetscCall(MatDenseGetArrayWrite(C, &out));
522: PetscCall(MatShellGetContext(product->A, &a));
523: switch (product->type) {
524: case MATPRODUCT_AB:
525: if (a->permutation == PETSC_TRUE) htool::add_distributed_operator_matrix_product_local_to_local<PetscScalar>('N', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, N, nullptr);
526: else htool::internal_add_distributed_operator_matrix_product_local_to_local<PetscScalar>('N', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, N, nullptr);
527: break;
528: case MATPRODUCT_AtB:
529: if (a->permutation == PETSC_TRUE) htool::add_distributed_operator_matrix_product_local_to_local<PetscScalar>('T', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, N, nullptr);
530: else htool::internal_add_distributed_operator_matrix_product_local_to_local<PetscScalar>('T', 1.0, a->distributed_operator_holder->distributed_operator, in, 0.0, out, N, nullptr);
531: break;
532: default:
533: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType %s is not supported", MatProductTypes[product->type]);
534: }
535: PetscCall(MatDenseRestoreArrayWrite(C, &out));
536: PetscCall(MatDenseRestoreArrayRead(product->B, &in));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
541: {
542: Mat_Product *product = C->product;
543: Mat A, B;
544: PetscBool flg;
546: PetscFunctionBegin;
547: MatCheckProduct(C, 1);
548: A = product->A;
549: B = product->B;
550: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, ""));
551: PetscCheck(flg && (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB), PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "ProductType %s not supported for %s", MatProductTypes[product->type], ((PetscObject)product->B)->type_name);
552: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
553: if (product->type == MATPRODUCT_AB) PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
554: else PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
555: }
556: PetscCall(MatSetType(C, MATDENSE));
557: PetscCall(MatSetUp(C));
558: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
559: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
560: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
561: C->ops->productsymbolic = nullptr;
562: C->ops->productnumeric = MatProductNumeric_Htool;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
567: {
568: PetscFunctionBegin;
569: MatCheckProduct(C, 1);
570: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A, const htool::DistributedOperator<PetscScalar> **distributed_operator)
575: {
576: Mat_Htool *a;
578: PetscFunctionBegin;
579: PetscCall(MatShellGetContext(A, &a));
580: *distributed_operator = &a->distributed_operator_holder->distributed_operator;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: /*@C
585: MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a `MATHTOOL`.
587: No Fortran Support, No C Support
589: Input Parameter:
590: . A - hierarchical matrix
592: Output Parameter:
593: . distributed_operator - opaque pointer to a Htool virtual matrix
595: Level: advanced
597: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`
598: @*/
599: PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A, const htool::DistributedOperator<PetscScalar> **distributed_operator)
600: {
601: PetscFunctionBegin;
603: PetscAssertPointer(distributed_operator, 2);
604: PetscTryMethod(A, "MatHtoolGetHierarchicalMat_C", (Mat, const htool::DistributedOperator<PetscScalar> **), (A, distributed_operator));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
609: {
610: Mat_Htool *a;
612: PetscFunctionBegin;
613: PetscCall(MatShellGetContext(A, &a));
614: a->kernel = kernel;
615: a->kernelctx = kernelctx;
616: delete a->wrapper;
617: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /*@C
622: MatHtoolSetKernel - Sets the kernel and context used for the assembly of a `MATHTOOL`.
624: Collective, No Fortran Support
626: Input Parameters:
627: + A - hierarchical matrix
628: . kernel - computational kernel (or `NULL`)
629: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
631: Level: advanced
633: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatCreateHtoolFromKernel()`
634: @*/
635: PetscErrorCode MatHtoolSetKernel(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
636: {
637: PetscFunctionBegin;
640: if (!kernel) PetscAssertPointer(kernelctx, 3);
641: PetscTryMethod(A, "MatHtoolSetKernel_C", (Mat, MatHtoolKernelFn *, void *), (A, kernel, kernelctx));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A, IS *is)
646: {
647: Mat_Htool *a;
648: PetscMPIInt rank;
649: const std::vector<PetscInt> *source;
650: const htool::Cluster<PetscReal> *local_source_cluster;
652: PetscFunctionBegin;
653: PetscCall(MatShellGetContext(A, &a));
654: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
655: local_source_cluster = a->source_cluster ? &a->source_cluster->get_cluster_on_partition(rank) : &a->target_cluster->get_cluster_on_partition(rank);
656: source = &local_source_cluster->get_permutation();
657: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), local_source_cluster->get_size(), source->data() + local_source_cluster->get_offset(), PETSC_COPY_VALUES, is));
658: PetscCall(ISSetPermutation(*is));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster for a `MATHTOOL` matrix.
665: Input Parameter:
666: . A - hierarchical matrix
668: Output Parameter:
669: . is - permutation
671: Level: advanced
673: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationTarget()`, `MatHtoolUsePermutation()`
674: @*/
675: PetscErrorCode MatHtoolGetPermutationSource(Mat A, IS *is)
676: {
677: PetscFunctionBegin;
679: if (!is) PetscAssertPointer(is, 2);
680: PetscTryMethod(A, "MatHtoolGetPermutationSource_C", (Mat, IS *), (A, is));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A, IS *is)
685: {
686: Mat_Htool *a;
687: const std::vector<PetscInt> *target;
688: PetscMPIInt rank;
690: PetscFunctionBegin;
691: PetscCall(MatShellGetContext(A, &a));
692: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
693: target = &a->target_cluster->get_permutation();
694: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), a->target_cluster->get_cluster_on_partition(rank).get_size(), target->data() + a->target_cluster->get_cluster_on_partition(rank).get_offset(), PETSC_COPY_VALUES, is));
695: PetscCall(ISSetPermutation(*is));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster for a `MATHTOOL` matrix.
702: Input Parameter:
703: . A - hierarchical matrix
705: Output Parameter:
706: . is - permutation
708: Level: advanced
710: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolUsePermutation()`
711: @*/
712: PetscErrorCode MatHtoolGetPermutationTarget(Mat A, IS *is)
713: {
714: PetscFunctionBegin;
716: if (!is) PetscAssertPointer(is, 2);
717: PetscTryMethod(A, "MatHtoolGetPermutationTarget_C", (Mat, IS *), (A, is));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A, PetscBool use)
722: {
723: Mat_Htool *a;
725: PetscFunctionBegin;
726: PetscCall(MatShellGetContext(A, &a));
727: a->permutation = use;
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*@
732: MatHtoolUsePermutation - Sets whether a `MATHTOOL` matrix should permute input (resp. output) vectors following its internal source (resp. target) permutation.
734: Input Parameters:
735: + A - hierarchical matrix
736: - use - Boolean value
738: Level: advanced
740: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolGetPermutationTarget()`
741: @*/
742: PetscErrorCode MatHtoolUsePermutation(Mat A, PetscBool use)
743: {
744: PetscFunctionBegin;
747: PetscTryMethod(A, "MatHtoolUsePermutation_C", (Mat, PetscBool), (A, use));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: static PetscErrorCode MatHtoolUseRecompression_Htool(Mat A, PetscBool use)
752: {
753: Mat_Htool *a;
755: PetscFunctionBegin;
756: PetscCall(MatShellGetContext(A, &a));
757: a->recompression = use;
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@
762: MatHtoolUseRecompression - Sets whether a `MATHTOOL` matrix should use recompression.
764: Input Parameters:
765: + A - hierarchical matrix
766: - use - Boolean value
768: Level: advanced
770: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`
771: @*/
772: PetscErrorCode MatHtoolUseRecompression(Mat A, PetscBool use)
773: {
774: PetscFunctionBegin;
777: PetscTryMethod(A, "MatHtoolUseRecompression_C", (Mat, PetscBool), (A, use));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: static PetscErrorCode MatConvert_Htool_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
782: {
783: Mat C;
784: Mat_Htool *a;
785: PetscScalar *array, shift, scale;
786: PetscInt lda;
788: PetscFunctionBegin;
789: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
790: PetscCall(MatShellGetContext(A, &a));
791: if (reuse == MAT_REUSE_MATRIX) {
792: C = *B;
793: PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible dimensions");
794: PetscCall(MatDenseGetLDA(C, &lda));
795: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
796: } else {
797: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
798: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
799: PetscCall(MatSetType(C, MATDENSE));
800: PetscCall(MatSetUp(C));
801: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
802: }
803: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
804: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
805: PetscCall(MatDenseGetArrayWrite(C, &array));
806: htool::copy_to_dense_in_user_numbering(a->distributed_operator_holder->hmatrix, array);
807: PetscCall(MatDenseRestoreArrayWrite(C, &array));
808: PetscCall(MatShift(C, shift));
809: PetscCall(MatScale(C, scale));
810: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &C));
811: else *B = C;
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: static PetscErrorCode GenEntriesTranspose(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *rows, const PetscInt *cols, PetscScalar *ptr, PetscCtx ctx)
816: {
817: MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose *)ctx;
818: PetscScalar *tmp;
820: PetscFunctionBegin;
821: PetscCall(generator->kernel(sdim, N, M, cols, rows, ptr, generator->kernelctx));
822: PetscCall(PetscMalloc1(M * N, &tmp));
823: PetscCall(PetscArraycpy(tmp, ptr, M * N));
824: for (PetscInt i = 0; i < M; ++i) {
825: for (PetscInt j = 0; j < N; ++j) ptr[i + j * M] = tmp[j + i * N];
826: }
827: PetscCall(PetscFree(tmp));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /* naive implementation which keeps a reference to the original Mat */
832: static PetscErrorCode MatTranspose_Htool(Mat A, MatReuse reuse, Mat *B)
833: {
834: Mat C;
835: Mat_Htool *a, *c;
836: PetscScalar shift, scale;
837: PetscInt M = A->rmap->N, N = A->cmap->N, m = A->rmap->n, n = A->cmap->n;
838: PetscContainer container;
839: MatHtoolKernelTranspose *kernelt;
841: PetscFunctionBegin;
842: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
843: PetscCall(MatShellGetContext(A, &a));
844: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
845: PetscCheck(reuse != MAT_INPLACE_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatTranspose() with MAT_INPLACE_MATRIX not supported");
846: if (reuse == MAT_INITIAL_MATRIX) {
847: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
848: PetscCall(MatSetSizes(C, n, m, N, M));
849: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
850: PetscCall(MatSetUp(C));
851: PetscCall(PetscNew(&kernelt));
852: PetscCall(PetscObjectContainerCompose((PetscObject)C, "KernelTranspose", kernelt, PetscCtxDestroyDefault));
853: } else {
854: C = *B;
855: PetscCall(PetscObjectQuery((PetscObject)C, "KernelTranspose", (PetscObject *)&container));
856: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatTranspose() with MAT_INITIAL_MATRIX first");
857: PetscCall(PetscContainerGetPointer(container, &kernelt));
858: }
859: PetscCall(MatShellGetContext(C, &c));
860: c->dim = a->dim;
861: PetscCall(MatShift(C, shift));
862: PetscCall(MatScale(C, scale));
863: c->kernel = GenEntriesTranspose;
864: if (kernelt->A != A) {
865: PetscCall(MatDestroy(&kernelt->A));
866: kernelt->A = A;
867: PetscCall(PetscObjectReference((PetscObject)A));
868: }
869: kernelt->kernel = a->kernel;
870: kernelt->kernelctx = a->kernelctx;
871: c->kernelctx = kernelt;
872: c->max_cluster_leaf_size = a->max_cluster_leaf_size;
873: c->epsilon = a->epsilon;
874: c->eta = a->eta;
875: c->block_tree_consistency = a->block_tree_consistency;
876: c->permutation = a->permutation;
877: c->recompression = a->recompression;
878: c->compressor = a->compressor;
879: c->clustering = a->clustering;
880: if (reuse == MAT_INITIAL_MATRIX) {
881: PetscCall(PetscMalloc1(N * c->dim, &c->gcoords_target));
882: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_source, N * c->dim));
883: if (a->gcoords_target != a->gcoords_source) {
884: PetscCall(PetscMalloc1(M * c->dim, &c->gcoords_source));
885: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_target, M * c->dim));
886: } else c->gcoords_source = c->gcoords_target;
887: }
888: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
889: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
890: if (reuse == MAT_INITIAL_MATRIX) *B = C;
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: static PetscErrorCode MatDestroy_Factor(Mat F)
895: {
896: PetscContainer container;
897: htool::HMatrix<PetscScalar> *A;
899: PetscFunctionBegin;
900: PetscCall(PetscObjectQuery((PetscObject)F, "HMatrix", (PetscObject *)&container));
901: if (container) {
902: PetscCall(PetscContainerGetPointer(container, &A));
903: delete A;
904: PetscCall(PetscObjectCompose((PetscObject)F, "HMatrix", nullptr));
905: }
906: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatFactorGetSolverType_C", nullptr));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: static PetscErrorCode MatFactorGetSolverType_Htool(Mat, MatSolverType *type)
911: {
912: PetscFunctionBegin;
913: *type = MATSOLVERHTOOL;
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: template <char trans>
918: static inline PetscErrorCode MatSolve_Private(Mat A, htool::Matrix<PetscScalar> &X)
919: {
920: PetscContainer container;
921: htool::HMatrix<PetscScalar> *B;
923: PetscFunctionBegin;
924: PetscCheck(A->factortype == MAT_FACTOR_LU || A->factortype == MAT_FACTOR_CHOLESKY, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_UNKNOWN_TYPE, "Only MAT_LU_FACTOR and MAT_CHOLESKY_FACTOR are supported");
925: PetscCall(PetscObjectQuery((PetscObject)A, "HMatrix", (PetscObject *)&container));
926: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call Mat%sFactorNumeric() before Mat%sSolve%s()", A->factortype == MAT_FACTOR_LU ? "LU" : "Cholesky", X.nb_cols() == 1 ? "" : "Mat", trans == 'N' ? "" : "Transpose");
927: PetscCall(PetscContainerGetPointer(container, &B));
928: if (A->factortype == MAT_FACTOR_LU) htool::lu_solve(trans, *B, X);
929: else htool::cholesky_solve('L', *B, X);
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
934: static PetscErrorCode MatSolve_Htool(Mat A, Type b, Type x)
935: {
936: PetscInt n;
937: htool::Matrix<PetscScalar> v;
938: PetscScalar *array;
940: PetscFunctionBegin;
941: PetscCall(VecGetLocalSize(b, &n));
942: PetscCall(VecCopy(b, x));
943: PetscCall(VecGetArrayWrite(x, &array));
944: v.assign(n, 1, array, false);
945: PetscCall(VecRestoreArrayWrite(x, &array));
946: PetscCall(MatSolve_Private<trans>(A, v));
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
951: static PetscErrorCode MatSolve_Htool(Mat A, Type B, Type X)
952: {
953: PetscInt m, N;
954: htool::Matrix<PetscScalar> v;
955: PetscScalar *array;
957: PetscFunctionBegin;
958: PetscCall(MatGetLocalSize(B, &m, nullptr));
959: PetscCall(MatGetLocalSize(B, nullptr, &N));
960: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
961: PetscCall(MatDenseGetArrayWrite(X, &array));
962: v.assign(m, N, array, false);
963: PetscCall(MatDenseRestoreArrayWrite(X, &array));
964: PetscCall(MatSolve_Private<trans>(A, v));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: template <MatFactorType ftype>
969: static PetscErrorCode MatFactorNumeric_Htool(Mat F, Mat A, const MatFactorInfo *)
970: {
971: Mat_Htool *a;
972: htool::HMatrix<PetscScalar> *B;
974: PetscFunctionBegin;
975: PetscCall(MatShellGetContext(A, &a));
976: B = new htool::HMatrix<PetscScalar>(a->distributed_operator_holder->hmatrix);
977: if (ftype == MAT_FACTOR_LU) htool::sequential_lu_factorization(*B);
978: else htool::sequential_cholesky_factorization('L', *B);
979: PetscCall(PetscObjectContainerCompose((PetscObject)F, "HMatrix", B, nullptr));
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: template <MatFactorType ftype>
984: PetscErrorCode MatFactorSymbolic_Htool(Mat F, Mat)
985: {
986: PetscFunctionBegin;
987: F->preallocated = PETSC_TRUE;
988: F->assembled = PETSC_TRUE;
989: F->ops->solve = MatSolve_Htool<'N', Vec>;
990: F->ops->matsolve = MatSolve_Htool<'N', Mat>;
991: if (!PetscDefined(USE_COMPLEX) || ftype == MAT_FACTOR_LU) {
992: F->ops->solvetranspose = MatSolve_Htool<'T', Vec>;
993: F->ops->matsolvetranspose = MatSolve_Htool<'T', Mat>;
994: }
995: F->ops->destroy = MatDestroy_Factor;
996: if (ftype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_LU>;
997: else F->ops->choleskyfactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_CHOLESKY>;
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: static PetscErrorCode MatLUFactorSymbolic_Htool(Mat F, Mat A, IS, IS, const MatFactorInfo *)
1002: {
1003: PetscFunctionBegin;
1004: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_LU>(F, A));
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: static PetscErrorCode MatCholeskyFactorSymbolic_Htool(Mat F, Mat A, IS, const MatFactorInfo *)
1009: {
1010: PetscFunctionBegin;
1011: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_CHOLESKY>(F, A));
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: static PetscErrorCode MatGetFactor_htool_htool(Mat A, MatFactorType ftype, Mat *F)
1016: {
1017: Mat B;
1018: Mat_Htool *a;
1019: PetscMPIInt size;
1021: PetscFunctionBegin;
1022: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1023: PetscCall(MatShellGetContext(A, &a));
1024: PetscCheck(size == 1, PetscObjectComm((PetscObject)A), PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel MatGetFactor()");
1025: PetscCheck(a->block_tree_consistency, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot factor a MatHtool with inconsistent block tree");
1026: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1027: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1028: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &((PetscObject)B)->type_name));
1029: PetscCall(MatSetUp(B));
1031: B->ops->getinfo = MatGetInfo_External;
1032: B->factortype = ftype;
1033: B->trivialsymbolic = PETSC_TRUE;
1035: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_Htool;
1036: else if (ftype == MAT_FACTOR_CHOLESKY) B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_Htool;
1038: PetscCall(PetscFree(B->solvertype));
1039: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &B->solvertype));
1041: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_Htool));
1042: *F = B;
1043: PetscFunctionReturn(PETSC_SUCCESS);
1044: }
1046: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Htool(void)
1047: {
1048: PetscFunctionBegin;
1049: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_LU, MatGetFactor_htool_htool));
1050: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_CHOLESKY, MatGetFactor_htool_htool));
1051: PetscFunctionReturn(PETSC_SUCCESS);
1052: }
1054: /*@C
1055: MatCreateHtoolFromKernel - Creates a `MATHTOOL` from a user-supplied kernel.
1057: Collective, No Fortran Support
1059: Input Parameters:
1060: + comm - MPI communicator
1061: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1062: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1063: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1064: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1065: . spacedim - dimension of the space coordinates
1066: . coords_target - coordinates of the target
1067: . coords_source - coordinates of the source
1068: . kernel - computational kernel (or `NULL`)
1069: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
1071: Output Parameter:
1072: . B - matrix
1074: Options Database Keys:
1075: + -mat_htool_max_cluster_leaf_size <`PetscInt`> - maximal leaf size in cluster tree
1076: . -mat_htool_epsilon <`PetscReal`> - relative error in Frobenius norm when approximating a block
1077: . -mat_htool_eta <`PetscReal`> - admissibility condition tolerance
1078: . -mat_htool_min_target_depth <`PetscInt`> - minimal cluster tree depth associated with the rows
1079: . -mat_htool_min_source_depth <`PetscInt`> - minimal cluster tree depth associated with the columns
1080: . -mat_htool_block_tree_consistency <`PetscBool`> - block tree consistency
1081: . -mat_htool_recompression <`PetscBool`> - use recompression
1082: . -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression
1083: - -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering
1085: Level: intermediate
1087: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATHTOOL`, `PCSetCoordinates()`, `MatHtoolSetKernel()`, `MatHtoolCompressorType`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
1088: @*/
1089: PetscErrorCode MatCreateHtoolFromKernel(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt spacedim, const PetscReal coords_target[], const PetscReal coords_source[], MatHtoolKernelFn *kernel, void *kernelctx, Mat *B)
1090: {
1091: Mat A;
1092: Mat_Htool *a;
1094: PetscFunctionBegin;
1095: PetscCall(MatCreate(comm, &A));
1097: PetscAssertPointer(coords_target, 7);
1098: PetscAssertPointer(coords_source, 8);
1100: if (!kernel) PetscAssertPointer(kernelctx, 10);
1101: PetscCall(MatSetSizes(A, m, n, M, N));
1102: PetscCall(MatSetType(A, MATHTOOL));
1103: PetscCall(MatSetUp(A));
1104: PetscCall(MatShellGetContext(A, &a));
1105: a->dim = spacedim;
1106: a->kernel = kernel;
1107: a->kernelctx = kernelctx;
1108: PetscCall(PetscCalloc1(A->rmap->N * spacedim, &a->gcoords_target));
1109: PetscCall(PetscArraycpy(a->gcoords_target + A->rmap->rstart * spacedim, coords_target, A->rmap->n * spacedim));
1110: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_target, A->rmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global target coordinates */
1111: if (coords_target != coords_source) {
1112: PetscCall(PetscCalloc1(A->cmap->N * spacedim, &a->gcoords_source));
1113: PetscCall(PetscArraycpy(a->gcoords_source + A->cmap->rstart * spacedim, coords_source, A->cmap->n * spacedim));
1114: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_source, A->cmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global source coordinates */
1115: } else a->gcoords_source = a->gcoords_target;
1116: *B = A;
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*MC
1121: MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.
1123: Use `./configure --download-htool` to install PETSc to use Htool.
1125: Options Database Key:
1126: . -mat_type htool - matrix type to `MATHTOOL`
1128: Level: beginner
1130: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
1131: M*/
1132: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
1133: {
1134: Mat_Htool *a;
1136: PetscFunctionBegin;
1137: PetscCall(MatSetType(A, MATSHELL));
1138: PetscCall(PetscNew(&a));
1139: PetscCall(MatShellSetContext(A, a));
1140: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Htool));
1141: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL_BLOCK, (PetscErrorCodeFn *)MatGetDiagonalBlock_Htool));
1142: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Htool));
1143: PetscCall(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1144: if (!PetscDefined(USE_COMPLEX)) PetscCall(MatShellSetOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1145: A->ops->increaseoverlap = MatIncreaseOverlap_Htool;
1146: A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
1147: PetscCall(MatShellSetOperation(A, MATOP_VIEW, (PetscErrorCodeFn *)MatView_Htool));
1148: PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (PetscErrorCodeFn *)MatSetFromOptions_Htool));
1149: PetscCall(MatShellSetOperation(A, MATOP_GET_ROW, (PetscErrorCodeFn *)MatGetRow_Htool));
1150: PetscCall(MatShellSetOperation(A, MATOP_RESTORE_ROW, (PetscErrorCodeFn *)MatRestoreRow_Htool));
1151: PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (PetscErrorCodeFn *)MatAssemblyEnd_Htool));
1152: PetscCall(MatShellSetOperation(A, MATOP_TRANSPOSE, (PetscErrorCodeFn *)MatTranspose_Htool));
1153: PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Htool));
1154: a->dim = 0;
1155: a->gcoords_target = nullptr;
1156: a->gcoords_source = nullptr;
1157: a->max_cluster_leaf_size = 10;
1158: a->epsilon = PetscSqrtReal(PETSC_SMALL);
1159: a->eta = 10.0;
1160: a->depth[0] = 0;
1161: a->depth[1] = 0;
1162: a->block_tree_consistency = PETSC_TRUE;
1163: a->permutation = PETSC_TRUE;
1164: a->recompression = PETSC_FALSE;
1165: a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
1166: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", MatProductSetFromOptions_Htool));
1167: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", MatProductSetFromOptions_Htool));
1168: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", MatConvert_Htool_Dense));
1169: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", MatConvert_Htool_Dense));
1170: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", MatHtoolGetHierarchicalMat_Htool));
1171: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", MatHtoolSetKernel_Htool));
1172: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", MatHtoolGetPermutationSource_Htool));
1173: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", MatHtoolGetPermutationTarget_Htool));
1174: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", MatHtoolUsePermutation_Htool));
1175: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUseRecompression_C", MatHtoolUseRecompression_Htool));
1176: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Immutable));
1177: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
1178: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
1179: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATHTOOL));
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }