Actual source code: htool.cxx
1: #include <../src/mat/impls/htool/htool.hpp>
2: #include <petscblaslapack.h>
3: #include <set>
5: const char *const MatHtoolCompressorTypes[] = {"sympartialACA", "fullACA", "SVD"};
6: const char *const MatHtoolClusteringTypes[] = {"PCARegular", "PCAGeometric", "BoundingBox1Regular", "BoundingBox1Geometric"};
7: const char HtoolCitation[] = "@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: static PetscBool HtoolCite = PETSC_FALSE;
19: static PetscErrorCode MatGetDiagonal_Htool(Mat A, Vec v)
20: {
21: Mat_Htool *a = (Mat_Htool *)A->data;
22: PetscScalar *x;
23: PetscBool flg;
25: PetscFunctionBegin;
26: PetscCall(MatHasCongruentLayouts(A, &flg));
27: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent layouts supported");
28: PetscCall(VecGetArrayWrite(v, &x));
29: a->hmatrix->copy_local_diagonal(x);
30: PetscCall(VecRestoreArrayWrite(v, &x));
31: PetscCall(VecScale(v, a->s));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode MatGetDiagonalBlock_Htool(Mat A, Mat *b)
36: {
37: Mat_Htool *a = (Mat_Htool *)A->data;
38: Mat B;
39: PetscScalar *ptr;
40: PetscBool flg;
42: PetscFunctionBegin;
43: PetscCall(MatHasCongruentLayouts(A, &flg));
44: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent layouts supported");
45: PetscCall(PetscObjectQuery((PetscObject)A, "DiagonalBlock", (PetscObject *)&B)); /* same logic as in MatGetDiagonalBlock_MPIDense() */
46: if (!B) {
47: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, A->rmap->n, A->rmap->n, A->rmap->n, nullptr, &B));
48: PetscCall(MatDenseGetArrayWrite(B, &ptr));
49: a->hmatrix->copy_local_diagonal_block(ptr);
50: PetscCall(MatDenseRestoreArrayWrite(B, &ptr));
51: PetscCall(MatPropagateSymmetryOptions(A, B));
52: PetscCall(MatScale(B, a->s));
53: PetscCall(PetscObjectCompose((PetscObject)A, "DiagonalBlock", (PetscObject)B));
54: *b = B;
55: PetscCall(MatDestroy(&B));
56: } else *b = B;
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode MatMult_Htool(Mat A, Vec x, Vec y)
61: {
62: Mat_Htool *a = (Mat_Htool *)A->data;
63: const PetscScalar *in;
64: PetscScalar *out;
66: PetscFunctionBegin;
67: PetscCall(VecGetArrayRead(x, &in));
68: PetscCall(VecGetArrayWrite(y, &out));
69: a->hmatrix->mvprod_local_to_local(in, out);
70: PetscCall(VecRestoreArrayRead(x, &in));
71: PetscCall(VecRestoreArrayWrite(y, &out));
72: PetscCall(VecScale(y, a->s));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /* naive implementation of MatMultAdd() needed for FEM-BEM coupling via MATNEST */
77: static PetscErrorCode MatMultAdd_Htool(Mat A, Vec v1, Vec v2, Vec v3)
78: {
79: Mat_Htool *a = (Mat_Htool *)A->data;
80: Vec tmp;
81: const PetscScalar scale = a->s;
83: PetscFunctionBegin;
84: PetscCall(VecDuplicate(v2, &tmp));
85: PetscCall(VecCopy(v2, v3)); /* no-op in MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]) since VecCopy() checks for x == y */
86: a->s = 1.0; /* set s to 1.0 since VecAXPY() may be used to scale the MatMult() output Vec */
87: PetscCall(MatMult_Htool(A, v1, tmp));
88: PetscCall(VecAXPY(v3, scale, tmp));
89: PetscCall(VecDestroy(&tmp));
90: a->s = scale; /* set s back to its original value */
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: static PetscErrorCode MatMultTranspose_Htool(Mat A, Vec x, Vec y)
95: {
96: Mat_Htool *a = (Mat_Htool *)A->data;
97: const PetscScalar *in;
98: PetscScalar *out;
100: PetscFunctionBegin;
101: PetscCall(VecGetArrayRead(x, &in));
102: PetscCall(VecGetArrayWrite(y, &out));
103: a->hmatrix->mvprod_transp_local_to_local(in, out);
104: PetscCall(VecRestoreArrayRead(x, &in));
105: PetscCall(VecRestoreArrayWrite(y, &out));
106: PetscCall(VecScale(y, a->s));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode MatIncreaseOverlap_Htool(Mat A, PetscInt is_max, IS is[], PetscInt ov)
111: {
112: std::set<PetscInt> set;
113: const PetscInt *idx;
114: PetscInt *oidx, size, bs[2];
115: PetscMPIInt csize;
117: PetscFunctionBegin;
118: PetscCall(MatGetBlockSizes(A, bs, bs + 1));
119: if (bs[0] != bs[1]) bs[0] = 1;
120: for (PetscInt i = 0; i < is_max; ++i) {
121: /* basic implementation that adds indices by shifting an IS by -ov, -ov+1..., -1, 1..., ov-1, ov */
122: /* needed to avoid subdomain matrices to replicate A since it is dense */
123: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is[i]), &csize));
124: PetscCheck(csize == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported parallel IS");
125: PetscCall(ISGetSize(is[i], &size));
126: PetscCall(ISGetIndices(is[i], &idx));
127: for (PetscInt j = 0; j < size; ++j) {
128: set.insert(idx[j]);
129: for (PetscInt k = 1; k <= ov; ++k) { /* for each layer of overlap */
130: if (idx[j] - k >= 0) set.insert(idx[j] - k); /* do not insert negative indices */
131: 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 */
132: }
133: }
134: PetscCall(ISRestoreIndices(is[i], &idx));
135: PetscCall(ISDestroy(is + i));
136: if (bs[0] > 1) {
137: for (std::set<PetscInt>::iterator it = set.cbegin(); it != set.cend(); it++) {
138: std::vector<PetscInt> block(bs[0]);
139: std::iota(block.begin(), block.end(), (*it / bs[0]) * bs[0]);
140: set.insert(block.cbegin(), block.cend());
141: }
142: }
143: size = set.size(); /* size with overlap */
144: PetscCall(PetscMalloc1(size, &oidx));
145: for (const PetscInt j : set) *oidx++ = j;
146: oidx -= size;
147: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, oidx, PETSC_OWN_POINTER, is + i));
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode MatCreateSubMatrices_Htool(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
153: {
154: Mat_Htool *a = (Mat_Htool *)A->data;
155: Mat D, B, BT;
156: const PetscScalar *copy;
157: PetscScalar *ptr;
158: const PetscInt *idxr, *idxc, *it;
159: PetscInt nrow, m, i;
160: PetscBool flg;
162: PetscFunctionBegin;
163: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
164: for (i = 0; i < n; ++i) {
165: PetscCall(ISGetLocalSize(irow[i], &nrow));
166: PetscCall(ISGetLocalSize(icol[i], &m));
167: PetscCall(ISGetIndices(irow[i], &idxr));
168: PetscCall(ISGetIndices(icol[i], &idxc));
169: if (scall != MAT_REUSE_MATRIX) PetscCall(MatCreateDense(PETSC_COMM_SELF, nrow, m, nrow, m, nullptr, (*submat) + i));
170: PetscCall(MatDenseGetArrayWrite((*submat)[i], &ptr));
171: if (irow[i] == icol[i]) { /* same row and column IS? */
172: PetscCall(MatHasCongruentLayouts(A, &flg));
173: if (flg) {
174: PetscCall(ISSorted(irow[i], &flg));
175: if (flg) { /* sorted IS? */
176: it = std::lower_bound(idxr, idxr + nrow, A->rmap->rstart);
177: if (it != idxr + nrow && *it == A->rmap->rstart) { /* rmap->rstart in IS? */
178: if (std::distance(idxr, it) + A->rmap->n <= nrow) { /* long enough IS to store the local diagonal block? */
179: for (PetscInt j = 0; j < A->rmap->n && flg; ++j)
180: if (PetscUnlikely(it[j] != A->rmap->rstart + j)) flg = PETSC_FALSE;
181: if (flg) { /* complete local diagonal block in IS? */
182: /* fast extraction when the local diagonal block is part of the submatrix, e.g., for PCASM or PCHPDDM
183: * [ B C E ]
184: * A = [ B D E ]
185: * [ B F E ]
186: */
187: m = std::distance(idxr, it); /* shift of the coefficient (0,0) of block D from above */
188: PetscCall(MatGetDiagonalBlock_Htool(A, &D));
189: PetscCall(MatDenseGetArrayRead(D, ©));
190: for (PetscInt k = 0; k < A->rmap->n; ++k) { PetscCall(PetscArraycpy(ptr + (m + k) * nrow + m, copy + k * A->rmap->n, A->rmap->n)); /* block D from above */ }
191: PetscCall(MatDenseRestoreArrayRead(D, ©));
192: if (m) {
193: a->wrapper->copy_submatrix(nrow, m, idxr, idxc, ptr); /* vertical block B from above */
194: /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
195: if (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE) {
196: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, m, A->rmap->n, m, ptr + m, &B));
197: PetscCall(MatDenseSetLDA(B, nrow));
198: PetscCall(MatCreateDense(PETSC_COMM_SELF, m, A->rmap->n, m, A->rmap->n, ptr + m * nrow, &BT));
199: PetscCall(MatDenseSetLDA(BT, nrow));
200: if (A->hermitian == PETSC_BOOL3_TRUE && PetscDefined(USE_COMPLEX)) {
201: PetscCall(MatHermitianTranspose(B, MAT_REUSE_MATRIX, &BT));
202: } else {
203: PetscCall(MatTransposeSetPrecursor(B, BT));
204: PetscCall(MatTranspose(B, MAT_REUSE_MATRIX, &BT));
205: }
206: PetscCall(MatDestroy(&B));
207: PetscCall(MatDestroy(&BT));
208: } else {
209: for (PetscInt k = 0; k < A->rmap->n; ++k) { /* block C from above */
210: a->wrapper->copy_submatrix(m, 1, idxr, idxc + m + k, ptr + (m + k) * nrow);
211: }
212: }
213: }
214: if (m + A->rmap->n != nrow) {
215: a->wrapper->copy_submatrix(nrow, std::distance(it + A->rmap->n, idxr + nrow), idxr, idxc + m + A->rmap->n, ptr + (m + A->rmap->n) * nrow); /* vertical block E from above */
216: /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
217: if (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE) {
218: PetscCall(MatCreateDense(PETSC_COMM_SELF, A->rmap->n, nrow - (m + A->rmap->n), A->rmap->n, nrow - (m + A->rmap->n), ptr + (m + A->rmap->n) * nrow + m, &B));
219: PetscCall(MatDenseSetLDA(B, nrow));
220: PetscCall(MatCreateDense(PETSC_COMM_SELF, nrow - (m + A->rmap->n), A->rmap->n, nrow - (m + A->rmap->n), A->rmap->n, ptr + m * nrow + m + A->rmap->n, &BT));
221: PetscCall(MatDenseSetLDA(BT, nrow));
222: if (A->hermitian == PETSC_BOOL3_TRUE && PetscDefined(USE_COMPLEX)) {
223: PetscCall(MatHermitianTranspose(B, MAT_REUSE_MATRIX, &BT));
224: } else {
225: PetscCall(MatTransposeSetPrecursor(B, BT));
226: PetscCall(MatTranspose(B, MAT_REUSE_MATRIX, &BT));
227: }
228: PetscCall(MatDestroy(&B));
229: PetscCall(MatDestroy(&BT));
230: } else {
231: for (PetscInt k = 0; k < A->rmap->n; ++k) { /* block F from above */
232: a->wrapper->copy_submatrix(std::distance(it + A->rmap->n, idxr + nrow), 1, it + A->rmap->n, idxc + m + k, ptr + (m + k) * nrow + m + A->rmap->n);
233: }
234: }
235: }
236: } /* complete local diagonal block not in IS */
237: } else flg = PETSC_FALSE; /* IS not long enough to store the local diagonal block */
238: } else flg = PETSC_FALSE; /* rmap->rstart not in IS */
239: } /* unsorted IS */
240: }
241: } else flg = PETSC_FALSE; /* different row and column IS */
242: if (!flg) a->wrapper->copy_submatrix(nrow, m, idxr, idxc, ptr); /* reassemble everything */
243: PetscCall(ISRestoreIndices(irow[i], &idxr));
244: PetscCall(ISRestoreIndices(icol[i], &idxc));
245: PetscCall(MatDenseRestoreArrayWrite((*submat)[i], &ptr));
246: PetscCall(MatScale((*submat)[i], a->s));
247: }
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode MatDestroy_Htool(Mat A)
252: {
253: Mat_Htool *a = (Mat_Htool *)A->data;
254: PetscContainer container;
255: MatHtoolKernelTranspose *kernelt;
257: PetscFunctionBegin;
258: PetscCall(PetscObjectChangeTypeName((PetscObject)A, nullptr));
259: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", nullptr));
260: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", nullptr));
261: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", nullptr));
262: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", nullptr));
263: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", nullptr));
264: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", nullptr));
265: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", nullptr));
266: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", nullptr));
267: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", nullptr));
268: PetscCall(PetscObjectQuery((PetscObject)A, "KernelTranspose", (PetscObject *)&container));
269: if (container) { /* created in MatTranspose_Htool() */
270: PetscCall(PetscContainerGetPointer(container, (void **)&kernelt));
271: PetscCall(MatDestroy(&kernelt->A));
272: PetscCall(PetscFree(kernelt));
273: PetscCall(PetscContainerDestroy(&container));
274: PetscCall(PetscObjectCompose((PetscObject)A, "KernelTranspose", nullptr));
275: }
276: if (a->gcoords_source != a->gcoords_target) PetscCall(PetscFree(a->gcoords_source));
277: PetscCall(PetscFree(a->gcoords_target));
278: PetscCall(PetscFree2(a->work_source, a->work_target));
279: delete a->wrapper;
280: delete a->hmatrix;
281: PetscCall(PetscFree(A->data));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode MatView_Htool(Mat A, PetscViewer pv)
286: {
287: Mat_Htool *a = (Mat_Htool *)A->data;
288: PetscBool flg;
290: PetscFunctionBegin;
291: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &flg));
292: if (flg) {
293: PetscCall(PetscViewerASCIIPrintf(pv, "symmetry: %c\n", a->hmatrix->get_symmetry_type()));
294: if (PetscAbsScalar(a->s - 1.0) > PETSC_MACHINE_EPSILON) {
295: #if defined(PETSC_USE_COMPLEX)
296: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g+%gi\n", (double)PetscRealPart(a->s), (double)PetscImaginaryPart(a->s)));
297: #else
298: PetscCall(PetscViewerASCIIPrintf(pv, "scaling: %g\n", (double)a->s));
299: #endif
300: }
301: PetscCall(PetscViewerASCIIPrintf(pv, "minimum cluster size: %" PetscInt_FMT "\n", a->bs[0]));
302: PetscCall(PetscViewerASCIIPrintf(pv, "maximum block size: %" PetscInt_FMT "\n", a->bs[1]));
303: PetscCall(PetscViewerASCIIPrintf(pv, "epsilon: %g\n", (double)a->epsilon));
304: PetscCall(PetscViewerASCIIPrintf(pv, "eta: %g\n", (double)a->eta));
305: PetscCall(PetscViewerASCIIPrintf(pv, "minimum target depth: %" PetscInt_FMT "\n", a->depth[0]));
306: PetscCall(PetscViewerASCIIPrintf(pv, "minimum source depth: %" PetscInt_FMT "\n", a->depth[1]));
307: PetscCall(PetscViewerASCIIPrintf(pv, "compressor: %s\n", MatHtoolCompressorTypes[a->compressor]));
308: PetscCall(PetscViewerASCIIPrintf(pv, "clustering: %s\n", MatHtoolClusteringTypes[a->clustering]));
309: PetscCall(PetscViewerASCIIPrintf(pv, "compression ratio: %s\n", a->hmatrix->get_infos("Compression_ratio").c_str()));
310: PetscCall(PetscViewerASCIIPrintf(pv, "space saving: %s\n", a->hmatrix->get_infos("Space_saving").c_str()));
311: PetscCall(PetscViewerASCIIPrintf(pv, "number of dense (resp. low rank) matrices: %s (resp. %s)\n", a->hmatrix->get_infos("Number_of_dmat").c_str(), a->hmatrix->get_infos("Number_of_lrmat").c_str()));
312: PetscCall(PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) dense block sizes: (%s, %s, %s)\n", a->hmatrix->get_infos("Dense_block_size_min").c_str(), a->hmatrix->get_infos("Dense_block_size_mean").c_str(),
313: a->hmatrix->get_infos("Dense_block_size_max").c_str()));
314: PetscCall(PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) low rank block sizes: (%s, %s, %s)\n", a->hmatrix->get_infos("Low_rank_block_size_min").c_str(), a->hmatrix->get_infos("Low_rank_block_size_mean").c_str(),
315: a->hmatrix->get_infos("Low_rank_block_size_max").c_str()));
316: PetscCall(PetscViewerASCIIPrintf(pv, "(minimum, mean, maximum) ranks: (%s, %s, %s)\n", a->hmatrix->get_infos("Rank_min").c_str(), a->hmatrix->get_infos("Rank_mean").c_str(), a->hmatrix->get_infos("Rank_max").c_str()));
317: }
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode MatScale_Htool(Mat A, PetscScalar s)
322: {
323: Mat_Htool *a = (Mat_Htool *)A->data;
325: PetscFunctionBegin;
326: a->s *= s;
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */
331: static PetscErrorCode MatGetRow_Htool(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
332: {
333: Mat_Htool *a = (Mat_Htool *)A->data;
334: PetscInt *idxc;
335: PetscBLASInt one = 1, bn;
337: PetscFunctionBegin;
338: if (nz) *nz = A->cmap->N;
339: if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */
340: PetscCall(PetscMalloc1(A->cmap->N, &idxc));
341: for (PetscInt i = 0; i < A->cmap->N; ++i) idxc[i] = i;
342: }
343: if (idx) *idx = idxc;
344: if (v) {
345: PetscCall(PetscMalloc1(A->cmap->N, v));
346: if (a->wrapper) a->wrapper->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
347: else reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx)->copy_submatrix(1, A->cmap->N, &row, idxc, *v);
348: PetscCall(PetscBLASIntCast(A->cmap->N, &bn));
349: PetscCallBLAS("BLASscal", BLASscal_(&bn, &a->s, *v, &one));
350: }
351: if (!idx) PetscCall(PetscFree(idxc));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode MatRestoreRow_Htool(Mat, PetscInt, PetscInt *, PetscInt **idx, PetscScalar **v)
356: {
357: PetscFunctionBegin;
358: if (idx) PetscCall(PetscFree(*idx));
359: if (v) PetscCall(PetscFree(*v));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode MatSetFromOptions_Htool(Mat A, PetscOptionItems *PetscOptionsObject)
364: {
365: Mat_Htool *a = (Mat_Htool *)A->data;
366: PetscInt n;
367: PetscBool flg;
369: PetscFunctionBegin;
370: PetscOptionsHeadBegin(PetscOptionsObject, "Htool options");
371: PetscCall(PetscOptionsInt("-mat_htool_min_cluster_size", "Minimal leaf size in cluster tree", nullptr, a->bs[0], a->bs, nullptr));
372: PetscCall(PetscOptionsInt("-mat_htool_max_block_size", "Maximal number of coefficients in a dense block", nullptr, a->bs[1], a->bs + 1, nullptr));
373: PetscCall(PetscOptionsReal("-mat_htool_epsilon", "Relative error in Frobenius norm when approximating a block", nullptr, a->epsilon, &a->epsilon, nullptr));
374: PetscCall(PetscOptionsReal("-mat_htool_eta", "Admissibility condition tolerance", nullptr, a->eta, &a->eta, nullptr));
375: PetscCall(PetscOptionsInt("-mat_htool_min_target_depth", "Minimal cluster tree depth associated with the rows", nullptr, a->depth[0], a->depth, nullptr));
376: PetscCall(PetscOptionsInt("-mat_htool_min_source_depth", "Minimal cluster tree depth associated with the columns", nullptr, a->depth[1], a->depth + 1, nullptr));
377: n = 0;
378: PetscCall(PetscOptionsEList("-mat_htool_compressor", "Type of compression", "MatHtoolCompressorType", MatHtoolCompressorTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes), MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA], &n, &flg));
379: if (flg) a->compressor = MatHtoolCompressorType(n);
380: n = 0;
381: PetscCall(PetscOptionsEList("-mat_htool_clustering", "Type of clustering", "MatHtoolClusteringType", MatHtoolClusteringTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes), MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR], &n, &flg));
382: if (flg) a->clustering = MatHtoolClusteringType(n);
383: PetscOptionsHeadEnd();
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: static PetscErrorCode MatAssemblyEnd_Htool(Mat A, MatAssemblyType)
388: {
389: Mat_Htool *a = (Mat_Htool *)A->data;
390: const PetscInt *ranges;
391: PetscInt *offset;
392: PetscMPIInt size;
393: char S = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 'H' : (A->symmetric == PETSC_BOOL3_TRUE ? 'S' : 'N'), uplo = S == 'N' ? 'N' : 'U';
394: htool::VirtualGenerator<PetscScalar> *generator = nullptr;
395: std::shared_ptr<htool::VirtualCluster> t, s = nullptr;
396: std::shared_ptr<htool::VirtualLowRankGenerator<PetscScalar>> compressor = nullptr;
398: PetscFunctionBegin;
399: PetscCall(PetscCitationsRegister(HtoolCitation, &HtoolCite));
400: delete a->wrapper;
401: delete a->hmatrix;
402: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
403: PetscCall(PetscMalloc1(2 * size, &offset));
404: PetscCall(MatGetOwnershipRanges(A, &ranges));
405: for (PetscInt i = 0; i < size; ++i) {
406: offset[2 * i] = ranges[i];
407: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
408: }
409: switch (a->clustering) {
410: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
411: t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
412: break;
413: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
414: t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
415: break;
416: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
417: t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
418: break;
419: default:
420: t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
421: }
422: t->set_minclustersize(a->bs[0]);
423: t->build(A->rmap->N, a->gcoords_target, offset, -1, PetscObjectComm((PetscObject)A));
424: if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N, A->cmap->N, a->dim, a->kernel, a->kernelctx);
425: else {
426: a->wrapper = nullptr;
427: generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx);
428: }
429: if (a->gcoords_target != a->gcoords_source) {
430: PetscCall(MatGetOwnershipRangesColumn(A, &ranges));
431: for (PetscInt i = 0; i < size; ++i) {
432: offset[2 * i] = ranges[i];
433: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
434: }
435: switch (a->clustering) {
436: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
437: s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
438: break;
439: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
440: s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
441: break;
442: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
443: s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
444: break;
445: default:
446: s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
447: }
448: s->set_minclustersize(a->bs[0]);
449: s->build(A->cmap->N, a->gcoords_source, offset, -1, PetscObjectComm((PetscObject)A));
450: S = uplo = 'N';
451: }
452: PetscCall(PetscFree(offset));
453: switch (a->compressor) {
454: case MAT_HTOOL_COMPRESSOR_FULL_ACA:
455: compressor = std::make_shared<htool::fullACA<PetscScalar>>();
456: break;
457: case MAT_HTOOL_COMPRESSOR_SVD:
458: compressor = std::make_shared<htool::SVD<PetscScalar>>();
459: break;
460: default:
461: compressor = std::make_shared<htool::sympartialACA<PetscScalar>>();
462: }
463: a->hmatrix = dynamic_cast<htool::VirtualHMatrix<PetscScalar> *>(new htool::HMatrix<PetscScalar>(t, s ? s : t, a->epsilon, a->eta, S, uplo, -1, PetscObjectComm((PetscObject)A)));
464: a->hmatrix->set_compression(compressor);
465: a->hmatrix->set_maxblocksize(a->bs[1]);
466: a->hmatrix->set_mintargetdepth(a->depth[0]);
467: a->hmatrix->set_minsourcedepth(a->depth[1]);
468: if (s) a->hmatrix->build(a->wrapper ? *a->wrapper : *generator, a->gcoords_target, a->gcoords_source);
469: else a->hmatrix->build(a->wrapper ? *a->wrapper : *generator, a->gcoords_target);
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: static PetscErrorCode MatProductNumeric_Htool(Mat C)
474: {
475: Mat_Product *product = C->product;
476: Mat_Htool *a = (Mat_Htool *)product->A->data;
477: const PetscScalar *in;
478: PetscScalar *out;
479: PetscInt N, lda;
481: PetscFunctionBegin;
482: MatCheckProduct(C, 1);
483: PetscCall(MatGetSize(C, nullptr, &N));
484: PetscCall(MatDenseGetLDA(C, &lda));
485: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
486: PetscCall(MatDenseGetArrayRead(product->B, &in));
487: PetscCall(MatDenseGetArrayWrite(C, &out));
488: switch (product->type) {
489: case MATPRODUCT_AB:
490: a->hmatrix->mvprod_local_to_local(in, out, N);
491: break;
492: case MATPRODUCT_AtB:
493: a->hmatrix->mvprod_transp_local_to_local(in, out, N);
494: break;
495: default:
496: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType %s is not supported", MatProductTypes[product->type]);
497: }
498: PetscCall(MatDenseRestoreArrayWrite(C, &out));
499: PetscCall(MatDenseRestoreArrayRead(product->B, &in));
500: PetscCall(MatScale(C, a->s));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
505: {
506: Mat_Product *product = C->product;
507: Mat A, B;
508: PetscBool flg;
510: PetscFunctionBegin;
511: MatCheckProduct(C, 1);
512: A = product->A;
513: B = product->B;
514: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, ""));
515: 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);
516: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
517: if (product->type == MATPRODUCT_AB) PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
518: else PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
519: }
520: PetscCall(MatSetType(C, MATDENSE));
521: PetscCall(MatSetUp(C));
522: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
523: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
524: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
525: C->ops->productsymbolic = nullptr;
526: C->ops->productnumeric = MatProductNumeric_Htool;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
531: {
532: PetscFunctionBegin;
533: MatCheckProduct(C, 1);
534: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A, const htool::VirtualHMatrix<PetscScalar> **hmatrix)
539: {
540: Mat_Htool *a = (Mat_Htool *)A->data;
542: PetscFunctionBegin;
543: *hmatrix = a->hmatrix;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@C
548: MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a `MATHTOOL`.
550: Input Parameter:
551: . A - hierarchical matrix
553: Output Parameter:
554: . hmatrix - opaque pointer to a Htool virtual matrix
556: Level: advanced
558: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`
559: @*/
560: PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A, const htool::VirtualHMatrix<PetscScalar> **hmatrix)
561: {
562: PetscFunctionBegin;
564: PetscAssertPointer(hmatrix, 2);
565: PetscTryMethod(A, "MatHtoolGetHierarchicalMat_C", (Mat, const htool::VirtualHMatrix<PetscScalar> **), (A, hmatrix));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
570: {
571: Mat_Htool *a = (Mat_Htool *)A->data;
573: PetscFunctionBegin;
574: a->kernel = kernel;
575: a->kernelctx = kernelctx;
576: delete a->wrapper;
577: if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N, A->cmap->N, a->dim, a->kernel, a->kernelctx);
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@C
582: MatHtoolSetKernel - Sets the kernel and context used for the assembly of a `MATHTOOL`.
584: Input Parameters:
585: + A - hierarchical matrix
586: . kernel - computational kernel (or `NULL`)
587: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
589: Level: advanced
591: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatCreateHtoolFromKernel()`
592: @*/
593: PETSC_EXTERN PetscErrorCode MatHtoolSetKernel(Mat A, MatHtoolKernelFn *kernel, void *kernelctx)
594: {
595: PetscFunctionBegin;
598: if (!kernel) PetscAssertPointer(kernelctx, 3);
599: PetscTryMethod(A, "MatHtoolSetKernel_C", (Mat, MatHtoolKernelFn *, void *), (A, kernel, kernelctx));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A, IS *is)
604: {
605: Mat_Htool *a = (Mat_Htool *)A->data;
606: std::vector<PetscInt> source;
608: PetscFunctionBegin;
609: source = a->hmatrix->get_source_cluster()->get_local_perm();
610: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), source.size(), source.data(), PETSC_COPY_VALUES, is));
611: PetscCall(ISSetPermutation(*is));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@C
616: MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster for a `MATHTOOL` matrix.
618: Input Parameter:
619: . A - hierarchical matrix
621: Output Parameter:
622: . is - permutation
624: Level: advanced
626: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationTarget()`, `MatHtoolUsePermutation()`
627: @*/
628: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationSource(Mat A, IS *is)
629: {
630: PetscFunctionBegin;
632: if (!is) PetscAssertPointer(is, 2);
633: PetscTryMethod(A, "MatHtoolGetPermutationSource_C", (Mat, IS *), (A, is));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A, IS *is)
638: {
639: Mat_Htool *a = (Mat_Htool *)A->data;
640: std::vector<PetscInt> target;
642: PetscFunctionBegin;
643: target = a->hmatrix->get_target_cluster()->get_local_perm();
644: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), target.size(), target.data(), PETSC_COPY_VALUES, is));
645: PetscCall(ISSetPermutation(*is));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@C
650: MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster for a `MATHTOOL` matrix.
652: Input Parameter:
653: . A - hierarchical matrix
655: Output Parameter:
656: . is - permutation
658: Level: advanced
660: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolUsePermutation()`
661: @*/
662: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationTarget(Mat A, IS *is)
663: {
664: PetscFunctionBegin;
666: if (!is) PetscAssertPointer(is, 2);
667: PetscTryMethod(A, "MatHtoolGetPermutationTarget_C", (Mat, IS *), (A, is));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A, PetscBool use)
672: {
673: Mat_Htool *a = (Mat_Htool *)A->data;
675: PetscFunctionBegin;
676: a->hmatrix->set_use_permutation(use);
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: /*@C
681: MatHtoolUsePermutation - Sets whether a `MATHTOOL` matrix should permute input (resp. output) vectors following its internal source (resp. target) permutation.
683: Input Parameters:
684: + A - hierarchical matrix
685: - use - Boolean value
687: Level: advanced
689: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolGetPermutationTarget()`
690: @*/
691: PETSC_EXTERN PetscErrorCode MatHtoolUsePermutation(Mat A, PetscBool use)
692: {
693: PetscFunctionBegin;
696: PetscTryMethod(A, "MatHtoolUsePermutation_C", (Mat, PetscBool), (A, use));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode MatConvert_Htool_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
701: {
702: Mat C;
703: Mat_Htool *a = (Mat_Htool *)A->data;
704: PetscInt lda;
705: PetscScalar *array;
707: PetscFunctionBegin;
708: if (reuse == MAT_REUSE_MATRIX) {
709: C = *B;
710: PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible dimensions");
711: PetscCall(MatDenseGetLDA(C, &lda));
712: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
713: } else {
714: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
715: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
716: PetscCall(MatSetType(C, MATDENSE));
717: PetscCall(MatSetUp(C));
718: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
719: }
720: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
721: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
722: PetscCall(MatDenseGetArrayWrite(C, &array));
723: a->hmatrix->copy_local_dense_perm(array);
724: PetscCall(MatDenseRestoreArrayWrite(C, &array));
725: PetscCall(MatScale(C, a->s));
726: if (reuse == MAT_INPLACE_MATRIX) {
727: PetscCall(MatHeaderReplace(A, &C));
728: } else *B = C;
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: static PetscErrorCode GenEntriesTranspose(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *rows, const PetscInt *cols, PetscScalar *ptr, void *ctx)
733: {
734: MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose *)ctx;
735: PetscScalar *tmp;
737: PetscFunctionBegin;
738: PetscCall(generator->kernel(sdim, N, M, cols, rows, ptr, generator->kernelctx));
739: PetscCall(PetscMalloc1(M * N, &tmp));
740: PetscCall(PetscArraycpy(tmp, ptr, M * N));
741: for (PetscInt i = 0; i < M; ++i) {
742: for (PetscInt j = 0; j < N; ++j) ptr[i + j * M] = tmp[j + i * N];
743: }
744: PetscCall(PetscFree(tmp));
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: /* naive implementation which keeps a reference to the original Mat */
749: static PetscErrorCode MatTranspose_Htool(Mat A, MatReuse reuse, Mat *B)
750: {
751: Mat C;
752: Mat_Htool *a = (Mat_Htool *)A->data, *c;
753: PetscInt M = A->rmap->N, N = A->cmap->N, m = A->rmap->n, n = A->cmap->n;
754: PetscContainer container;
755: MatHtoolKernelTranspose *kernelt;
757: PetscFunctionBegin;
758: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
759: PetscCheck(reuse != MAT_INPLACE_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatTranspose() with MAT_INPLACE_MATRIX not supported");
760: if (reuse == MAT_INITIAL_MATRIX) {
761: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
762: PetscCall(MatSetSizes(C, n, m, N, M));
763: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
764: PetscCall(MatSetUp(C));
765: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)C), &container));
766: PetscCall(PetscNew(&kernelt));
767: PetscCall(PetscContainerSetPointer(container, kernelt));
768: PetscCall(PetscObjectCompose((PetscObject)C, "KernelTranspose", (PetscObject)container));
769: } else {
770: C = *B;
771: PetscCall(PetscObjectQuery((PetscObject)C, "KernelTranspose", (PetscObject *)&container));
772: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatTranspose() with MAT_INITIAL_MATRIX first");
773: PetscCall(PetscContainerGetPointer(container, (void **)&kernelt));
774: }
775: c = (Mat_Htool *)C->data;
776: c->dim = a->dim;
777: c->s = a->s;
778: c->kernel = GenEntriesTranspose;
779: if (kernelt->A != A) {
780: PetscCall(MatDestroy(&kernelt->A));
781: kernelt->A = A;
782: PetscCall(PetscObjectReference((PetscObject)A));
783: }
784: kernelt->kernel = a->kernel;
785: kernelt->kernelctx = a->kernelctx;
786: c->kernelctx = kernelt;
787: if (reuse == MAT_INITIAL_MATRIX) {
788: PetscCall(PetscMalloc1(N * c->dim, &c->gcoords_target));
789: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_source, N * c->dim));
790: if (a->gcoords_target != a->gcoords_source) {
791: PetscCall(PetscMalloc1(M * c->dim, &c->gcoords_source));
792: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_target, M * c->dim));
793: } else c->gcoords_source = c->gcoords_target;
794: PetscCall(PetscCalloc2(M, &c->work_source, N, &c->work_target));
795: }
796: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
797: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
798: if (reuse == MAT_INITIAL_MATRIX) *B = C;
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: /*@C
803: MatCreateHtoolFromKernel - Creates a `MATHTOOL` from a user-supplied kernel.
805: Input Parameters:
806: + comm - MPI communicator
807: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
808: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
809: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
810: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
811: . spacedim - dimension of the space coordinates
812: . coords_target - coordinates of the target
813: . coords_source - coordinates of the source
814: . kernel - computational kernel (or `NULL`)
815: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
817: Output Parameter:
818: . B - matrix
820: Options Database Keys:
821: + -mat_htool_min_cluster_size <`PetscInt`> - minimal leaf size in cluster tree
822: . -mat_htool_max_block_size <`PetscInt`> - maximal number of coefficients in a dense block
823: . -mat_htool_epsilon <`PetscReal`> - relative error in Frobenius norm when approximating a block
824: . -mat_htool_eta <`PetscReal`> - admissibility condition tolerance
825: . -mat_htool_min_target_depth <`PetscInt`> - minimal cluster tree depth associated with the rows
826: . -mat_htool_min_source_depth <`PetscInt`> - minimal cluster tree depth associated with the columns
827: . -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression
828: - -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering
830: Level: intermediate
832: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATHTOOL`, `PCSetCoordinates()`, `MatHtoolSetKernel()`, `MatHtoolCompressorType`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
833: @*/
834: PetscErrorCode MatCreateHtoolFromKernel(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt spacedim, const PetscReal coords_target[], const PetscReal coords_source[], MatHtoolKernelFn *kernel, void *kernelctx, Mat *B)
835: {
836: Mat A;
837: Mat_Htool *a;
839: PetscFunctionBegin;
840: PetscCall(MatCreate(comm, &A));
842: PetscAssertPointer(coords_target, 7);
843: PetscAssertPointer(coords_source, 8);
845: if (!kernel) PetscAssertPointer(kernelctx, 10);
846: PetscCall(MatSetSizes(A, m, n, M, N));
847: PetscCall(MatSetType(A, MATHTOOL));
848: PetscCall(MatSetUp(A));
849: a = (Mat_Htool *)A->data;
850: a->dim = spacedim;
851: a->s = 1.0;
852: a->kernel = kernel;
853: a->kernelctx = kernelctx;
854: PetscCall(PetscCalloc1(A->rmap->N * spacedim, &a->gcoords_target));
855: PetscCall(PetscArraycpy(a->gcoords_target + A->rmap->rstart * spacedim, coords_target, A->rmap->n * spacedim));
856: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_target, A->rmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global target coordinates */
857: if (coords_target != coords_source) {
858: PetscCall(PetscCalloc1(A->cmap->N * spacedim, &a->gcoords_source));
859: PetscCall(PetscArraycpy(a->gcoords_source + A->cmap->rstart * spacedim, coords_source, A->cmap->n * spacedim));
860: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_source, A->cmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global source coordinates */
861: } else a->gcoords_source = a->gcoords_target;
862: PetscCall(PetscCalloc2(A->cmap->N, &a->work_source, A->rmap->N, &a->work_target));
863: *B = A;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*MC
868: MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.
870: Use `./configure --download-htool` to install PETSc to use Htool.
872: Options Database Key:
873: . -mat_type htool - matrix type to `MATHTOOL`
875: Level: beginner
877: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
878: M*/
879: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
880: {
881: Mat_Htool *a;
883: PetscFunctionBegin;
884: PetscCall(PetscNew(&a));
885: A->data = (void *)a;
886: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATHTOOL));
887: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
888: A->ops->getdiagonal = MatGetDiagonal_Htool;
889: A->ops->getdiagonalblock = MatGetDiagonalBlock_Htool;
890: A->ops->mult = MatMult_Htool;
891: A->ops->multadd = MatMultAdd_Htool;
892: A->ops->multtranspose = MatMultTranspose_Htool;
893: if (!PetscDefined(USE_COMPLEX)) A->ops->multhermitiantranspose = MatMultTranspose_Htool;
894: A->ops->increaseoverlap = MatIncreaseOverlap_Htool;
895: A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
896: A->ops->transpose = MatTranspose_Htool;
897: A->ops->destroy = MatDestroy_Htool;
898: A->ops->view = MatView_Htool;
899: A->ops->setfromoptions = MatSetFromOptions_Htool;
900: A->ops->scale = MatScale_Htool;
901: A->ops->getrow = MatGetRow_Htool;
902: A->ops->restorerow = MatRestoreRow_Htool;
903: A->ops->assemblyend = MatAssemblyEnd_Htool;
904: a->dim = 0;
905: a->gcoords_target = nullptr;
906: a->gcoords_source = nullptr;
907: a->s = 1.0;
908: a->bs[0] = 10;
909: a->bs[1] = 1000000;
910: a->epsilon = PetscSqrtReal(PETSC_SMALL);
911: a->eta = 10.0;
912: a->depth[0] = 0;
913: a->depth[1] = 0;
914: a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
915: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", MatProductSetFromOptions_Htool));
916: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", MatProductSetFromOptions_Htool));
917: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", MatConvert_Htool_Dense));
918: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", MatConvert_Htool_Dense));
919: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", MatHtoolGetHierarchicalMat_Htool));
920: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", MatHtoolSetKernel_Htool));
921: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", MatHtoolGetPermutationSource_Htool));
922: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", MatHtoolGetPermutationTarget_Htool));
923: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", MatHtoolUsePermutation_Htool));
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }