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 *nz, PetscInt **idx, PetscScalar **v)
356: {
357: PetscFunctionBegin;
358: if (nz) *nz = 0;
359: if (idx) PetscCall(PetscFree(*idx));
360: if (v) PetscCall(PetscFree(*v));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode MatSetFromOptions_Htool(Mat A, PetscOptionItems *PetscOptionsObject)
365: {
366: Mat_Htool *a = (Mat_Htool *)A->data;
367: PetscInt n;
368: PetscBool flg;
370: PetscFunctionBegin;
371: PetscOptionsHeadBegin(PetscOptionsObject, "Htool options");
372: PetscCall(PetscOptionsInt("-mat_htool_min_cluster_size", "Minimal leaf size in cluster tree", nullptr, a->bs[0], a->bs, nullptr));
373: PetscCall(PetscOptionsInt("-mat_htool_max_block_size", "Maximal number of coefficients in a dense block", nullptr, a->bs[1], a->bs + 1, nullptr));
374: PetscCall(PetscOptionsReal("-mat_htool_epsilon", "Relative error in Frobenius norm when approximating a block", nullptr, a->epsilon, &a->epsilon, nullptr));
375: PetscCall(PetscOptionsReal("-mat_htool_eta", "Admissibility condition tolerance", nullptr, a->eta, &a->eta, nullptr));
376: PetscCall(PetscOptionsInt("-mat_htool_min_target_depth", "Minimal cluster tree depth associated with the rows", nullptr, a->depth[0], a->depth, nullptr));
377: PetscCall(PetscOptionsInt("-mat_htool_min_source_depth", "Minimal cluster tree depth associated with the columns", nullptr, a->depth[1], a->depth + 1, nullptr));
378: n = 0;
379: PetscCall(PetscOptionsEList("-mat_htool_compressor", "Type of compression", "MatHtoolCompressorType", MatHtoolCompressorTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes), MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA], &n, &flg));
380: if (flg) a->compressor = MatHtoolCompressorType(n);
381: n = 0;
382: PetscCall(PetscOptionsEList("-mat_htool_clustering", "Type of clustering", "MatHtoolClusteringType", MatHtoolClusteringTypes, PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes), MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR], &n, &flg));
383: if (flg) a->clustering = MatHtoolClusteringType(n);
384: PetscOptionsHeadEnd();
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode MatAssemblyEnd_Htool(Mat A, MatAssemblyType)
389: {
390: Mat_Htool *a = (Mat_Htool *)A->data;
391: const PetscInt *ranges;
392: PetscInt *offset;
393: PetscMPIInt size;
394: char S = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 'H' : (A->symmetric == PETSC_BOOL3_TRUE ? 'S' : 'N'), uplo = S == 'N' ? 'N' : 'U';
395: htool::VirtualGenerator<PetscScalar> *generator = nullptr;
396: std::shared_ptr<htool::VirtualCluster> t, s = nullptr;
397: std::shared_ptr<htool::VirtualLowRankGenerator<PetscScalar>> compressor = nullptr;
399: PetscFunctionBegin;
400: PetscCall(PetscCitationsRegister(HtoolCitation, &HtoolCite));
401: delete a->wrapper;
402: delete a->hmatrix;
403: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
404: PetscCall(PetscMalloc1(2 * size, &offset));
405: PetscCall(MatGetOwnershipRanges(A, &ranges));
406: for (PetscInt i = 0; i < size; ++i) {
407: offset[2 * i] = ranges[i];
408: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
409: }
410: switch (a->clustering) {
411: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
412: t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
413: break;
414: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
415: t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
416: break;
417: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
418: t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
419: break;
420: default:
421: t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
422: }
423: t->set_minclustersize(a->bs[0]);
424: t->build(A->rmap->N, a->gcoords_target, offset, -1, PetscObjectComm((PetscObject)A));
425: if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N, A->cmap->N, a->dim, a->kernel, a->kernelctx);
426: else {
427: a->wrapper = nullptr;
428: generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar> *>(a->kernelctx);
429: }
430: if (a->gcoords_target != a->gcoords_source) {
431: PetscCall(MatGetOwnershipRangesColumn(A, &ranges));
432: for (PetscInt i = 0; i < size; ++i) {
433: offset[2 * i] = ranges[i];
434: offset[2 * i + 1] = ranges[i + 1] - ranges[i];
435: }
436: switch (a->clustering) {
437: case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
438: s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
439: break;
440: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
441: s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
442: break;
443: case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
444: s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
445: break;
446: default:
447: s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
448: }
449: s->set_minclustersize(a->bs[0]);
450: s->build(A->cmap->N, a->gcoords_source, offset, -1, PetscObjectComm((PetscObject)A));
451: S = uplo = 'N';
452: }
453: PetscCall(PetscFree(offset));
454: switch (a->compressor) {
455: case MAT_HTOOL_COMPRESSOR_FULL_ACA:
456: compressor = std::make_shared<htool::fullACA<PetscScalar>>();
457: break;
458: case MAT_HTOOL_COMPRESSOR_SVD:
459: compressor = std::make_shared<htool::SVD<PetscScalar>>();
460: break;
461: default:
462: compressor = std::make_shared<htool::sympartialACA<PetscScalar>>();
463: }
464: 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)));
465: a->hmatrix->set_compression(compressor);
466: a->hmatrix->set_maxblocksize(a->bs[1]);
467: a->hmatrix->set_mintargetdepth(a->depth[0]);
468: a->hmatrix->set_minsourcedepth(a->depth[1]);
469: if (s) a->hmatrix->build(a->wrapper ? *a->wrapper : *generator, a->gcoords_target, a->gcoords_source);
470: else a->hmatrix->build(a->wrapper ? *a->wrapper : *generator, a->gcoords_target);
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode MatProductNumeric_Htool(Mat C)
475: {
476: Mat_Product *product = C->product;
477: Mat_Htool *a = (Mat_Htool *)product->A->data;
478: const PetscScalar *in;
479: PetscScalar *out;
480: PetscInt N, lda;
482: PetscFunctionBegin;
483: MatCheckProduct(C, 1);
484: PetscCall(MatGetSize(C, nullptr, &N));
485: PetscCall(MatDenseGetLDA(C, &lda));
486: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
487: PetscCall(MatDenseGetArrayRead(product->B, &in));
488: PetscCall(MatDenseGetArrayWrite(C, &out));
489: switch (product->type) {
490: case MATPRODUCT_AB:
491: a->hmatrix->mvprod_local_to_local(in, out, N);
492: break;
493: case MATPRODUCT_AtB:
494: a->hmatrix->mvprod_transp_local_to_local(in, out, N);
495: break;
496: default:
497: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType %s is not supported", MatProductTypes[product->type]);
498: }
499: PetscCall(MatDenseRestoreArrayWrite(C, &out));
500: PetscCall(MatDenseRestoreArrayRead(product->B, &in));
501: PetscCall(MatScale(C, a->s));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
506: {
507: Mat_Product *product = C->product;
508: Mat A, B;
509: PetscBool flg;
511: PetscFunctionBegin;
512: MatCheckProduct(C, 1);
513: A = product->A;
514: B = product->B;
515: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, ""));
516: PetscCheck(flg, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "MatProduct_AB not supported for %s", ((PetscObject)product->B)->type_name);
517: switch (product->type) {
518: case MATPRODUCT_AB:
519: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
520: break;
521: case MATPRODUCT_AtB:
522: if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
523: break;
524: default:
525: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
526: }
527: PetscCall(MatSetType(C, MATDENSE));
528: PetscCall(MatSetUp(C));
529: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
530: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
531: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
532: C->ops->productsymbolic = nullptr;
533: C->ops->productnumeric = MatProductNumeric_Htool;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
538: {
539: PetscFunctionBegin;
540: MatCheckProduct(C, 1);
541: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A, const htool::VirtualHMatrix<PetscScalar> **hmatrix)
546: {
547: Mat_Htool *a = (Mat_Htool *)A->data;
549: PetscFunctionBegin;
550: *hmatrix = a->hmatrix;
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: /*@C
555: MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a `MATHTOOL`.
557: Input Parameter:
558: . A - hierarchical matrix
560: Output Parameter:
561: . hmatrix - opaque pointer to a Htool virtual matrix
563: Level: advanced
565: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`
566: @*/
567: PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A, const htool::VirtualHMatrix<PetscScalar> **hmatrix)
568: {
569: PetscFunctionBegin;
572: PetscTryMethod(A, "MatHtoolGetHierarchicalMat_C", (Mat, const htool::VirtualHMatrix<PetscScalar> **), (A, hmatrix));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A, MatHtoolKernel kernel, void *kernelctx)
577: {
578: Mat_Htool *a = (Mat_Htool *)A->data;
580: PetscFunctionBegin;
581: a->kernel = kernel;
582: a->kernelctx = kernelctx;
583: delete a->wrapper;
584: if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N, A->cmap->N, a->dim, a->kernel, a->kernelctx);
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: /*@C
589: MatHtoolSetKernel - Sets the kernel and context used for the assembly of a `MATHTOOL`.
591: Input Parameters:
592: + A - hierarchical matrix
593: . kernel - computational kernel (or `NULL`)
594: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
596: Level: advanced
598: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatCreateHtoolFromKernel()`
599: @*/
600: PETSC_EXTERN PetscErrorCode MatHtoolSetKernel(Mat A, MatHtoolKernel kernel, void *kernelctx)
601: {
602: PetscFunctionBegin;
606: PetscTryMethod(A, "MatHtoolSetKernel_C", (Mat, MatHtoolKernel, void *), (A, kernel, kernelctx));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A, IS *is)
611: {
612: Mat_Htool *a = (Mat_Htool *)A->data;
613: std::vector<PetscInt> source;
615: PetscFunctionBegin;
616: source = a->hmatrix->get_source_cluster()->get_local_perm();
617: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), source.size(), source.data(), PETSC_COPY_VALUES, is));
618: PetscCall(ISSetPermutation(*is));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /*@C
623: MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster for a `MATHTOOL` matrix.
625: Input Parameter:
626: . A - hierarchical matrix
628: Output Parameter:
629: . is - permutation
631: Level: advanced
633: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationTarget()`, `MatHtoolUsePermutation()`
634: @*/
635: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationSource(Mat A, IS *is)
636: {
637: PetscFunctionBegin;
640: PetscTryMethod(A, "MatHtoolGetPermutationSource_C", (Mat, IS *), (A, is));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A, IS *is)
645: {
646: Mat_Htool *a = (Mat_Htool *)A->data;
647: std::vector<PetscInt> target;
649: PetscFunctionBegin;
650: target = a->hmatrix->get_target_cluster()->get_local_perm();
651: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), target.size(), target.data(), PETSC_COPY_VALUES, is));
652: PetscCall(ISSetPermutation(*is));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*@C
657: MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster for a `MATHTOOL` matrix.
659: Input Parameter:
660: . A - hierarchical matrix
662: Output Parameter:
663: . is - permutation
665: Level: advanced
667: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolUsePermutation()`
668: @*/
669: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationTarget(Mat A, IS *is)
670: {
671: PetscFunctionBegin;
674: PetscTryMethod(A, "MatHtoolGetPermutationTarget_C", (Mat, IS *), (A, is));
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A, PetscBool use)
679: {
680: Mat_Htool *a = (Mat_Htool *)A->data;
682: PetscFunctionBegin;
683: a->hmatrix->set_use_permutation(use);
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@C
688: MatHtoolUsePermutation - Sets whether a `MATHTOOL` matrix should permute input (resp. output) vectors following its internal source (resp. target) permutation.
690: Input Parameters:
691: + A - hierarchical matrix
692: - use - Boolean value
694: Level: advanced
696: .seealso: [](ch_matrices), `Mat`, `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolGetPermutationTarget()`
697: @*/
698: PETSC_EXTERN PetscErrorCode MatHtoolUsePermutation(Mat A, PetscBool use)
699: {
700: PetscFunctionBegin;
703: PetscTryMethod(A, "MatHtoolUsePermutation_C", (Mat, PetscBool), (A, use));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: static PetscErrorCode MatConvert_Htool_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
708: {
709: Mat C;
710: Mat_Htool *a = (Mat_Htool *)A->data;
711: PetscInt lda;
712: PetscScalar *array;
714: PetscFunctionBegin;
715: if (reuse == MAT_REUSE_MATRIX) {
716: C = *B;
717: PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible dimensions");
718: PetscCall(MatDenseGetLDA(C, &lda));
719: PetscCheck(lda == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")", lda, C->rmap->n);
720: } else {
721: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
722: PetscCall(MatSetSizes(C, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
723: PetscCall(MatSetType(C, MATDENSE));
724: PetscCall(MatSetUp(C));
725: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
726: }
727: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
728: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
729: PetscCall(MatDenseGetArrayWrite(C, &array));
730: a->hmatrix->copy_local_dense_perm(array);
731: PetscCall(MatDenseRestoreArrayWrite(C, &array));
732: PetscCall(MatScale(C, a->s));
733: if (reuse == MAT_INPLACE_MATRIX) {
734: PetscCall(MatHeaderReplace(A, &C));
735: } else *B = C;
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode GenEntriesTranspose(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *rows, const PetscInt *cols, PetscScalar *ptr, void *ctx)
740: {
741: MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose *)ctx;
742: PetscScalar *tmp;
744: PetscFunctionBegin;
745: PetscCall(generator->kernel(sdim, N, M, cols, rows, ptr, generator->kernelctx));
746: PetscCall(PetscMalloc1(M * N, &tmp));
747: PetscCall(PetscArraycpy(tmp, ptr, M * N));
748: for (PetscInt i = 0; i < M; ++i) {
749: for (PetscInt j = 0; j < N; ++j) ptr[i + j * M] = tmp[j + i * N];
750: }
751: PetscCall(PetscFree(tmp));
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /* naive implementation which keeps a reference to the original Mat */
756: static PetscErrorCode MatTranspose_Htool(Mat A, MatReuse reuse, Mat *B)
757: {
758: Mat C;
759: Mat_Htool *a = (Mat_Htool *)A->data, *c;
760: PetscInt M = A->rmap->N, N = A->cmap->N, m = A->rmap->n, n = A->cmap->n;
761: PetscContainer container;
762: MatHtoolKernelTranspose *kernelt;
764: PetscFunctionBegin;
765: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
766: PetscCheck(reuse != MAT_INPLACE_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatTranspose() with MAT_INPLACE_MATRIX not supported");
767: if (reuse == MAT_INITIAL_MATRIX) {
768: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
769: PetscCall(MatSetSizes(C, n, m, N, M));
770: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
771: PetscCall(MatSetUp(C));
772: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)C), &container));
773: PetscCall(PetscNew(&kernelt));
774: PetscCall(PetscContainerSetPointer(container, kernelt));
775: PetscCall(PetscObjectCompose((PetscObject)C, "KernelTranspose", (PetscObject)container));
776: } else {
777: C = *B;
778: PetscCall(PetscObjectQuery((PetscObject)C, "KernelTranspose", (PetscObject *)&container));
779: PetscCheck(container, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatTranspose() with MAT_INITIAL_MATRIX first");
780: PetscCall(PetscContainerGetPointer(container, (void **)&kernelt));
781: }
782: c = (Mat_Htool *)C->data;
783: c->dim = a->dim;
784: c->s = a->s;
785: c->kernel = GenEntriesTranspose;
786: if (kernelt->A != A) {
787: PetscCall(MatDestroy(&kernelt->A));
788: kernelt->A = A;
789: PetscCall(PetscObjectReference((PetscObject)A));
790: }
791: kernelt->kernel = a->kernel;
792: kernelt->kernelctx = a->kernelctx;
793: c->kernelctx = kernelt;
794: if (reuse == MAT_INITIAL_MATRIX) {
795: PetscCall(PetscMalloc1(N * c->dim, &c->gcoords_target));
796: PetscCall(PetscArraycpy(c->gcoords_target, a->gcoords_source, N * c->dim));
797: if (a->gcoords_target != a->gcoords_source) {
798: PetscCall(PetscMalloc1(M * c->dim, &c->gcoords_source));
799: PetscCall(PetscArraycpy(c->gcoords_source, a->gcoords_target, M * c->dim));
800: } else c->gcoords_source = c->gcoords_target;
801: PetscCall(PetscCalloc2(M, &c->work_source, N, &c->work_target));
802: }
803: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
804: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
805: if (reuse == MAT_INITIAL_MATRIX) *B = C;
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }
809: /*@C
810: MatCreateHtoolFromKernel - Creates a `MATHTOOL` from a user-supplied kernel.
812: Input Parameters:
813: + comm - MPI communicator
814: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
815: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
816: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
817: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
818: . spacedim - dimension of the space coordinates
819: . coords_target - coordinates of the target
820: . coords_source - coordinates of the source
821: . kernel - computational kernel (or `NULL`)
822: - kernelctx - kernel context (if kernel is `NULL`, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)
824: Output Parameter:
825: . B - matrix
827: Options Database Keys:
828: + -mat_htool_min_cluster_size <`PetscInt`> - minimal leaf size in cluster tree
829: . -mat_htool_max_block_size <`PetscInt`> - maximal number of coefficients in a dense block
830: . -mat_htool_epsilon <`PetscReal`> - relative error in Frobenius norm when approximating a block
831: . -mat_htool_eta <`PetscReal`> - admissibility condition tolerance
832: . -mat_htool_min_target_depth <`PetscInt`> - minimal cluster tree depth associated with the rows
833: . -mat_htool_min_source_depth <`PetscInt`> - minimal cluster tree depth associated with the columns
834: . -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression
835: - -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering
837: Level: intermediate
839: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATHTOOL`, `PCSetCoordinates()`, `MatHtoolSetKernel()`, `MatHtoolCompressorType`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
840: @*/
841: PetscErrorCode MatCreateHtoolFromKernel(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt spacedim, const PetscReal coords_target[], const PetscReal coords_source[], MatHtoolKernel kernel, void *kernelctx, Mat *B)
842: {
843: Mat A;
844: Mat_Htool *a;
846: PetscFunctionBegin;
847: PetscCall(MatCreate(comm, &A));
853: PetscCall(MatSetSizes(A, m, n, M, N));
854: PetscCall(MatSetType(A, MATHTOOL));
855: PetscCall(MatSetUp(A));
856: a = (Mat_Htool *)A->data;
857: a->dim = spacedim;
858: a->s = 1.0;
859: a->kernel = kernel;
860: a->kernelctx = kernelctx;
861: PetscCall(PetscCalloc1(A->rmap->N * spacedim, &a->gcoords_target));
862: PetscCall(PetscArraycpy(a->gcoords_target + A->rmap->rstart * spacedim, coords_target, A->rmap->n * spacedim));
863: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_target, A->rmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global target coordinates */
864: if (coords_target != coords_source) {
865: PetscCall(PetscCalloc1(A->cmap->N * spacedim, &a->gcoords_source));
866: PetscCall(PetscArraycpy(a->gcoords_source + A->cmap->rstart * spacedim, coords_source, A->cmap->n * spacedim));
867: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, a->gcoords_source, A->cmap->N * spacedim, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)A))); /* global source coordinates */
868: } else a->gcoords_source = a->gcoords_target;
869: PetscCall(PetscCalloc2(A->cmap->N, &a->work_source, A->rmap->N, &a->work_target));
870: *B = A;
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: /*MC
875: MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.
877: Use `./configure --download-htool` to install PETSc to use Htool.
879: Options Database Key:
880: . -mat_type htool - matrix type to `MATHTOOL`
882: Level: beginner
884: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
885: M*/
886: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
887: {
888: Mat_Htool *a;
890: PetscFunctionBegin;
891: PetscCall(PetscNew(&a));
892: A->data = (void *)a;
893: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATHTOOL));
894: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
895: A->ops->getdiagonal = MatGetDiagonal_Htool;
896: A->ops->getdiagonalblock = MatGetDiagonalBlock_Htool;
897: A->ops->mult = MatMult_Htool;
898: A->ops->multadd = MatMultAdd_Htool;
899: A->ops->multtranspose = MatMultTranspose_Htool;
900: if (!PetscDefined(USE_COMPLEX)) A->ops->multhermitiantranspose = MatMultTranspose_Htool;
901: A->ops->increaseoverlap = MatIncreaseOverlap_Htool;
902: A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
903: A->ops->transpose = MatTranspose_Htool;
904: A->ops->destroy = MatDestroy_Htool;
905: A->ops->view = MatView_Htool;
906: A->ops->setfromoptions = MatSetFromOptions_Htool;
907: A->ops->scale = MatScale_Htool;
908: A->ops->getrow = MatGetRow_Htool;
909: A->ops->restorerow = MatRestoreRow_Htool;
910: A->ops->assemblyend = MatAssemblyEnd_Htool;
911: a->dim = 0;
912: a->gcoords_target = nullptr;
913: a->gcoords_source = nullptr;
914: a->s = 1.0;
915: a->bs[0] = 10;
916: a->bs[1] = 1000000;
917: a->epsilon = PetscSqrtReal(PETSC_SMALL);
918: a->eta = 10.0;
919: a->depth[0] = 0;
920: a->depth[1] = 0;
921: a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
922: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_seqdense_C", MatProductSetFromOptions_Htool));
923: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_htool_mpidense_C", MatProductSetFromOptions_Htool));
924: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_seqdense_C", MatConvert_Htool_Dense));
925: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_htool_mpidense_C", MatConvert_Htool_Dense));
926: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetHierarchicalMat_C", MatHtoolGetHierarchicalMat_Htool));
927: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolSetKernel_C", MatHtoolSetKernel_Htool));
928: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationSource_C", MatHtoolGetPermutationSource_Htool));
929: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolGetPermutationTarget_C", MatHtoolGetPermutationTarget_Htool));
930: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHtoolUsePermutation_C", MatHtoolUsePermutation_Htool));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }