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