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