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, &copy));
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, &copy));
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: }