Actual source code: ihtool.cxx
1: #include <../src/mat/impls/htool/htool.hpp>
2: #include <petscdraw.h>
3: #include <set>
5: const char *const MatHtoolCompressorTypes[] = {"sympartialACA", "fullACA", "SVD"};
6: const char *const MatHtoolClusteringTypes[] = {"PCARegular", "PCAGeometric", "BoundingBox1Regular", "BoundingBox1Geometric"};
7: const char *HtoolCitations[2] = {"@article{marchand2020two,\n"
8: " Author = {Marchand, Pierre and Claeys, Xavier and Jolivet, Pierre and Nataf, Fr\\'ed\\'eric and Tournier, Pierre-Henri},\n"
9: " Title = {Two-level preconditioning for $h$-version boundary element approximation of hypersingular operator with {GenEO}},\n"
10: " Year = {2020},\n"
11: " Publisher = {Elsevier},\n"
12: " Journal = {Numerische Mathematik},\n"
13: " Volume = {146},\n"
14: " Pages = {597--628},\n"
15: " Url = {https://github.com/htool-ddm/htool}\n"
16: "}\n",
17: "@article{Marchand2026,\n"
18: " Author = {Marchand, Pierre and Tournier, Pierre-Henri and Jolivet, Pierre},\n"
19: " Title = {{Htool-DDM}: A {C++} library for parallel solvers and compressed linear systems},\n"
20: " Year = {2026},\n"
21: " Publisher = {The Open Journal},\n"
22: " Journal = {Journal of Open Source Software},\n"
23: " Volume = {11},\n"
24: " Number = {118},\n"
25: " Pages = {9279},\n"
26: " Url = {https://doi.org/10.21105/joss.09279}\n"
27: "}\n"};
28: static PetscBool HtoolCite[2] = {PETSC_FALSE, PETSC_FALSE};
30: static PetscErrorCode MatGetDiagonal_Htool(Mat A, Vec v)
31: {
32: Mat_Htool *a;
33: PetscScalar *x;
34: PetscBool flg;
36: PetscFunctionBegin;
37: PetscCall(MatHasCongruentLayouts(A, &flg));
38: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent layouts supported");
39: PetscCall(MatShellGetContext(A, &a));
40: PetscCheck(a->block_diagonal_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Block diagonal htool::HMatrix not found");
41: PetscCall(VecGetArrayWrite(v, &x));
42: PetscCallExternalVoid("copy_diagonal_in_user_numbering", htool::copy_diagonal_in_user_numbering(*a->block_diagonal_hmatrix, x));
43: PetscCall(VecRestoreArrayWrite(v, &x));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode MatGetDiagonalBlock_Htool(Mat A, Mat *b)
48: {
49: Mat_Htool *a, *c;
50: Mat B;
51: PetscScalar shift[2], scale[2];
52: PetscBool flg;
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: 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));
60: if (B) {
61: PetscCall(MatShellGetScalingShifts(B, shift + 1, scale + 1, (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));
62: /* invalidate cache when scale or shift changed; PetscObjectCompose() below releases the old entry */
63: if (scale[0] != scale[1] || shift[0] != shift[1]) B = nullptr;
64: }
65: if (!B) {
66: PetscCheck(a->block_diagonal_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Block diagonal htool::HMatrix not found");
67: PetscCall(MatCreate(PETSC_COMM_SELF, &B));
68: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
69: PetscCall(MatSetType(B, MATHTOOL));
70: PetscCall(MatSetUp(B));
71: PetscCall(MatShellGetContext(B, &c));
72: c->dim = a->dim;
73: c->max_cluster_leaf_size = a->max_cluster_leaf_size;
74: c->epsilon = a->epsilon;
75: c->eta = a->eta;
76: c->depth[0] = a->depth[0];
77: c->depth[1] = a->depth[1];
78: c->block_tree_consistency = a->block_tree_consistency;
79: c->permutation = a->permutation;
80: c->recompression = a->recompression;
81: c->compressor = a->compressor;
82: c->clustering = a->clustering;
83: c->kernel = a->kernel;
84: c->kernelctx = a->kernelctx;
85: c->local_to_local_operator = std::make_unique<htool::LocalToLocalHMatrix<PetscScalar>>(*a->block_diagonal_hmatrix);
86: c->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(a->block_diagonal_hmatrix->get_target_cluster(), a->block_diagonal_hmatrix->get_source_cluster(), PetscObjectComm((PetscObject)A), *c->local_to_local_operator);
87: c->distributed_operator = &c->distributed_operator_holder->distributed_operator;
88: c->block_diagonal_hmatrix = a->block_diagonal_hmatrix;
89: c->local_hmatrix_view = a->block_diagonal_hmatrix;
90: B->assembled = PETSC_TRUE;
91: PetscCall(MatPropagateSymmetryOptions(A, B));
92: PetscCall(PetscObjectCompose((PetscObject)A, "DiagonalBlock", (PetscObject)B));
93: *b = B;
94: PetscCall(MatDestroy(&B));
95: PetscCall(MatScale(*b, *scale));
96: PetscCall(MatShift(*b, *shift));
97: } else *b = B;
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode MatDuplicate_Htool(Mat A, MatDuplicateOption op, Mat *B)
102: {
103: Mat C;
104: Mat_Htool *a, *c;
105: PetscMPIInt rank;
106: PetscScalar shift, scale;
107: htool::Cluster<PetscReal> *source_cluster;
109: PetscFunctionBegin;
110: 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));
111: PetscCall(MatShellGetContext(A, &a));
112: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
113: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
114: PetscCall(MatSetBlockSizesFromMats(C, A, A));
115: PetscCall(MatSetType(C, MATHTOOL));
116: PetscCall(MatSetUp(C));
117: PetscCall(MatPropagateSymmetryOptions(A, C));
118: PetscCall(MatShellGetContext(C, &c));
119: c->dim = a->dim;
120: if (a->gcoords_target) {
121: PetscCall(PetscMalloc1(A->rmap->N * c->dim, &c->gcoords_target));
122: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_target, A->rmap->N * c->dim));
123: }
124: if (a->gcoords_source == a->gcoords_target) c->gcoords_source = c->gcoords_target;
125: else if (a->gcoords_source) {
126: PetscCall(PetscMalloc1(A->cmap->N * c->dim, &c->gcoords_source));
127: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_source, A->cmap->N * c->dim));
128: }
129: c->max_cluster_leaf_size = a->max_cluster_leaf_size;
130: c->epsilon = a->epsilon;
131: c->eta = a->eta;
132: c->depth[0] = a->depth[0];
133: c->depth[1] = a->depth[1];
134: c->block_tree_consistency = a->block_tree_consistency;
135: c->permutation = a->permutation;
136: c->recompression = a->recompression;
137: c->compressor = a->compressor;
138: c->clustering = a->clustering;
139: c->kernel = a->kernel;
140: c->kernelctx = a->kernelctx;
141: // no copy of wrapper because it can be created within MatAssemblyEnd() or useless
142: if (a->local_hmatrix) {
143: c->target_cluster = a->target_cluster;
144: if (a->source_cluster) {
145: c->source_cluster = a->source_cluster;
146: source_cluster = c->source_cluster.get();
147: } else source_cluster = c->target_cluster.get();
148: c->local_hmatrix = std::make_unique<htool::HMatrix<PetscScalar>>(*a->local_hmatrix);
149: if (a->global_to_local_operator) {
150: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
151: c->global_to_local_operator = std::make_unique<htool::RestrictedGlobalToLocalHMatrix<PetscScalar>>(*c->local_hmatrix, c->local_hmatrix->get_target_cluster(), c->local_hmatrix->get_source_cluster(), false, false);
152: c->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(*c->target_cluster, *source_cluster, PetscObjectComm((PetscObject)C), *c->global_to_local_operator);
153: c->distributed_operator = &c->distributed_operator_holder->distributed_operator;
154: c->block_diagonal_hmatrix = c->local_hmatrix->get_sub_hmatrix(c->target_cluster->get_cluster_on_partition(rank), source_cluster->get_cluster_on_partition(rank));
155: c->local_hmatrix_view = c->local_hmatrix.get();
156: } else if (a->local_to_local_operator) {
157: c->local_to_local_operator = std::make_unique<htool::LocalToLocalHMatrix<PetscScalar>>(*a->block_diagonal_hmatrix);
158: c->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(a->block_diagonal_hmatrix->get_target_cluster(), a->block_diagonal_hmatrix->get_source_cluster(), PetscObjectComm((PetscObject)A), *c->local_to_local_operator);
159: c->distributed_operator = &c->distributed_operator_holder->distributed_operator;
160: c->block_diagonal_hmatrix = c->local_hmatrix.get();
161: c->local_hmatrix_view = c->local_hmatrix.get();
162: }
163: C->assembled = PETSC_TRUE;
164: } else if (op != MAT_DO_NOT_COPY_VALUES) {
165: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
166: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
167: }
168: if (C->assembled == PETSC_TRUE) {
169: PetscCall(MatScale(C, scale));
170: PetscCall(MatShift(C, shift));
171: }
172: *B = C;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode MatMult_Htool(Mat A, Vec x, Vec y)
177: {
178: Mat_Htool *a;
179: const PetscScalar *in;
180: PetscScalar *out;
182: PetscFunctionBegin;
183: PetscCall(MatShellGetContext(A, &a));
184: PetscCall(VecGetArrayRead(x, &in));
185: PetscCall(VecGetArrayWrite(y, &out));
186: if (a->permutation) PetscCallExternalVoid("add_distributed_operator_vector_product_local_to_local", htool::add_distributed_operator_vector_product_local_to_local<PetscScalar>('N', 1.0, *a->distributed_operator, in, 0.0, out, nullptr));
187: else PetscCallExternalVoid("internal_add_distributed_operator_vector_product_local_to_local", htool::internal_add_distributed_operator_vector_product_local_to_local<PetscScalar>('N', 1.0, *a->distributed_operator, in, 0.0, out, nullptr));
188: PetscCall(VecRestoreArrayRead(x, &in));
189: PetscCall(VecRestoreArrayWrite(y, &out));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode MatMultTranspose_Htool(Mat A, Vec x, Vec y)
194: {
195: Mat_Htool *a;
196: const PetscScalar *in;
197: PetscScalar *out;
199: PetscFunctionBegin;
200: PetscCall(MatShellGetContext(A, &a));
201: PetscCall(VecGetArrayRead(x, &in));
202: PetscCall(VecGetArrayWrite(y, &out));
203: if (a->permutation) PetscCallExternalVoid("add_distributed_operator_vector_product_local_to_local", htool::add_distributed_operator_vector_product_local_to_local<PetscScalar>('T', 1.0, *a->distributed_operator, in, 0.0, out, nullptr));
204: else PetscCallExternalVoid("internal_add_distributed_operator_vector_product_local_to_local", htool::internal_add_distributed_operator_vector_product_local_to_local<PetscScalar>('T', 1.0, *a->distributed_operator, in, 0.0, out, nullptr));
205: PetscCall(VecRestoreArrayRead(x, &in));
206: PetscCall(VecRestoreArrayWrite(y, &out));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static PetscErrorCode MatIncreaseOverlap_Htool(Mat A, PetscInt is_max, IS is[], PetscInt ov)
211: {
212: std::set<PetscInt> set;
213: const PetscInt *idx;
214: PetscInt *oidx, size, bs[2];
215: PetscMPIInt csize;
217: PetscFunctionBegin;
218: PetscCall(MatGetBlockSizes(A, bs, bs + 1));
219: if (bs[0] != bs[1]) bs[0] = 1;
220: for (PetscInt i = 0; i < is_max; ++i) {
221: /* basic implementation that adds indices by shifting an IS by -ov, -ov+1..., -1, 1..., ov-1, ov */
222: /* needed to avoid subdomain matrices to replicate A since it is dense */
223: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is[i]), &csize));
224: PetscCheck(csize == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel IS");
225: PetscCall(ISGetSize(is[i], &size));
226: PetscCall(ISGetIndices(is[i], &idx));
227: for (PetscInt j = 0; j < size; ++j) {
228: set.insert(idx[j]);
229: for (PetscInt k = 1; k <= ov; ++k) { /* for each layer of overlap */
230: if (idx[j] - k >= 0) set.insert(idx[j] - k); /* do not insert negative indices */
231: 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 */
232: }
233: }
234: PetscCall(ISRestoreIndices(is[i], &idx));
235: PetscCall(ISDestroy(is + i));
236: if (bs[0] > 1) {
237: for (std::set<PetscInt>::iterator it = set.cbegin(); it != set.cend(); it++) {
238: std::vector<PetscInt> block(bs[0]);
239: std::iota(block.begin(), block.end(), (*it / bs[0]) * bs[0]);
240: set.insert(block.cbegin(), block.cend());
241: }
242: }
243: size = set.size(); /* size with overlap */
244: PetscCall(PetscMalloc1(size, &oidx));
245: for (const PetscInt j : set) *oidx++ = j;
246: oidx -= size;
247: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, oidx, PETSC_OWN_POINTER, is + i));
248: }
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode MatCreateSubMatrices_Htool(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
253: {
254: Mat_Htool *a, *d;
255: PetscScalar *ptr;
256: PetscScalar shift, scale;
257: const PetscInt *idxr, *idxc, *it;
258: PetscInt nrow, m;
259: PetscBool flg;
261: PetscFunctionBegin;
262: 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));
263: PetscCall(MatShellGetContext(A, &a));
264: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
265: for (PetscInt i = 0; i < n; ++i) {
266: PetscCall(ISGetLocalSize(irow[i], &nrow));
267: PetscCall(ISGetLocalSize(icol[i], &m));
268: PetscCall(ISGetIndices(irow[i], &idxr));
269: PetscCall(ISGetIndices(icol[i], &idxc));
270: flg = PETSC_FALSE;
271: if (irow[i] == icol[i]) { /* same row and column IS? */
272: PetscCall(MatHasCongruentLayouts(A, &flg));
273: if (flg) {
274: PetscCall(ISSorted(irow[i], &flg));
275: if (flg) { /* sorted IS? */
276: it = std::lower_bound(idxr, idxr + nrow, A->rmap->rstart);
277: if (it != idxr + nrow && *it == A->rmap->rstart) { /* rmap->rstart in IS? */
278: if (std::distance(idxr, it) + A->rmap->n <= nrow) { /* long enough IS to store the local diagonal block? */
279: for (PetscInt j = 0; j < A->rmap->n && flg; ++j)
280: if (PetscUnlikely(it[j] != A->rmap->rstart + j)) flg = PETSC_FALSE;
281: if (flg) { /* complete local diagonal block in IS? */
282: PetscInt nb, nd, na, nblk, didx, bidx = -1, aidx = -1;
283: PetscInt blk_sz[3], blk_off[3];
284: Mat submats[9] = {}, D, B, BT;
285: PetscBool sym = (PetscBool)(A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE);
287: /* fast extraction when the local diagonal block is part of the submatrix, e.g., for PCASM or PCHPDDM:
288: * returns a MATNEST with the local MATHTOOL block D and dense off-diagonal blocks
289: * [ B C E ]
290: * A = [ B D E ]
291: * [ B F E ]
292: */
293: nb = (PetscInt)std::distance(idxr, it); /* size of "before" partition (may be 0) */
294: nd = A->rmap->n; /* size of local diagonal block */
295: na = nrow - nb - nd; /* size of "after" partition (may be 0) */
296: nblk = (nb > 0 ? 1 : 0) + 1 + (na > 0 ? 1 : 0);
297: didx = (nb > 0 ? 1 : 0); /* row/column index of D in the MATNEST */
298: if (nb > 0) bidx = 0;
299: if (na > 0) aidx = nblk - 1;
301: /* block sizes and offsets indexed by MATNEST position */
302: if (nb > 0) {
303: blk_sz[bidx] = nb;
304: blk_off[bidx] = 0;
305: }
306: blk_sz[didx] = nd;
307: blk_off[didx] = nb;
308: if (na > 0) {
309: blk_sz[aidx] = na;
310: blk_off[aidx] = nb + nd;
311: }
313: PetscCall(MatCreate(PETSC_COMM_SELF, &D));
314: PetscCall(MatSetSizes(D, nd, nd, nd, nd));
315: PetscCall(MatSetType(D, MATHTOOL));
316: PetscCall(MatSetUp(D));
317: PetscCall(MatPropagateSymmetryOptions(A, D));
318: PetscCall(MatShellGetContext(D, &d));
319: d->dim = a->dim;
320: if (a->gcoords_target) {
321: PetscCall(PetscMalloc1(A->rmap->N * d->dim, &d->gcoords_target));
322: PetscCall(PetscArraycpy(d->gcoords_target, a->gcoords_target, A->rmap->N * d->dim));
323: }
324: if (a->gcoords_source == a->gcoords_target) d->gcoords_source = d->gcoords_target;
325: else if (a->gcoords_source) {
326: PetscCall(PetscMalloc1(A->cmap->N * d->dim, &d->gcoords_source));
327: PetscCall(PetscArraycpy(d->gcoords_source, a->gcoords_source, A->cmap->N * d->dim));
328: }
329: d->max_cluster_leaf_size = a->max_cluster_leaf_size;
330: d->epsilon = a->epsilon;
331: d->eta = a->eta;
332: d->depth[0] = a->depth[0];
333: d->depth[1] = a->depth[1];
334: d->block_tree_consistency = a->block_tree_consistency;
335: d->permutation = a->permutation;
336: d->recompression = a->recompression;
337: d->compressor = a->compressor;
338: d->clustering = a->clustering;
339: d->kernel = a->kernel;
340: d->kernelctx = a->kernelctx;
341: d->target_cluster = a->target_cluster;
342: if (a->source_cluster) d->source_cluster = a->source_cluster;
343: d->local_hmatrix = std::make_unique<htool::HMatrix<PetscScalar>>(*a->block_diagonal_hmatrix);
344: d->local_to_local_operator = std::make_unique<htool::LocalToLocalHMatrix<PetscScalar>>(*d->local_hmatrix);
345: d->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(a->block_diagonal_hmatrix->get_target_cluster(), a->block_diagonal_hmatrix->get_source_cluster(), PetscObjectComm((PetscObject)A), *d->local_to_local_operator);
346: d->distributed_operator = &d->distributed_operator_holder->distributed_operator;
347: d->block_diagonal_hmatrix = d->local_hmatrix.get();
348: d->local_hmatrix_view = d->local_hmatrix.get();
349: D->assembled = PETSC_TRUE;
351: if (scall != MAT_REUSE_MATRIX) {
352: /* diagonal MATHTOOL block and dense off-diagonal blocks */
353: submats[didx * nblk + didx] = D;
354: PetscCall(PetscObjectReference((PetscObject)D));
355: for (PetscInt kr = 0; kr < nblk; kr++) {
356: for (PetscInt kc = (sym ? kr : 0); kc < nblk; kc++) {
357: if (kr == didx && kc == didx) continue;
358: PetscCall(MatCreateDense(PETSC_COMM_SELF, blk_sz[kr], blk_sz[kc], blk_sz[kr], blk_sz[kc], nullptr, &submats[kr * nblk + kc]));
359: }
360: }
361: PetscCall(MatCreateNest(PETSC_COMM_SELF, nblk, nullptr, nblk, nullptr, submats, (*submat) + i));
362: for (PetscInt k = 0; k < nblk * nblk; k++) PetscCall(MatDestroy(&submats[k]));
363: } else PetscCall(MatNestSetSubMat((*submat)[i], didx, didx, D));
364: PetscCall(MatDestroy(&D));
366: /* fill MATDENSE off-diagonal blocks; upper-triangle (kc >= kr) first, then exploit symmetry */
367: for (PetscInt kr = 0; kr < nblk; kr++) {
368: for (PetscInt kc = (sym ? kr : 0); kc < nblk; kc++) {
369: if (kr == didx && kc == didx) continue; /* MATHTOOL diagonal, skip */
370: PetscCall(MatNestGetSubMat((*submat)[i], kr, kc, &B));
371: PetscCall(MatDenseGetArrayWrite(B, &ptr));
372: a->wrapper->copy_submatrix(blk_sz[kr], blk_sz[kc], idxr + blk_off[kr], idxc + blk_off[kc], ptr);
373: PetscCall(MatDenseRestoreArrayWrite(B, &ptr));
374: }
375: }
376: if (sym && scall == MAT_REUSE_MATRIX) { /* need to reset the lower-triangular blocks to avoid having them being MatScale() and MatShift() twice */
377: for (PetscInt kr = 0; kr < nblk; kr++) {
378: for (PetscInt kc = 0; kc < kr; kc++) PetscCall(MatNestSetSubMat((*submat)[i], kr, kc, nullptr));
379: }
380: }
381: PetscCall(MatScale((*submat)[i], scale));
382: PetscCall(MatShift((*submat)[i], shift)); /* MatScale() and MatShift() before filling the lower-triangular blocks */
383: /* exploit symmetry: lower-triangular blocks are (conjugate) transposes of the upper ones */
384: if (sym) {
385: for (PetscInt kr = 0; kr < nblk; kr++) {
386: for (PetscInt kc = 0; kc < kr; kc++) {
387: PetscCall(MatNestGetSubMat((*submat)[i], kc, kr, &B)); /* upper block (already filled) */
388: if (A->hermitian == PETSC_BOOL3_TRUE && PetscDefined(USE_COMPLEX)) PetscCall(MatCreateHermitianTranspose(B, &BT));
389: else PetscCall(MatCreateTranspose(B, &BT));
390: PetscCall(MatNestSetSubMat((*submat)[i], kr, kc, BT));
391: PetscCall(MatDestroy(&BT));
392: }
393: }
394: }
395: } /* complete local diagonal block in IS */
396: } else flg = PETSC_FALSE; /* IS not long enough to store the local diagonal block */
397: } else flg = PETSC_FALSE; /* rmap->rstart not in IS */
398: } /* unsorted IS */
399: }
400: } else flg = PETSC_FALSE; /* different row and column IS */
401: if (!flg) { /* dense fallback: reassemble everything */
402: if (scall != MAT_REUSE_MATRIX) PetscCall(MatCreateDense(PETSC_COMM_SELF, nrow, m, nrow, m, nullptr, (*submat) + i));
403: PetscCall(MatDenseGetArrayWrite((*submat)[i], &ptr));
404: a->wrapper->copy_submatrix(nrow, m, idxr, idxc, ptr);
405: PetscCall(MatDenseRestoreArrayWrite((*submat)[i], &ptr));
406: }
407: PetscCall(ISRestoreIndices(irow[i], &idxr));
408: PetscCall(ISRestoreIndices(icol[i], &idxc));
409: if (!flg) {
410: PetscCall(MatScale((*submat)[i], scale));
411: PetscCall(MatShift((*submat)[i], shift));
412: }
413: }
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode MatDestroy_Htool(Mat A)
418: {
419: Mat_Htool *a;
420: PetscContainer container;
421: MatHtoolKernelTranspose *kernelt;
423: PetscFunctionBegin;
424: PetscCall(MatShellGetContext(A, &a));
425: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", nullptr));
426: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", nullptr));
427: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", nullptr));
428: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", nullptr));
429: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", nullptr));
430: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", nullptr));
431: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", nullptr));
432: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", nullptr));
433: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", nullptr));
434: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUseRecompression_C", nullptr));
435: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEpsilon_C", nullptr));
436: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEpsilon_C", nullptr));
437: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEta_C", nullptr));
438: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEta_C", nullptr));
439: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMaxClusterLeafSize_C", nullptr));
440: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMaxClusterLeafSize_C", nullptr));
441: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinTargetDepth_C", nullptr));
442: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinTargetDepth_C", nullptr));
443: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinSourceDepth_C", nullptr));
444: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinSourceDepth_C", nullptr));
445: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetBlockTreeConsistency_C", nullptr));
446: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetBlockTreeConsistency_C", nullptr));
447: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetCompressorType_C", nullptr));
448: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetCompressorType_C", nullptr));
449: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetClusteringType_C", nullptr));
450: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetClusteringType_C", nullptr));
451: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolCreateFromKernel_C", nullptr));
452: PetscCall(PetscObjectQuery((PetscObject)A, "KernelTranspose", (PetscObject *)&container));
453: if (container) { /* created in MatTranspose_Htool() */
454: PetscCall(PetscContainerGetPointer(container, &kernelt));
455: PetscCall(MatDestroy(&kernelt->A));
456: PetscCall(PetscObjectCompose((PetscObject)A, "KernelTranspose", nullptr));
457: }
458: if (a->gcoords_source != a->gcoords_target) PetscCall(PetscFree(a->gcoords_source));
459: PetscCall(PetscFree(a->gcoords_target));
460: a->distributed_operator_holder.reset();
461: a->global_to_local_operator.reset();
462: a->local_to_local_operator.reset();
463: a->local_hmatrix.reset();
464: a->source_cluster.reset();
465: a->target_cluster.reset();
466: delete a->wrapper;
467: PetscCall(PetscFree(a));
468: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", nullptr)); // needed to avoid a call to MatShellSetContext_Immutable()
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: static PetscErrorCode MatView_Htool_Draw_Zoom(PetscDraw draw, void *ptr)
473: {
474: Mat A = (Mat)ptr;
475: Mat_Htool *a;
476: PetscReal x_r, y_r, x_l, y_l, w, h;
477: PetscInt min_max[2] = {PETSC_INT_MAX, 0};
478: const int greens[] = {PETSC_DRAW_LIMEGREEN, PETSC_DRAW_FORESTGREEN, PETSC_DRAW_DARKGREEN};
479: int color;
480: char str[16];
481: std::vector<const htool::HMatrix<PetscScalar, PetscReal> *> dense_blocks, low_rank_blocks;
483: PetscFunctionBegin;
484: PetscCall(MatShellGetContext(A, &a));
485: PetscCallExternalVoid("get_leaves", htool::get_leaves(*a->local_hmatrix_view, dense_blocks, low_rank_blocks));
486: for (const htool::HMatrix<PetscScalar, PetscReal> *block : low_rank_blocks) {
487: const PetscInt rank = block->get_rank();
489: if (rank < min_max[0]) min_max[0] = rank;
490: if (rank > min_max[1]) min_max[1] = rank;
491: }
492: PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)A), min_max, min_max));
493: if (min_max[0] == PETSC_INT_MAX) min_max[0] = min_max[1];
494: PetscCall(PetscDrawStringGetSize(draw, &w, &h));
495: PetscDrawCollectiveBegin(draw);
496: for (const htool::HMatrix<PetscScalar, PetscReal> *block : dense_blocks) {
497: x_l = x_r = (PetscReal)block->get_source_cluster().get_offset();
498: x_r += (PetscReal)block->get_source_cluster().get_size();
499: y_l = y_r = (PetscReal)(A->rmap->N - block->get_target_cluster().get_offset());
500: y_l -= (PetscReal)block->get_target_cluster().get_size();
501: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, PETSC_DRAW_RED, PETSC_DRAW_RED, PETSC_DRAW_RED, PETSC_DRAW_RED));
502: PetscCall(PetscDrawLine(draw, x_l, y_l, x_r, y_l, PETSC_DRAW_BLACK));
503: PetscCall(PetscDrawLine(draw, x_r, y_l, x_r, y_r, PETSC_DRAW_BLACK));
504: PetscCall(PetscDrawLine(draw, x_r, y_r, x_l, y_r, PETSC_DRAW_BLACK));
505: PetscCall(PetscDrawLine(draw, x_l, y_r, x_l, y_l, PETSC_DRAW_BLACK));
506: }
507: for (const htool::HMatrix<PetscScalar, PetscReal> *block : low_rank_blocks) {
508: PetscReal th;
509: const PetscInt rank = block->get_rank();
511: x_l = x_r = (PetscReal)block->get_source_cluster().get_offset();
512: x_r += (PetscReal)block->get_source_cluster().get_size();
513: y_l = y_r = (PetscReal)(A->rmap->N - block->get_target_cluster().get_offset());
514: y_l -= (PetscReal)block->get_target_cluster().get_size();
515: if (min_max[1] > min_max[0]) color = greens[(int)((PetscReal)(rank - min_max[0]) / (PetscReal)(min_max[1] - min_max[0]) * (PETSC_STATIC_ARRAY_LENGTH(greens) - 1) + 0.5)];
516: else color = greens[PETSC_STATIC_ARRAY_LENGTH(greens) - 1];
517: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
518: PetscCall(PetscDrawLine(draw, x_l, y_l, x_r, y_l, PETSC_DRAW_BLACK));
519: PetscCall(PetscDrawLine(draw, x_r, y_l, x_r, y_r, PETSC_DRAW_BLACK));
520: PetscCall(PetscDrawLine(draw, x_r, y_r, x_l, y_r, PETSC_DRAW_BLACK));
521: PetscCall(PetscDrawLine(draw, x_l, y_r, x_l, y_l, PETSC_DRAW_BLACK));
522: PetscCall(PetscSNPrintf(str, sizeof(str), "%d", rank));
523: PetscCall(PetscDrawStringGetSize(draw, nullptr, &th));
524: if (x_r - x_l > 4 * w && y_r - y_l > 4 * h) PetscCall(PetscDrawStringCentered(draw, 0.5 * (x_l + x_r), 0.5 * (y_l + y_r) - th / 2, PETSC_DRAW_BLACK, str));
525: }
526: PetscDrawCollectiveEnd(draw);
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: static PetscErrorCode MatView_Htool_Draw(Mat A, PetscViewer viewer)
531: {
532: PetscDraw draw;
533: PetscReal x_r = (PetscReal)A->cmap->N, y_r = (PetscReal)A->rmap->N, w, h;
534: PetscBool flg;
536: PetscFunctionBegin;
537: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
538: PetscCall(PetscDrawIsNull(draw, &flg));
539: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
540: w = x_r / 10.0;
541: h = y_r / 10.0;
542: PetscCall(PetscDrawSetCoordinates(draw, -w, -h, x_r + w, y_r + h));
543: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
544: PetscCall(PetscDrawZoom(draw, MatView_Htool_Draw_Zoom, A));
545: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", nullptr));
546: PetscCall(PetscDrawSave(draw));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode MatView_Htool(Mat A, PetscViewer pv)
551: {
552: Mat_Htool *a;
553: PetscScalar shift, scale;
554: PetscBool flg;
555: std::map<std::string, std::string> hmatrix_information;
557: PetscFunctionBegin;
558: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERDRAW, &flg));
559: if (flg) PetscCall(MatView_Htool_Draw(A, pv));
560: else {
561: PetscCall(MatShellGetContext(A, &a));
562: hmatrix_information = htool::get_distributed_hmatrix_information(*a->local_hmatrix_view, PetscObjectComm((PetscObject)A));
563: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &flg));
564: if (flg) {
565: 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));
566: PetscCall(PetscViewerASCIIPrintf(pv, "symmetry: %c\n", a->block_diagonal_hmatrix ? a->block_diagonal_hmatrix->get_symmetry() : 'N'));
567: if (PetscAbsScalar(scale - 1.0) > PETSC_MACHINE_EPSILON) {
568: #if PetscDefined(USE_COMPLEX)
569: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g+%gi\n", (double)PetscRealPart(scale), (double)PetscImaginaryPart(scale)));
570: #else
571: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g\n", (double)scale));
572: #endif
573: }
574: if (PetscAbsScalar(shift) > PETSC_MACHINE_EPSILON) {
575: #if PetscDefined(USE_COMPLEX)
576: PetscCall(PetscViewerASCIIPrintf(pv, "shift: %g+%gi\n", (double)PetscRealPart(shift), (double)PetscImaginaryPart(shift)));
577: #else
578: PetscCall(PetscViewerASCIIPrintf(pv, "shift: %g\n", (double)shift));
579: #endif
580: }
581: PetscCall(PetscViewerASCIIPrintf(pv, "maximal cluster leaf size: %" PetscInt_FMT "\n", a->max_cluster_leaf_size));
582: PetscCall(PetscViewerASCIIPrintf(pv, "epsilon: %g\n", (double)a->epsilon));
583: PetscCall(PetscViewerASCIIPrintf(pv, "eta: %g\n", (double)a->eta));
584: PetscCall(PetscViewerASCIIPrintf(pv, "minimum target depth: %" PetscInt_FMT "\n", a->depth[0]));
585: PetscCall(PetscViewerASCIIPrintf(pv, "minimum source depth: %" PetscInt_FMT "\n", a->depth[1]));
586: PetscCall(PetscViewerASCIIPrintf(pv, "compressor: %s\n", MatHtoolCompressorTypes[a->compressor]));
587: PetscCall(PetscViewerASCIIPrintf(pv, "clustering: %s\n", MatHtoolClusteringTypes[a->clustering]));
588: PetscCall(PetscViewerASCIIPrintf(pv, "compression ratio: %s\n", hmatrix_information["Compression_ratio"].c_str()));
589: PetscCall(PetscViewerASCIIPrintf(pv, "space saving: %s\n", hmatrix_information["Space_saving"].c_str()));
590: PetscCall(PetscViewerASCIIPrintf(pv, "block tree consistency: %s\n", PetscBools[a->local_hmatrix_view->is_block_tree_consistent()]));
591: PetscCall(PetscViewerASCIIPrintf(pv, "recompression: %s\n", PetscBools[a->recompression]));
592: 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()));
593: PetscCall(
594: 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()));
595: 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(),
596: hmatrix_information["Low_rank_block_size_max"].c_str()));
597: 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()));
598: }
599: }
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */
604: static PetscErrorCode MatGetRow_Htool(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
605: {
606: Mat_Htool *a;
607: PetscScalar shift, scale;
608: PetscInt *idxc;
609: PetscBLASInt one = 1, bn;
611: PetscFunctionBegin;
612: 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));
613: PetscCall(MatShellGetContext(A, &a));
614: if (nz) *nz = A->cmap->N;
615: if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */
616: PetscCall(PetscMalloc1(A->cmap->N, &idxc));
617: for (PetscInt i = 0; i < A->cmap->N; ++i) idxc[i] = i;
618: }
619: if (idx) *idx = idxc;
620: if (v) {
621: PetscCall(PetscMalloc1(A->cmap->N, v));
622: if (a->wrapper) a->wrapper->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
623: else reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx)->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
624: PetscCall(PetscBLASIntCast(A->cmap->N, &bn));
625: PetscCallExternalVoid("scal", htool::Blas<PetscScalar>::scal(&bn, &scale, *v, &one));
626: if (row < A->cmap->N) (*v)[row] += shift;
627: }
628: if (!idx) PetscCall(PetscFree(idxc));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode MatRestoreRow_Htool(Mat, PetscInt, PetscInt *, PetscInt **idx, PetscScalar **v)
633: {
634: PetscFunctionBegin;
635: if (idx) PetscCall(PetscFree(*idx));
636: if (v) PetscCall(PetscFree(*v));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: static PetscErrorCode MatSetFromOptions_Htool(Mat A, PetscOptionItems PetscOptionsObject)
641: {
642: Mat_Htool *a;
643: PetscReal r;
644: PetscInt n;
645: PetscBool b, flg, changed = PETSC_FALSE;
647: PetscFunctionBegin;
648: PetscCall(MatShellGetContext(A, &a));
649: PetscOptionsHeadBegin(PetscOptionsObject, "Htool options");
650: n = a->max_cluster_leaf_size;
651: PetscCall(PetscOptionsBoundedInt("-mat_htool_max_cluster_leaf_size", "Maximal leaf size in cluster tree", nullptr, a->max_cluster_leaf_size, &n, &flg, 0));
652: if (flg) {
653: if (n != a->max_cluster_leaf_size) changed = PETSC_TRUE;
654: a->max_cluster_leaf_size = n;
655: }
656: r = a->epsilon;
657: PetscCall(PetscOptionsBoundedReal("-mat_htool_epsilon", "Relative error in Frobenius norm when approximating a block", nullptr, a->epsilon, &r, &flg, 0.0));
658: if (flg) {
659: if (r != a->epsilon) changed = PETSC_TRUE;
660: a->epsilon = r;
661: }
662: r = a->eta;
663: PetscCall(PetscOptionsReal("-mat_htool_eta", "Admissibility condition tolerance", nullptr, a->eta, &r, &flg));
664: if (flg) {
665: if (r != a->eta) changed = PETSC_TRUE;
666: a->eta = r;
667: }
668: n = a->depth[0];
669: PetscCall(PetscOptionsBoundedInt("-mat_htool_min_target_depth", "Minimal cluster tree depth associated with the rows", nullptr, a->depth[0], &n, &flg, 0));
670: if (flg) {
671: if (n != a->depth[0]) changed = PETSC_TRUE;
672: a->depth[0] = n;
673: }
674: n = a->depth[1];
675: PetscCall(PetscOptionsBoundedInt("-mat_htool_min_source_depth", "Minimal cluster tree depth associated with the columns", nullptr, a->depth[1], &n, &flg, 0));
676: if (flg) {
677: if (n != a->depth[1]) changed = PETSC_TRUE;
678: a->depth[1] = n;
679: }
680: b = a->block_tree_consistency;
681: PetscCall(PetscOptionsBool("-mat_htool_block_tree_consistency", "Block tree consistency", nullptr, a->block_tree_consistency, &b, &flg));
682: if (flg) {
683: if (b != a->block_tree_consistency) changed = PETSC_TRUE;
684: a->block_tree_consistency = b;
685: }
686: b = a->recompression;
687: PetscCall(PetscOptionsBool("-mat_htool_recompression", "Use recompression", nullptr, a->recompression, &b, &flg));
688: if (flg) {
689: if (b != a->recompression) changed = PETSC_TRUE;
690: a->recompression = b;
691: }
692: n = static_cast<PetscInt>(a->compressor);
693: PetscCall(PetscOptionsEList("-mat_htool_compressor", "Type of compression", "MatHtoolCompressorType", MatHtoolCompressorTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes), MatHtoolCompressorTypes[a->compressor], &n, &flg));
694: if (flg) {
695: if (n != static_cast<PetscInt>(a->compressor)) changed = PETSC_TRUE;
696: a->compressor = MatHtoolCompressorType(n);
697: }
698: n = static_cast<PetscInt>(a->clustering);
699: PetscCall(PetscOptionsEList("-mat_htool_clustering", "Type of clustering", "MatHtoolClusteringType", MatHtoolClusteringTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes), MatHtoolClusteringTypes[a->clustering], &n, &flg));
700: if (flg) {
701: if (n != static_cast<PetscInt>(a->clustering)) changed = PETSC_TRUE;
702: a->clustering = MatHtoolClusteringType(n);
703: }
704: PetscOptionsHeadEnd();
705: if (changed) A->assembled = PETSC_FALSE;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /* Htool symmetry tag matching the current flags of A: 'H' (Hermitian, complex only), 'S' (symmetric), or 'N' */
710: static inline char MatHtoolGetSymmetry_Private(Mat A)
711: {
712: return PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 'H' : (A->symmetric == PETSC_BOOL3_TRUE ? 'S' : 'N');
713: }
715: static PetscErrorCode MatAssemblyEnd_Htool(Mat A, MatAssemblyType)
716: {
717: Mat_Htool *a;
718: const PetscInt *ranges;
719: PetscInt *offset;
720: PetscMPIInt size, rank;
721: char S = MatHtoolGetSymmetry_Private(A), uplo = S == 'N' ? 'N' : 'U';
722: htool::VirtualGenerator<PetscScalar> *generator = nullptr;
723: htool::ClusterTreeBuilder<PetscReal> recursive_build_strategy;
724: htool::Cluster<PetscReal> *source_cluster;
725: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> compressor;
727: PetscFunctionBegin;
728: for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(HtoolCite); ++i) PetscCall(PetscCitationsRegister(HtoolCitations[i], HtoolCite + i));
729: PetscCall(MatShellGetContext(A, &a));
730: if (A->was_assembled != PETSC_TRUE) {
731: delete a->wrapper;
732: a->target_cluster.reset();
733: a->source_cluster.reset();
734: a->local_hmatrix.reset();
735: a->local_to_local_operator.reset();
736: a->global_to_local_operator.reset();
737: a->distributed_operator_holder.reset();
738: // clustering
739: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
740: PetscCall(PetscMalloc1(2 * size, &offset));
741: PetscCall(MatGetOwnershipRanges(A, &ranges));
742: for (PetscInt i = 0; i < size; ++i) {
743: offset[2 * i] = ranges[i];
744: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
745: }
746: switch (a->clustering) {
747: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
748: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
749: break;
750: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
751: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
752: break;
753: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
754: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
755: break;
756: default:
757: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
758: }
759: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
760: a->target_cluster = std::make_shared<htool::Cluster<PetscReal>>(recursive_build_strategy.create_cluster_tree_from_local_partition(A->rmap->N, a->dim, a->gcoords_target, 2, size, offset));
761: if (a->gcoords_target != a->gcoords_source) {
762: PetscCall(MatGetOwnershipRangesColumn(A, &ranges));
763: for (PetscInt i = 0; i < size; ++i) {
764: offset[2 * i] = ranges[i];
765: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
766: }
767: switch (a->clustering) {
768: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
769: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
770: break;
771: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
772: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
773: break;
774: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
775: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
776: break;
777: default:
778: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
779: }
780: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
781: a->source_cluster = std::make_shared<htool::Cluster<PetscReal>>(recursive_build_strategy.create_cluster_tree_from_local_partition(A->cmap->N, a->dim, a->gcoords_source, 2, size, offset));
782: S = uplo = 'N';
783: source_cluster = a->source_cluster.get();
784: } else source_cluster = a->target_cluster.get();
785: PetscCall(PetscFree(offset));
786: // generator
787: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
788: else {
789: a->wrapper = nullptr;
790: generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx);
791: }
792: // compressor
793: switch (a->compressor) {
794: case MAT_HTOOL_COMPRESSOR_FULL_ACA:
795: compressor = std::make_shared<htool::fullACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
796: break;
797: case MAT_HTOOL_COMPRESSOR_SVD:
798: compressor = std::make_shared<htool::SVD<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
799: break;
800: default:
801: compressor = std::make_shared<htool::sympartialACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
802: }
803: // local hierarchical matrix
804: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
805: auto hmatrix_builder = htool::HMatrixTreeBuilder<PetscScalar>(a->epsilon, a->eta, S, uplo);
806: if (a->recompression) {
807: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> RecompressedLowRankGenerator = std::make_shared<htool::RecompressedLowRankGenerator<PetscScalar>>(*compressor, std::function<void(htool::LowRankMatrix<PetscScalar> &)>(htool::SVD_recompression<PetscScalar>));
808: hmatrix_builder.set_low_rank_generator(RecompressedLowRankGenerator);
809: } else hmatrix_builder.set_low_rank_generator(compressor);
810: hmatrix_builder.set_minimal_target_depth(a->depth[0]);
811: hmatrix_builder.set_minimal_source_depth(a->depth[1]);
812: 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");
813: hmatrix_builder.set_block_tree_consistency(a->block_tree_consistency);
814: a->local_hmatrix = std::make_unique<htool::HMatrix<PetscScalar>>(hmatrix_builder.build(a->wrapper ? *a->wrapper : *generator, *a->target_cluster, *source_cluster, rank, rank));
815: a->global_to_local_operator = std::make_unique<htool::RestrictedGlobalToLocalHMatrix<PetscScalar>>(*a->local_hmatrix, a->local_hmatrix->get_target_cluster(), a->local_hmatrix->get_source_cluster(), false, false);
816: a->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(*a->target_cluster, *source_cluster, PetscObjectComm((PetscObject)A), *a->global_to_local_operator);
817: a->distributed_operator = &a->distributed_operator_holder->distributed_operator;
818: a->block_diagonal_hmatrix = a->local_hmatrix->get_sub_hmatrix(a->target_cluster->get_cluster_on_partition(rank), source_cluster->get_cluster_on_partition(rank));
819: a->local_hmatrix_view = a->local_hmatrix.get();
820: }
821: A->was_assembled = PETSC_FALSE;
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode MatProductNumeric_Htool(Mat C)
826: {
827: Mat_Product *product = C->product;
828: Mat_Htool *a;
829: const PetscScalar *in;
830: PetscScalar *out;
831: PetscInt N, lda;
833: PetscFunctionBegin;
834: MatCheckProduct(C, 1);
835: PetscCall(MatGetSize(C, nullptr, &N));
836: PetscCall(MatDenseGetLDA(C, &lda));
837: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
838: PetscCall(MatDenseGetArrayRead(product->B, &in));
839: PetscCall(MatDenseGetArrayWrite(C, &out));
840: PetscCall(MatShellGetContext(product->A, &a));
841: switch (product->type) {
842: case MATPRODUCT_AB:
843: if (a->permutation) PetscCallExternalVoid("add_distributed_operator_matrix_product_local_to_local", htool::add_distributed_operator_matrix_product_local_to_local<PetscScalar>('N', 1.0, *a->distributed_operator, in, 0.0, out, N, nullptr));
844: else PetscCallExternalVoid("internal_add_distributed_operator_matrix_product_local_to_local", htool::internal_add_distributed_operator_matrix_product_local_to_local<PetscScalar>('N', 1.0, *a->distributed_operator, in, 0.0, out, N, nullptr));
845: break;
846: case MATPRODUCT_AtB:
847: if (a->permutation) PetscCallExternalVoid("add_distributed_operator_matrix_product_local_to_local", htool::add_distributed_operator_matrix_product_local_to_local<PetscScalar>('T', 1.0, *a->distributed_operator, in, 0.0, out, N, nullptr));
848: else PetscCallExternalVoid("internal_add_distributed_operator_matrix_product_local_to_local", htool::internal_add_distributed_operator_matrix_product_local_to_local<PetscScalar>('T', 1.0, *a->distributed_operator, in, 0.0, out, N, nullptr));
849: break;
850: default:
851: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType %s is not supported", MatProductTypes[product->type]);
852: }
853: PetscCall(MatDenseRestoreArrayWrite(C, &out));
854: PetscCall(MatDenseRestoreArrayRead(product->B, &in));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
859: {
860: Mat_Product *product = C->product;
861: Mat A, B;
862: PetscBool flg;
864: PetscFunctionBegin;
865: MatCheckProduct(C, 1);
866: A = product->A;
867: B = product->B;
868: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, ""));
869: 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);
870: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
871: if (product->type == MATPRODUCT_AB) PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
872: else PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
873: }
874: PetscCall(MatSetType(C, MATDENSE));
875: PetscCall(MatSetVecType(C, B->defaultvectype));
876: PetscCall(MatSetUp(C));
877: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
878: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
879: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
880: C->ops->productsymbolic = nullptr;
881: C->ops->productnumeric = MatProductNumeric_Htool;
882: PetscFunctionReturn(PETSC_SUCCESS);
883: }
885: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
886: {
887: PetscFunctionBegin;
888: MatCheckProduct(C, 1);
889: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A, void *distributed_operator)
894: {
895: Mat_Htool *a;
897: PetscFunctionBegin;
898: PetscCall(MatShellGetContext(A, &a));
899: *(const void **)distributed_operator = static_cast<const void *>(a->distributed_operator);
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
904: {
905: Mat_Htool *a;
907: PetscFunctionBegin;
908: PetscCall(MatShellGetContext(A, &a));
909: if (kernel != a->kernel || a->kernelctx != kernelctx) A->assembled = PETSC_FALSE;
910: a->kernel = kernel;
911: a->kernelctx = kernelctx;
912: delete a->wrapper;
913: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A, IS *is)
918: {
919: Mat_Htool *a;
920: PetscMPIInt rank;
921: const std::vector<PetscInt> *source;
922: const htool::Cluster<PetscReal> *local_source_cluster;
924: PetscFunctionBegin;
925: PetscCall(MatShellGetContext(A, &a));
926: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
927: local_source_cluster = a->source_cluster ? &a->source_cluster->get_cluster_on_partition(rank) : &a->target_cluster->get_cluster_on_partition(rank);
928: source = &local_source_cluster->get_permutation();
929: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), local_source_cluster->get_size(), source->data() + local_source_cluster->get_offset(), PETSC_COPY_VALUES, is));
930: PetscCall(ISSetPermutation(*is));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A, IS *is)
935: {
936: Mat_Htool *a;
937: const std::vector<PetscInt> *target;
938: PetscMPIInt rank;
940: PetscFunctionBegin;
941: PetscCall(MatShellGetContext(A, &a));
942: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
943: target = &a->target_cluster->get_permutation();
944: 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));
945: PetscCall(ISSetPermutation(*is));
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A, PetscBool use)
950: {
951: Mat_Htool *a;
953: PetscFunctionBegin;
954: PetscCall(MatShellGetContext(A, &a));
955: a->permutation = use;
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: static PetscErrorCode MatHtoolUseRecompression_Htool(Mat A, PetscBool use)
960: {
961: Mat_Htool *a;
963: PetscFunctionBegin;
964: PetscCall(MatShellGetContext(A, &a));
965: if (a->recompression != use) A->assembled = PETSC_FALSE;
966: a->recompression = use;
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: #define PETSC_HTOOL_PARAMETER(Type, Name, member) \
971: static PetscErrorCode MatHtoolGet##Name##_Htool(Mat A, Type *v) \
972: { \
973: Mat_Htool *a; \
974: PetscFunctionBegin; \
975: PetscCall(MatShellGetContext(A, &a)); \
976: *v = a->member; \
977: PetscFunctionReturn(PETSC_SUCCESS); \
978: } \
979: static PetscErrorCode MatHtoolSet##Name##_Htool(Mat A, Type v) \
980: { \
981: Mat_Htool *a; \
982: PetscFunctionBegin; \
983: PetscCall(MatShellGetContext(A, &a)); \
984: if (a->member != v) A->assembled = PETSC_FALSE; \
985: a->member = v; \
986: PetscFunctionReturn(PETSC_SUCCESS); \
987: }
989: PETSC_HTOOL_PARAMETER(PetscReal, Epsilon, epsilon)
990: PETSC_HTOOL_PARAMETER(PetscReal, Eta, eta)
991: PETSC_HTOOL_PARAMETER(PetscInt, MaxClusterLeafSize, max_cluster_leaf_size)
992: PETSC_HTOOL_PARAMETER(PetscInt, MinTargetDepth, depth[0])
993: PETSC_HTOOL_PARAMETER(PetscInt, MinSourceDepth, depth[1])
994: PETSC_HTOOL_PARAMETER(PetscBool, BlockTreeConsistency, block_tree_consistency)
995: PETSC_HTOOL_PARAMETER(MatHtoolCompressorType, CompressorType, compressor)
996: PETSC_HTOOL_PARAMETER(MatHtoolClusteringType, ClusteringType, clustering)
998: static PetscErrorCode MatConvert_Htool_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
999: {
1000: Mat C;
1001: Mat_Htool *a;
1002: PetscScalar *array, shift, scale;
1003: PetscInt lda;
1005: PetscFunctionBegin;
1006: 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));
1007: PetscCall(MatShellGetContext(A, &a));
1008: if (reuse == MAT_REUSE_MATRIX) {
1009: C = *B;
1010: PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible dimensions");
1011: PetscCall(MatDenseGetLDA(C, &lda));
1012: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
1013: } else {
1014: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1015: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1016: PetscCall(MatSetType(C, MATDENSE));
1017: PetscCall(MatSetUp(C));
1018: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1019: }
1020: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1021: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1022: PetscCall(MatDenseGetArrayWrite(C, &array));
1023: PetscCallExternalVoid("copy_to_dense_in_user_numbering", htool::copy_to_dense_in_user_numbering(*a->local_hmatrix_view, array));
1024: PetscCall(MatDenseRestoreArrayWrite(C, &array));
1025: PetscCall(MatScale(C, scale));
1026: PetscCall(MatShift(C, shift));
1027: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &C));
1028: else *B = C;
1029: PetscFunctionReturn(PETSC_SUCCESS);
1030: }
1032: static PetscErrorCode GenEntriesTranspose(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *rows, const PetscInt *cols, PetscScalar *ptr, PetscCtx ctx)
1033: {
1034: MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose *)ctx;
1035: PetscScalar *tmp;
1037: PetscFunctionBegin;
1038: PetscCall(generator->kernel(sdim, N, M, cols, rows, ptr, generator->kernelctx));
1039: PetscCall(PetscMalloc1(M * N, &tmp));
1040: PetscCall(PetscArraycpy(tmp, ptr, M * N));
1041: for (PetscInt i = 0; i < M; ++i) {
1042: for (PetscInt j = 0; j < N; ++j) ptr[i + j * M] = tmp[j + i * N];
1043: }
1044: PetscCall(PetscFree(tmp));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: /* naive implementation which keeps a reference to the original Mat */
1049: static PetscErrorCode MatTranspose_Htool(Mat A, MatReuse reuse, Mat *B)
1050: {
1051: Mat C;
1052: Mat_Htool *a, *c;
1053: PetscScalar shift, scale;
1054: PetscInt M = A->rmap->N, N = A->cmap->N, m = A->rmap->n, n = A->cmap->n;
1055: PetscContainer container;
1056: MatHtoolKernelTranspose *kernelt;
1058: PetscFunctionBegin;
1059: 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));
1060: PetscCall(MatShellGetContext(A, &a));
1061: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1062: PetscCheck(reuse != MAT_INPLACE_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatTranspose() with MAT_INPLACE_MATRIX not supported");
1063: if (reuse == MAT_INITIAL_MATRIX) {
1064: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1065: PetscCall(MatSetSizes(C, n, m, N, M));
1066: PetscCall(MatSetType(C, MATHTOOL));
1067: PetscCall(MatSetUp(C));
1068: PetscCall(PetscNew(&kernelt));
1069: PetscCall(PetscObjectContainerCompose((PetscObject)C, "KernelTranspose", kernelt, PetscCtxDestroyDefault));
1070: } else {
1071: C = *B;
1072: PetscCall(PetscObjectQuery((PetscObject)C, "KernelTranspose", (PetscObject *)&container));
1073: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatTranspose() with MAT_INITIAL_MATRIX first");
1074: PetscCall(PetscContainerGetPointer(container, &kernelt));
1075: }
1076: PetscCall(MatShellGetContext(C, &c));
1077: c->dim = a->dim;
1078: c->kernel = GenEntriesTranspose;
1079: if (kernelt->A != A) {
1080: PetscCall(MatDestroy(&kernelt->A));
1081: kernelt->A = A;
1082: PetscCall(PetscObjectReference((PetscObject)A));
1083: }
1084: kernelt->kernel = a->kernel;
1085: kernelt->kernelctx = a->kernelctx;
1086: c->kernelctx = kernelt;
1087: c->max_cluster_leaf_size = a->max_cluster_leaf_size;
1088: c->epsilon = a->epsilon;
1089: c->eta = a->eta;
1090: c->block_tree_consistency = a->block_tree_consistency;
1091: c->permutation = a->permutation;
1092: c->recompression = a->recompression;
1093: c->compressor = a->compressor;
1094: c->clustering = a->clustering;
1095: if (reuse == MAT_INITIAL_MATRIX) {
1096: PetscCall(PetscMalloc1(N * c->dim, &c->gcoords_target));
1097: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_source, N * c->dim));
1098: if (a->gcoords_target != a->gcoords_source) {
1099: PetscCall(PetscMalloc1(M * c->dim, &c->gcoords_source));
1100: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_target, M * c->dim));
1101: } else c->gcoords_source = c->gcoords_target;
1102: }
1103: if (reuse != MAT_INITIAL_MATRIX) C->assembled = PETSC_FALSE; // so that C->was_assembled is not set to PETSC_TRUE
1104: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1105: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1106: PetscCall(MatScale(C, scale));
1107: PetscCall(MatShift(C, shift));
1108: if (reuse == MAT_INITIAL_MATRIX) *B = C;
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: struct MatFactorCtx {
1113: htool::HMatrix<PetscScalar> *hmatrix; /* factorized HMatrix filled by MatFactorNumeric_Htool() */
1114: PetscScalar scale; /* scaling factor from MatShellGetScalingShifts(), applied as inverse scaling in Mat[Mat]Solve() */
1115: char ldlt; /* symmetry ('S' or 'H') of an LDL^T factorization, or 0 for a Cholesky factorization */
1116: };
1118: static PetscErrorCode MatFactorCtxDestroy(PetscCtxRt ctx)
1119: {
1120: MatFactorCtx *data = *reinterpret_cast<MatFactorCtx **>(ctx);
1122: PetscFunctionBegin;
1123: delete data->hmatrix;
1124: PetscCall(PetscFree(data));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: static PetscErrorCode MatDestroy_Factor(Mat F)
1129: {
1130: PetscFunctionBegin;
1131: PetscCall(PetscObjectCompose((PetscObject)F, "HMatrix", nullptr));
1132: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatFactorGetSolverType_C", nullptr));
1133: PetscFunctionReturn(PETSC_SUCCESS);
1134: }
1136: static PetscErrorCode MatFactorGetSolverType_Htool(Mat, MatSolverType *type)
1137: {
1138: PetscFunctionBegin;
1139: *type = MATSOLVERHTOOL;
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: template <char trans>
1144: static inline PetscErrorCode MatSolve_Private(Mat A, htool::Matrix<PetscScalar> &X)
1145: {
1146: PetscContainer container;
1147: MatFactorCtx *data;
1149: PetscFunctionBegin;
1150: PetscCall(PetscObjectQuery((PetscObject)A, "HMatrix", (PetscObject *)&container));
1151: 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");
1152: PetscCall(PetscContainerGetPointer(container, &data));
1153: if (A->factortype == MAT_FACTOR_LU) PetscCallExternalVoid("lu_solve", htool::lu_solve(trans, *data->hmatrix, X));
1154: else if (data->ldlt) PetscCallExternalVoid("ldlt_solve", htool::ldlt_solve(data->ldlt, 'U', *data->hmatrix, X));
1155: else PetscCallExternalVoid("cholesky_solve", htool::cholesky_solve('U', *data->hmatrix, X));
1156: PetscCallExternalVoid("scale", htool::scale(1.0 / data->scale, X));
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
1161: static PetscErrorCode MatSolve_Htool(Mat A, Type b, Type x)
1162: {
1163: htool::Matrix<PetscScalar> v;
1164: PetscScalar *array;
1165: PetscInt n;
1167: PetscFunctionBegin;
1168: PetscCall(VecGetLocalSize(b, &n));
1169: PetscCall(VecCopy(b, x));
1170: PetscCall(VecGetArrayWrite(x, &array));
1171: PetscCallCXX(v.assign(n, 1, array, false));
1172: PetscCall(VecRestoreArrayWrite(x, &array));
1173: PetscCall(MatSolve_Private<trans>(A, v));
1174: PetscFunctionReturn(PETSC_SUCCESS);
1175: }
1177: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
1178: static PetscErrorCode MatSolve_Htool(Mat A, Type B, Type X)
1179: {
1180: htool::Matrix<PetscScalar> v;
1181: PetscScalar *array;
1182: PetscInt m, N, lda;
1184: PetscFunctionBegin;
1185: PetscCall(MatGetLocalSize(B, &m, nullptr));
1186: PetscCall(MatGetLocalSize(B, nullptr, &N));
1187: PetscCall(MatDenseGetLDA(X, &lda));
1188: PetscCheck(lda == X->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, X->rmap->n);
1189: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
1190: PetscCall(MatDenseGetArrayWrite(X, &array));
1191: v.assign(m, N, array, false);
1192: PetscCall(MatDenseRestoreArrayWrite(X, &array));
1193: PetscCall(MatSolve_Private<trans>(A, v));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: template <MatFactorType ftype>
1198: static PetscErrorCode MatFactorNumeric_Htool(Mat F, Mat A, const MatFactorInfo *)
1199: {
1200: Mat_Htool *a;
1201: PetscContainer container;
1202: MatFactorCtx *data;
1203: char symmetry;
1205: PetscFunctionBegin;
1206: PetscCall(MatShellGetContext(A, &a));
1207: PetscCall(PetscObjectQuery((PetscObject)F, "HMatrix", (PetscObject *)&container));
1208: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Mat%sFactorSymbolic() must be called before Mat%sFactorNumeric()", ftype == MAT_FACTOR_LU ? "LU" : "Cholesky", ftype == MAT_FACTOR_LU ? "LU" : "Cholesky");
1209: PetscCall(PetscContainerGetPointer(container, &data));
1210: if (ftype == MAT_FACTOR_LU) PetscCheck(a->local_hmatrix_view->get_UPLO() == 'N', PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "LU factorization requires a MATHTOOL with full storage");
1211: delete data->hmatrix;
1212: data->hmatrix = new htool::HMatrix<PetscScalar>(*a->local_hmatrix_view);
1213: PetscCall(MatShellGetScalingShifts(A, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &data->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));
1214: symmetry = data->hmatrix->get_symmetry();
1215: if (symmetry == 'N') symmetry = MatHtoolGetSymmetry_Private(A); // full storage, symmetry flags may have been set after MatAssemblyEnd()
1216: data->ldlt = ftype == MAT_FACTOR_CHOLESKY && (A->spd == PETSC_BOOL3_FALSE || (PetscDefined(USE_COMPLEX) && symmetry == 'S')) ? (symmetry == 'N' ? 'S' : symmetry) : 0;
1217: if (ftype == MAT_FACTOR_LU) PetscCallExternalVoid("sequential_lu_factorization", htool::sequential_lu_factorization(*data->hmatrix));
1218: else if (data->ldlt) PetscCallExternalVoid("sequential_ldlt_factorization", htool::sequential_ldlt_factorization(data->ldlt, 'U', *data->hmatrix));
1219: else PetscCallExternalVoid("sequential_cholesky_factorization", htool::sequential_cholesky_factorization('U', *data->hmatrix));
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: template <MatFactorType ftype>
1224: PetscErrorCode MatFactorSymbolic_Htool(Mat F, Mat)
1225: {
1226: PetscContainer container;
1227: MatFactorCtx *data;
1229: PetscFunctionBegin;
1230: F->preallocated = PETSC_TRUE;
1231: F->assembled = PETSC_TRUE;
1232: F->ops->solve = MatSolve_Htool<'N', Vec>;
1233: F->ops->matsolve = MatSolve_Htool<'N', Mat>;
1234: if (!PetscDefined(USE_COMPLEX) || ftype == MAT_FACTOR_LU) {
1235: F->ops->solvetranspose = MatSolve_Htool<'T', Vec>;
1236: F->ops->matsolvetranspose = MatSolve_Htool<'T', Mat>;
1237: }
1238: F->ops->destroy = MatDestroy_Factor;
1239: if (ftype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_LU>;
1240: else F->ops->choleskyfactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_CHOLESKY>;
1241: PetscCall(PetscObjectQuery((PetscObject)F, "HMatrix", (PetscObject *)&container));
1242: if (!container) {
1243: PetscCall(PetscNew(&data));
1244: PetscCall(PetscObjectContainerCompose((PetscObject)F, "HMatrix", data, MatFactorCtxDestroy));
1245: }
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: static PetscErrorCode MatLUFactorSymbolic_Htool(Mat F, Mat A, IS, IS, const MatFactorInfo *)
1250: {
1251: PetscFunctionBegin;
1252: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_LU>(F, A));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscErrorCode MatCholeskyFactorSymbolic_Htool(Mat F, Mat A, IS, const MatFactorInfo *)
1257: {
1258: PetscFunctionBegin;
1259: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_CHOLESKY>(F, A));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: static PetscErrorCode MatGetFactor_htool_htool(Mat A, MatFactorType ftype, Mat *F)
1264: {
1265: Mat B;
1266: Mat_Htool *a;
1267: PetscMPIInt size;
1269: PetscFunctionBegin;
1270: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1271: PetscCall(MatShellGetContext(A, &a));
1272: PetscCheck(size == 1, PetscObjectComm((PetscObject)A), PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel MatGetFactor()");
1273: PetscCheck(a->block_tree_consistency, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot factor a MatHtool with inconsistent block tree");
1274: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1275: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1276: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &((PetscObject)B)->type_name));
1277: PetscCall(MatSetUp(B));
1279: B->ops->getinfo = MatGetInfo_External;
1280: B->factortype = ftype;
1281: B->trivialsymbolic = PETSC_TRUE;
1283: PetscCheck(ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_CHOLESKY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_FACTOR_LU and MAT_FACTOR_CHOLESKY are supported");
1284: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_Htool;
1285: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_Htool;
1287: PetscCall(PetscFree(B->solvertype));
1288: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &B->solvertype));
1290: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_Htool));
1291: *F = B;
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }
1295: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Htool(void)
1296: {
1297: PetscFunctionBegin;
1298: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_LU, MatGetFactor_htool_htool));
1299: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_CHOLESKY, MatGetFactor_htool_htool));
1300: PetscFunctionReturn(PETSC_SUCCESS);
1301: }
1303: static PetscErrorCode MatHtoolCreateFromKernel_Htool(Mat A, PetscInt spacedim, const PetscReal coords_target[], const PetscReal coords_source[], MatHtoolKernelFn *kernel, void *kernelctx)
1304: {
1305: Mat_Htool *a;
1307: PetscFunctionBegin;
1308: PetscCall(MatShellGetContext(A, &a));
1309: a->dim = spacedim;
1310: a->kernel = kernel;
1311: a->kernelctx = kernelctx;
1312: PetscCall(PetscCalloc1(A->rmap->N * spacedim, &a->gcoords_target));
1313: PetscCall(PetscArraycpy(a->gcoords_target + A->rmap->rstart * spacedim, coords_target, A->rmap->n * spacedim));
1314: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_target, A->rmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A)));
1315: if (coords_target != coords_source) {
1316: PetscCall(PetscCalloc1(A->cmap->N * spacedim, &a->gcoords_source));
1317: PetscCall(PetscArraycpy(a->gcoords_source + A->cmap->rstart * spacedim, coords_source, A->cmap->n * spacedim));
1318: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_source, A->cmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A)));
1319: } else a->gcoords_source = a->gcoords_target;
1320: PetscFunctionReturn(PETSC_SUCCESS);
1321: }
1323: /*MC
1324: MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.
1326: Use `./configure --download-htool` to install PETSc to use Htool.
1328: Options Database Key:
1329: . -mat_type htool - matrix type to `MATHTOOL`
1331: Level: beginner
1333: Note:
1334: Sequential `MAT_FACTOR_LU` and `MAT_FACTOR_CHOLESKY` factorizations are available with `MATSOLVERHTOOL`. A `MAT_FACTOR_CHOLESKY`
1335: factorization uses the upper triangular part of the matrix. It is a true Cholesky factorization, unless the matrix is flagged
1336: as not symmetric positive definite with `MatSetOption()` and `MAT_SPD` set to `PETSC_FALSE`, or is complex symmetric and not
1337: flagged with `MAT_HERMITIAN`. In these cases, an $LDL^T$ factorization with Bunch-Kaufman pivoting within dense diagonal blocks is
1338: computed, with $L^T$ replaced by $L^H$ for a matrix flagged with `MAT_HERMITIAN`.
1340: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
1341: M*/
1342: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
1343: {
1344: Mat_Htool *a;
1346: PetscFunctionBegin;
1347: PetscCall(MatSetType(A, MATSHELL));
1348: PetscCall(PetscNew(&a));
1349: PetscCall(MatShellSetContext(A, a));
1350: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Htool));
1351: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL_BLOCK, (PetscErrorCodeFn *)MatGetDiagonalBlock_Htool));
1352: PetscCall(MatShellSetOperation(A, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Htool));
1353: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Htool));
1354: PetscCall(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1355: if (!PetscDefined(USE_COMPLEX)) PetscCall(MatShellSetOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1356: A->ops->increaseoverlap = MatIncreaseOverlap_Htool;
1357: A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
1358: PetscCall(MatShellSetOperation(A, MATOP_VIEW, (PetscErrorCodeFn *)MatView_Htool));
1359: PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (PetscErrorCodeFn *)MatSetFromOptions_Htool));
1360: PetscCall(MatShellSetOperation(A, MATOP_GET_ROW, (PetscErrorCodeFn *)MatGetRow_Htool));
1361: PetscCall(MatShellSetOperation(A, MATOP_RESTORE_ROW, (PetscErrorCodeFn *)MatRestoreRow_Htool));
1362: PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (PetscErrorCodeFn *)MatAssemblyEnd_Htool));
1363: PetscCall(MatShellSetOperation(A, MATOP_TRANSPOSE, (PetscErrorCodeFn *)MatTranspose_Htool));
1364: PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Htool));
1365: a->dim = 0;
1366: a->gcoords_target = nullptr;
1367: a->gcoords_source = nullptr;
1368: a->max_cluster_leaf_size = 10;
1369: a->epsilon = PetscSqrtReal(PETSC_SMALL);
1370: a->eta = 10.0;
1371: a->depth[0] = 0;
1372: a->depth[1] = 0;
1373: a->block_tree_consistency = PETSC_TRUE;
1374: a->permutation = PETSC_TRUE;
1375: a->recompression = PETSC_FALSE;
1376: a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
1377: A->assembled = PETSC_FALSE; // MatCreate_Shell() forces this value to PETSC_TRUE
1378: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", MatProductSetFromOptions_Htool));
1379: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", MatProductSetFromOptions_Htool));
1380: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", MatConvert_Htool_Dense));
1381: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", MatConvert_Htool_Dense));
1382: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", MatHtoolGetHierarchicalMat_Htool));
1383: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", MatHtoolSetKernel_Htool));
1384: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", MatHtoolGetPermutationSource_Htool));
1385: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", MatHtoolGetPermutationTarget_Htool));
1386: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", MatHtoolUsePermutation_Htool));
1387: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUseRecompression_C", MatHtoolUseRecompression_Htool));
1388: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEpsilon_C", MatHtoolGetEpsilon_Htool));
1389: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEpsilon_C", MatHtoolSetEpsilon_Htool));
1390: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEta_C", MatHtoolGetEta_Htool));
1391: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEta_C", MatHtoolSetEta_Htool));
1392: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMaxClusterLeafSize_C", MatHtoolGetMaxClusterLeafSize_Htool));
1393: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMaxClusterLeafSize_C", MatHtoolSetMaxClusterLeafSize_Htool));
1394: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinTargetDepth_C", MatHtoolGetMinTargetDepth_Htool));
1395: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinTargetDepth_C", MatHtoolSetMinTargetDepth_Htool));
1396: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinSourceDepth_C", MatHtoolGetMinSourceDepth_Htool));
1397: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinSourceDepth_C", MatHtoolSetMinSourceDepth_Htool));
1398: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetBlockTreeConsistency_C", MatHtoolGetBlockTreeConsistency_Htool));
1399: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetBlockTreeConsistency_C", MatHtoolSetBlockTreeConsistency_Htool));
1400: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetCompressorType_C", MatHtoolGetCompressorType_Htool));
1401: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetCompressorType_C", MatHtoolSetCompressorType_Htool));
1402: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetClusteringType_C", MatHtoolGetClusteringType_Htool));
1403: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetClusteringType_C", MatHtoolSetClusteringType_Htool));
1404: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolCreateFromKernel_C", MatHtoolCreateFromKernel_Htool));
1405: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Immutable));
1406: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
1407: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
1408: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATHTOOL));
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }