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: static PetscErrorCode MatAssemblyEnd_Htool(Mat A, MatAssemblyType)
710: {
711: Mat_Htool *a;
712: const PetscInt *ranges;
713: PetscInt *offset;
714: PetscMPIInt size, rank;
715: char S = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 'H' : (A->symmetric == PETSC_BOOL3_TRUE ? 'S' : 'N'), uplo = S == 'N' ? 'N' : 'U';
716: htool::VirtualGenerator<PetscScalar> *generator = nullptr;
717: htool::ClusterTreeBuilder<PetscReal> recursive_build_strategy;
718: htool::Cluster<PetscReal> *source_cluster;
719: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> compressor;
721: PetscFunctionBegin;
722: for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(HtoolCite); ++i) PetscCall(PetscCitationsRegister(HtoolCitations[i], HtoolCite + i));
723: PetscCall(MatShellGetContext(A, &a));
724: if (A->was_assembled != PETSC_TRUE) {
725: delete a->wrapper;
726: a->target_cluster.reset();
727: a->source_cluster.reset();
728: a->local_hmatrix.reset();
729: a->local_to_local_operator.reset();
730: a->global_to_local_operator.reset();
731: a->distributed_operator_holder.reset();
732: // clustering
733: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
734: PetscCall(PetscMalloc1(2 * size, &offset));
735: PetscCall(MatGetOwnershipRanges(A, &ranges));
736: for (PetscInt i = 0; i < size; ++i) {
737: offset[2 * i] = ranges[i];
738: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
739: }
740: switch (a->clustering) {
741: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
742: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
743: break;
744: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
745: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
746: break;
747: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
748: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
749: break;
750: default:
751: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
752: }
753: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
754: 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));
755: if (a->gcoords_target != a->gcoords_source) {
756: PetscCall(MatGetOwnershipRangesColumn(A, &ranges));
757: for (PetscInt i = 0; i < size; ++i) {
758: offset[2 * i] = ranges[i];
759: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
760: }
761: switch (a->clustering) {
762: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
763: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
764: break;
765: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
766: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::GeometricSplitting<PetscReal>>>());
767: break;
768: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
769: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeBoundingBox<PetscReal>, htool::RegularSplitting<PetscReal>>>());
770: break;
771: default:
772: recursive_build_strategy.set_partitioning_strategy(std::make_shared<htool::Partitioning<PetscReal, htool::ComputeLargestExtent<PetscReal>, htool::RegularSplitting<PetscReal>>>());
773: }
774: recursive_build_strategy.set_maximal_leaf_size(a->max_cluster_leaf_size);
775: 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));
776: S = uplo = 'N';
777: source_cluster = a->source_cluster.get();
778: } else source_cluster = a->target_cluster.get();
779: PetscCall(PetscFree(offset));
780: // generator
781: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
782: else {
783: a->wrapper = nullptr;
784: generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx);
785: }
786: // compressor
787: switch (a->compressor) {
788: case MAT_HTOOL_COMPRESSOR_FULL_ACA:
789: compressor = std::make_shared<htool::fullACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
790: break;
791: case MAT_HTOOL_COMPRESSOR_SVD:
792: compressor = std::make_shared<htool::SVD<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
793: break;
794: default:
795: compressor = std::make_shared<htool::sympartialACA<PetscScalar>>(a->wrapper ? *a->wrapper : *generator, a->target_cluster->get_permutation().data(), source_cluster->get_permutation().data());
796: }
797: // local hierarchical matrix
798: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
799: auto hmatrix_builder = htool::HMatrixTreeBuilder<PetscScalar>(a->epsilon, a->eta, S, uplo);
800: if (a->recompression) {
801: std::shared_ptr<htool::VirtualInternalLowRankGenerator<PetscScalar>> RecompressedLowRankGenerator = std::make_shared<htool::RecompressedLowRankGenerator<PetscScalar>>(*compressor, std::function<void(htool::LowRankMatrix<PetscScalar> &)>(htool::SVD_recompression<PetscScalar>));
802: hmatrix_builder.set_low_rank_generator(RecompressedLowRankGenerator);
803: } else hmatrix_builder.set_low_rank_generator(compressor);
804: hmatrix_builder.set_minimal_target_depth(a->depth[0]);
805: hmatrix_builder.set_minimal_source_depth(a->depth[1]);
806: 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");
807: hmatrix_builder.set_block_tree_consistency(a->block_tree_consistency);
808: a->local_hmatrix = std::make_unique<htool::HMatrix<PetscScalar>>(hmatrix_builder.build(a->wrapper ? *a->wrapper : *generator, *a->target_cluster, *source_cluster, rank, rank));
809: 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);
810: a->distributed_operator_holder = std::make_unique<htool::CustomApproximationBuilder<PetscScalar>>(*a->target_cluster, *source_cluster, PetscObjectComm((PetscObject)A), *a->global_to_local_operator);
811: a->distributed_operator = &a->distributed_operator_holder->distributed_operator;
812: 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));
813: a->local_hmatrix_view = a->local_hmatrix.get();
814: }
815: A->was_assembled = PETSC_FALSE;
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: static PetscErrorCode MatProductNumeric_Htool(Mat C)
820: {
821: Mat_Product *product = C->product;
822: Mat_Htool *a;
823: const PetscScalar *in;
824: PetscScalar *out;
825: PetscInt N, lda;
827: PetscFunctionBegin;
828: MatCheckProduct(C, 1);
829: PetscCall(MatGetSize(C, nullptr, &N));
830: PetscCall(MatDenseGetLDA(C, &lda));
831: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
832: PetscCall(MatDenseGetArrayRead(product->B, &in));
833: PetscCall(MatDenseGetArrayWrite(C, &out));
834: PetscCall(MatShellGetContext(product->A, &a));
835: switch (product->type) {
836: case MATPRODUCT_AB:
837: 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));
838: 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));
839: break;
840: case MATPRODUCT_AtB:
841: 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));
842: 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));
843: break;
844: default:
845: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType %s is not supported", MatProductTypes[product->type]);
846: }
847: PetscCall(MatDenseRestoreArrayWrite(C, &out));
848: PetscCall(MatDenseRestoreArrayRead(product->B, &in));
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
853: {
854: Mat_Product *product = C->product;
855: Mat A, B;
856: PetscBool flg;
858: PetscFunctionBegin;
859: MatCheckProduct(C, 1);
860: A = product->A;
861: B = product->B;
862: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, ""));
863: 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);
864: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
865: if (product->type == MATPRODUCT_AB) PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
866: else PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
867: }
868: PetscCall(MatSetType(C, MATDENSE));
869: PetscCall(MatSetUp(C));
870: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
871: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
872: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
873: C->ops->productsymbolic = nullptr;
874: C->ops->productnumeric = MatProductNumeric_Htool;
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
879: {
880: PetscFunctionBegin;
881: MatCheckProduct(C, 1);
882: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
886: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A, void *distributed_operator)
887: {
888: Mat_Htool *a;
890: PetscFunctionBegin;
891: PetscCall(MatShellGetContext(A, &a));
892: *(const void **)distributed_operator = static_cast<const void *>(a->distributed_operator);
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
897: {
898: Mat_Htool *a;
900: PetscFunctionBegin;
901: PetscCall(MatShellGetContext(A, &a));
902: if (kernel != a->kernel || a->kernelctx != kernelctx) A->assembled = PETSC_FALSE;
903: a->kernel = kernel;
904: a->kernelctx = kernelctx;
905: delete a->wrapper;
906: if (a->kernel) a->wrapper = new WrapperHtool(a->dim, a->kernel, a->kernelctx);
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A, IS *is)
911: {
912: Mat_Htool *a;
913: PetscMPIInt rank;
914: const std::vector<PetscInt> *source;
915: const htool::Cluster<PetscReal> *local_source_cluster;
917: PetscFunctionBegin;
918: PetscCall(MatShellGetContext(A, &a));
919: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
920: local_source_cluster = a->source_cluster ? &a->source_cluster->get_cluster_on_partition(rank) : &a->target_cluster->get_cluster_on_partition(rank);
921: source = &local_source_cluster->get_permutation();
922: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), local_source_cluster->get_size(), source->data() + local_source_cluster->get_offset(), PETSC_COPY_VALUES, is));
923: PetscCall(ISSetPermutation(*is));
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A, IS *is)
928: {
929: Mat_Htool *a;
930: const std::vector<PetscInt> *target;
931: PetscMPIInt rank;
933: PetscFunctionBegin;
934: PetscCall(MatShellGetContext(A, &a));
935: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
936: target = &a->target_cluster->get_permutation();
937: 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));
938: PetscCall(ISSetPermutation(*is));
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A, PetscBool use)
943: {
944: Mat_Htool *a;
946: PetscFunctionBegin;
947: PetscCall(MatShellGetContext(A, &a));
948: a->permutation = use;
949: PetscFunctionReturn(PETSC_SUCCESS);
950: }
952: static PetscErrorCode MatHtoolUseRecompression_Htool(Mat A, PetscBool use)
953: {
954: Mat_Htool *a;
956: PetscFunctionBegin;
957: PetscCall(MatShellGetContext(A, &a));
958: if (a->recompression != use) A->assembled = PETSC_FALSE;
959: a->recompression = use;
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: #define PETSC_HTOOL_PARAMETER(Type, Name, member) \
964: static PetscErrorCode MatHtoolGet##Name##_Htool(Mat A, Type *v) \
965: { \
966: Mat_Htool *a; \
967: PetscFunctionBegin; \
968: PetscCall(MatShellGetContext(A, &a)); \
969: *v = a->member; \
970: PetscFunctionReturn(PETSC_SUCCESS); \
971: } \
972: static PetscErrorCode MatHtoolSet##Name##_Htool(Mat A, Type v) \
973: { \
974: Mat_Htool *a; \
975: PetscFunctionBegin; \
976: PetscCall(MatShellGetContext(A, &a)); \
977: if (a->member != v) A->assembled = PETSC_FALSE; \
978: a->member = v; \
979: PetscFunctionReturn(PETSC_SUCCESS); \
980: }
982: PETSC_HTOOL_PARAMETER(PetscReal, Epsilon, epsilon)
983: PETSC_HTOOL_PARAMETER(PetscReal, Eta, eta)
984: PETSC_HTOOL_PARAMETER(PetscInt, MaxClusterLeafSize, max_cluster_leaf_size)
985: PETSC_HTOOL_PARAMETER(PetscInt, MinTargetDepth, depth[0])
986: PETSC_HTOOL_PARAMETER(PetscInt, MinSourceDepth, depth[1])
987: PETSC_HTOOL_PARAMETER(PetscBool, BlockTreeConsistency, block_tree_consistency)
988: PETSC_HTOOL_PARAMETER(MatHtoolCompressorType, CompressorType, compressor)
989: PETSC_HTOOL_PARAMETER(MatHtoolClusteringType, ClusteringType, clustering)
991: static PetscErrorCode MatConvert_Htool_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
992: {
993: Mat C;
994: Mat_Htool *a;
995: PetscScalar *array, shift, scale;
996: PetscInt lda;
998: PetscFunctionBegin;
999: 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));
1000: PetscCall(MatShellGetContext(A, &a));
1001: if (reuse == MAT_REUSE_MATRIX) {
1002: C = *B;
1003: PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible dimensions");
1004: PetscCall(MatDenseGetLDA(C, &lda));
1005: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
1006: } else {
1007: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1008: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1009: PetscCall(MatSetType(C, MATDENSE));
1010: PetscCall(MatSetUp(C));
1011: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1012: }
1013: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1014: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1015: PetscCall(MatDenseGetArrayWrite(C, &array));
1016: PetscCallExternalVoid("copy_to_dense_in_user_numbering", htool::copy_to_dense_in_user_numbering(*a->local_hmatrix_view, array));
1017: PetscCall(MatDenseRestoreArrayWrite(C, &array));
1018: PetscCall(MatScale(C, scale));
1019: PetscCall(MatShift(C, shift));
1020: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &C));
1021: else *B = C;
1022: PetscFunctionReturn(PETSC_SUCCESS);
1023: }
1025: static PetscErrorCode GenEntriesTranspose(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *rows, const PetscInt *cols, PetscScalar *ptr, PetscCtx ctx)
1026: {
1027: MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose *)ctx;
1028: PetscScalar *tmp;
1030: PetscFunctionBegin;
1031: PetscCall(generator->kernel(sdim, N, M, cols, rows, ptr, generator->kernelctx));
1032: PetscCall(PetscMalloc1(M * N, &tmp));
1033: PetscCall(PetscArraycpy(tmp, ptr, M * N));
1034: for (PetscInt i = 0; i < M; ++i) {
1035: for (PetscInt j = 0; j < N; ++j) ptr[i + j * M] = tmp[j + i * N];
1036: }
1037: PetscCall(PetscFree(tmp));
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /* naive implementation which keeps a reference to the original Mat */
1042: static PetscErrorCode MatTranspose_Htool(Mat A, MatReuse reuse, Mat *B)
1043: {
1044: Mat C;
1045: Mat_Htool *a, *c;
1046: PetscScalar shift, scale;
1047: PetscInt M = A->rmap->N, N = A->cmap->N, m = A->rmap->n, n = A->cmap->n;
1048: PetscContainer container;
1049: MatHtoolKernelTranspose *kernelt;
1051: PetscFunctionBegin;
1052: 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));
1053: PetscCall(MatShellGetContext(A, &a));
1054: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1055: PetscCheck(reuse != MAT_INPLACE_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatTranspose() with MAT_INPLACE_MATRIX not supported");
1056: if (reuse == MAT_INITIAL_MATRIX) {
1057: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1058: PetscCall(MatSetSizes(C, n, m, N, M));
1059: PetscCall(MatSetType(C, MATHTOOL));
1060: PetscCall(MatSetUp(C));
1061: PetscCall(PetscNew(&kernelt));
1062: PetscCall(PetscObjectContainerCompose((PetscObject)C, "KernelTranspose", kernelt, PetscCtxDestroyDefault));
1063: } else {
1064: C = *B;
1065: PetscCall(PetscObjectQuery((PetscObject)C, "KernelTranspose", (PetscObject *)&container));
1066: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatTranspose() with MAT_INITIAL_MATRIX first");
1067: PetscCall(PetscContainerGetPointer(container, &kernelt));
1068: }
1069: PetscCall(MatShellGetContext(C, &c));
1070: c->dim = a->dim;
1071: c->kernel = GenEntriesTranspose;
1072: if (kernelt->A != A) {
1073: PetscCall(MatDestroy(&kernelt->A));
1074: kernelt->A = A;
1075: PetscCall(PetscObjectReference((PetscObject)A));
1076: }
1077: kernelt->kernel = a->kernel;
1078: kernelt->kernelctx = a->kernelctx;
1079: c->kernelctx = kernelt;
1080: c->max_cluster_leaf_size = a->max_cluster_leaf_size;
1081: c->epsilon = a->epsilon;
1082: c->eta = a->eta;
1083: c->block_tree_consistency = a->block_tree_consistency;
1084: c->permutation = a->permutation;
1085: c->recompression = a->recompression;
1086: c->compressor = a->compressor;
1087: c->clustering = a->clustering;
1088: if (reuse == MAT_INITIAL_MATRIX) {
1089: PetscCall(PetscMalloc1(N * c->dim, &c->gcoords_target));
1090: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_source, N * c->dim));
1091: if (a->gcoords_target != a->gcoords_source) {
1092: PetscCall(PetscMalloc1(M * c->dim, &c->gcoords_source));
1093: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_target, M * c->dim));
1094: } else c->gcoords_source = c->gcoords_target;
1095: }
1096: if (reuse != MAT_INITIAL_MATRIX) C->assembled = PETSC_FALSE; // so that C->was_assembled is not set to PETSC_TRUE
1097: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1098: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1099: PetscCall(MatScale(C, scale));
1100: PetscCall(MatShift(C, shift));
1101: if (reuse == MAT_INITIAL_MATRIX) *B = C;
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: struct MatFactorCtx {
1106: htool::HMatrix<PetscScalar> *hmatrix; /* factorized HMatrix filled by MatFactorNumeric_Htool() */
1107: PetscScalar scale; /* scaling factor from MatShellGetScalingShifts(), applied as inverse scaling in Mat[Mat]Solve() */
1108: };
1110: static PetscErrorCode MatFactorCtxDestroy(PetscCtxRt ctx)
1111: {
1112: MatFactorCtx *data = *reinterpret_cast<MatFactorCtx **>(ctx);
1114: PetscFunctionBegin;
1115: delete data->hmatrix;
1116: PetscCall(PetscFree(data));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: static PetscErrorCode MatDestroy_Factor(Mat F)
1121: {
1122: PetscFunctionBegin;
1123: PetscCall(PetscObjectCompose((PetscObject)F, "HMatrix", nullptr));
1124: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatFactorGetSolverType_C", nullptr));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: static PetscErrorCode MatFactorGetSolverType_Htool(Mat, MatSolverType *type)
1129: {
1130: PetscFunctionBegin;
1131: *type = MATSOLVERHTOOL;
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: template <char trans>
1136: static inline PetscErrorCode MatSolve_Private(Mat A, htool::Matrix<PetscScalar> &X)
1137: {
1138: PetscContainer container;
1139: MatFactorCtx *data;
1141: PetscFunctionBegin;
1142: PetscCall(PetscObjectQuery((PetscObject)A, "HMatrix", (PetscObject *)&container));
1143: 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");
1144: PetscCall(PetscContainerGetPointer(container, &data));
1145: if (A->factortype == MAT_FACTOR_LU) PetscCallExternalVoid("lu_solve", htool::lu_solve(trans, *data->hmatrix, X));
1146: else PetscCallExternalVoid("cholesky_solve", htool::cholesky_solve('U', *data->hmatrix, X));
1147: PetscCallExternalVoid("scale", htool::scale(1.0 / data->scale, X));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
1152: static PetscErrorCode MatSolve_Htool(Mat A, Type b, Type x)
1153: {
1154: htool::Matrix<PetscScalar> v;
1155: PetscScalar *array;
1156: PetscInt n;
1158: PetscFunctionBegin;
1159: PetscCall(VecGetLocalSize(b, &n));
1160: PetscCall(VecCopy(b, x));
1161: PetscCall(VecGetArrayWrite(x, &array));
1162: PetscCallCXX(v.assign(n, 1, array, false));
1163: PetscCall(VecRestoreArrayWrite(x, &array));
1164: PetscCall(MatSolve_Private<trans>(A, v));
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: template <char trans, class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
1169: static PetscErrorCode MatSolve_Htool(Mat A, Type B, Type X)
1170: {
1171: htool::Matrix<PetscScalar> v;
1172: PetscScalar *array;
1173: PetscInt m, N, lda;
1175: PetscFunctionBegin;
1176: PetscCall(MatGetLocalSize(B, &m, nullptr));
1177: PetscCall(MatGetLocalSize(B, nullptr, &N));
1178: PetscCall(MatDenseGetLDA(X, &lda));
1179: PetscCheck(lda == X->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, X->rmap->n);
1180: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
1181: PetscCall(MatDenseGetArrayWrite(X, &array));
1182: v.assign(m, N, array, false);
1183: PetscCall(MatDenseRestoreArrayWrite(X, &array));
1184: PetscCall(MatSolve_Private<trans>(A, v));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: template <MatFactorType ftype>
1189: static PetscErrorCode MatFactorNumeric_Htool(Mat F, Mat A, const MatFactorInfo *)
1190: {
1191: Mat_Htool *a;
1192: PetscContainer container;
1193: MatFactorCtx *data;
1195: PetscFunctionBegin;
1196: PetscCall(MatShellGetContext(A, &a));
1197: PetscCall(PetscObjectQuery((PetscObject)F, "HMatrix", (PetscObject *)&container));
1198: 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");
1199: PetscCall(PetscContainerGetPointer(container, &data));
1200: 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");
1201: delete data->hmatrix;
1202: data->hmatrix = new htool::HMatrix<PetscScalar>(*a->local_hmatrix_view);
1203: 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));
1204: if (ftype == MAT_FACTOR_LU) PetscCallExternalVoid("sequential_lu_factorization", htool::sequential_lu_factorization(*data->hmatrix));
1205: else PetscCallExternalVoid("sequential_cholesky_factorization", htool::sequential_cholesky_factorization('U', *data->hmatrix));
1206: PetscFunctionReturn(PETSC_SUCCESS);
1207: }
1209: template <MatFactorType ftype>
1210: PetscErrorCode MatFactorSymbolic_Htool(Mat F, Mat)
1211: {
1212: PetscContainer container;
1213: MatFactorCtx *data;
1215: PetscFunctionBegin;
1216: F->preallocated = PETSC_TRUE;
1217: F->assembled = PETSC_TRUE;
1218: F->ops->solve = MatSolve_Htool<'N', Vec>;
1219: F->ops->matsolve = MatSolve_Htool<'N', Mat>;
1220: if (!PetscDefined(USE_COMPLEX) || ftype == MAT_FACTOR_LU) {
1221: F->ops->solvetranspose = MatSolve_Htool<'T', Vec>;
1222: F->ops->matsolvetranspose = MatSolve_Htool<'T', Mat>;
1223: }
1224: F->ops->destroy = MatDestroy_Factor;
1225: if (ftype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_LU>;
1226: else F->ops->choleskyfactornumeric = MatFactorNumeric_Htool<MAT_FACTOR_CHOLESKY>;
1227: PetscCall(PetscObjectQuery((PetscObject)F, "HMatrix", (PetscObject *)&container));
1228: if (!container) {
1229: PetscCall(PetscNew(&data));
1230: PetscCall(PetscObjectContainerCompose((PetscObject)F, "HMatrix", data, MatFactorCtxDestroy));
1231: }
1232: PetscFunctionReturn(PETSC_SUCCESS);
1233: }
1235: static PetscErrorCode MatLUFactorSymbolic_Htool(Mat F, Mat A, IS, IS, const MatFactorInfo *)
1236: {
1237: PetscFunctionBegin;
1238: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_LU>(F, A));
1239: PetscFunctionReturn(PETSC_SUCCESS);
1240: }
1242: static PetscErrorCode MatCholeskyFactorSymbolic_Htool(Mat F, Mat A, IS, const MatFactorInfo *)
1243: {
1244: PetscFunctionBegin;
1245: PetscCall(MatFactorSymbolic_Htool<MAT_FACTOR_CHOLESKY>(F, A));
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: static PetscErrorCode MatGetFactor_htool_htool(Mat A, MatFactorType ftype, Mat *F)
1250: {
1251: Mat B;
1252: Mat_Htool *a;
1253: PetscMPIInt size;
1255: PetscFunctionBegin;
1256: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1257: PetscCall(MatShellGetContext(A, &a));
1258: PetscCheck(size == 1, PetscObjectComm((PetscObject)A), PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel MatGetFactor()");
1259: PetscCheck(a->block_tree_consistency, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot factor a MatHtool with inconsistent block tree");
1260: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1261: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1262: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &((PetscObject)B)->type_name));
1263: PetscCall(MatSetUp(B));
1265: B->ops->getinfo = MatGetInfo_External;
1266: B->factortype = ftype;
1267: B->trivialsymbolic = PETSC_TRUE;
1269: 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");
1270: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_Htool;
1271: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_Htool;
1273: PetscCall(PetscFree(B->solvertype));
1274: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &B->solvertype));
1276: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_Htool));
1277: *F = B;
1278: PetscFunctionReturn(PETSC_SUCCESS);
1279: }
1281: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Htool(void)
1282: {
1283: PetscFunctionBegin;
1284: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_LU, MatGetFactor_htool_htool));
1285: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATHTOOL, MAT_FACTOR_CHOLESKY, MatGetFactor_htool_htool));
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: static PetscErrorCode MatHtoolCreateFromKernel_Htool(Mat A, PetscInt spacedim, const PetscReal coords_target[], const PetscReal coords_source[], MatHtoolKernelFn *kernel, void *kernelctx)
1290: {
1291: Mat_Htool *a;
1293: PetscFunctionBegin;
1294: PetscCall(MatShellGetContext(A, &a));
1295: a->dim = spacedim;
1296: a->kernel = kernel;
1297: a->kernelctx = kernelctx;
1298: PetscCall(PetscCalloc1(A->rmap->N * spacedim, &a->gcoords_target));
1299: PetscCall(PetscArraycpy(a->gcoords_target + A->rmap->rstart * spacedim, coords_target, A->rmap->n * spacedim));
1300: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_target, A->rmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A)));
1301: if (coords_target != coords_source) {
1302: PetscCall(PetscCalloc1(A->cmap->N * spacedim, &a->gcoords_source));
1303: PetscCall(PetscArraycpy(a->gcoords_source + A->cmap->rstart * spacedim, coords_source, A->cmap->n * spacedim));
1304: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_source, A->cmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A)));
1305: } else a->gcoords_source = a->gcoords_target;
1306: PetscFunctionReturn(PETSC_SUCCESS);
1307: }
1309: /*MC
1310: MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.
1312: Use `./configure --download-htool` to install PETSc to use Htool.
1314: Options Database Key:
1315: . -mat_type htool - matrix type to `MATHTOOL`
1317: Level: beginner
1319: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
1320: M*/
1321: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
1322: {
1323: Mat_Htool *a;
1325: PetscFunctionBegin;
1326: PetscCall(MatSetType(A, MATSHELL));
1327: PetscCall(PetscNew(&a));
1328: PetscCall(MatShellSetContext(A, a));
1329: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Htool));
1330: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL_BLOCK, (PetscErrorCodeFn *)MatGetDiagonalBlock_Htool));
1331: PetscCall(MatShellSetOperation(A, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Htool));
1332: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Htool));
1333: PetscCall(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1334: if (!PetscDefined(USE_COMPLEX)) PetscCall(MatShellSetOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Htool));
1335: A->ops->increaseoverlap = MatIncreaseOverlap_Htool;
1336: A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
1337: PetscCall(MatShellSetOperation(A, MATOP_VIEW, (PetscErrorCodeFn *)MatView_Htool));
1338: PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (PetscErrorCodeFn *)MatSetFromOptions_Htool));
1339: PetscCall(MatShellSetOperation(A, MATOP_GET_ROW, (PetscErrorCodeFn *)MatGetRow_Htool));
1340: PetscCall(MatShellSetOperation(A, MATOP_RESTORE_ROW, (PetscErrorCodeFn *)MatRestoreRow_Htool));
1341: PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (PetscErrorCodeFn *)MatAssemblyEnd_Htool));
1342: PetscCall(MatShellSetOperation(A, MATOP_TRANSPOSE, (PetscErrorCodeFn *)MatTranspose_Htool));
1343: PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Htool));
1344: a->dim = 0;
1345: a->gcoords_target = nullptr;
1346: a->gcoords_source = nullptr;
1347: a->max_cluster_leaf_size = 10;
1348: a->epsilon = PetscSqrtReal(PETSC_SMALL);
1349: a->eta = 10.0;
1350: a->depth[0] = 0;
1351: a->depth[1] = 0;
1352: a->block_tree_consistency = PETSC_TRUE;
1353: a->permutation = PETSC_TRUE;
1354: a->recompression = PETSC_FALSE;
1355: a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
1356: A->assembled = PETSC_FALSE; // MatCreate_Shell() forces this value to PETSC_TRUE
1357: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", MatProductSetFromOptions_Htool));
1358: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", MatProductSetFromOptions_Htool));
1359: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", MatConvert_Htool_Dense));
1360: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", MatConvert_Htool_Dense));
1361: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", MatHtoolGetHierarchicalMat_Htool));
1362: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", MatHtoolSetKernel_Htool));
1363: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", MatHtoolGetPermutationSource_Htool));
1364: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", MatHtoolGetPermutationTarget_Htool));
1365: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", MatHtoolUsePermutation_Htool));
1366: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUseRecompression_C", MatHtoolUseRecompression_Htool));
1367: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEpsilon_C", MatHtoolGetEpsilon_Htool));
1368: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEpsilon_C", MatHtoolSetEpsilon_Htool));
1369: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetEta_C", MatHtoolGetEta_Htool));
1370: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetEta_C", MatHtoolSetEta_Htool));
1371: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMaxClusterLeafSize_C", MatHtoolGetMaxClusterLeafSize_Htool));
1372: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMaxClusterLeafSize_C", MatHtoolSetMaxClusterLeafSize_Htool));
1373: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinTargetDepth_C", MatHtoolGetMinTargetDepth_Htool));
1374: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinTargetDepth_C", MatHtoolSetMinTargetDepth_Htool));
1375: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetMinSourceDepth_C", MatHtoolGetMinSourceDepth_Htool));
1376: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetMinSourceDepth_C", MatHtoolSetMinSourceDepth_Htool));
1377: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetBlockTreeConsistency_C", MatHtoolGetBlockTreeConsistency_Htool));
1378: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetBlockTreeConsistency_C", MatHtoolSetBlockTreeConsistency_Htool));
1379: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetCompressorType_C", MatHtoolGetCompressorType_Htool));
1380: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetCompressorType_C", MatHtoolSetCompressorType_Htool));
1381: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetClusteringType_C", MatHtoolGetClusteringType_Htool));
1382: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetClusteringType_C", MatHtoolSetClusteringType_Htool));
1383: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolCreateFromKernel_C", MatHtoolCreateFromKernel_Htool));
1384: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Immutable));
1385: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
1386: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
1387: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATHTOOL));
1388: PetscFunctionReturn(PETSC_SUCCESS);
1389: }