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) {
157:     PetscCalloc1(n,submat);
158:   }
159:   for (i=0; i<n; ++i) {
160:     ISGetLocalSize(irow[i],&nrow);
161:     ISGetLocalSize(icol[i],&m);
162:     ISGetIndices(irow[i],&idxr);
163:     ISGetIndices(icol[i],&idxc);
164:     if (scall != MAT_REUSE_MATRIX) {
165:       MatCreateDense(PETSC_COMM_SELF,nrow,m,nrow,m,NULL,(*submat)+i);
166:     }
167:     MatDenseGetArrayWrite((*submat)[i],&ptr);
168:     if (irow[i] == icol[i]) { /* same row and column IS? */
169:       MatHasCongruentLayouts(A,&flg);
170:       if (flg) {
171:         ISSorted(irow[i],&flg);
172:         if (flg) { /* sorted IS? */
173:           it = std::lower_bound(idxr,idxr+nrow,A->rmap->rstart);
174:           if (it != idxr+nrow && *it == A->rmap->rstart) { /* rmap->rstart in IS? */
175:             if (std::distance(idxr,it) + A->rmap->n <= nrow) { /* long enough IS to store the local diagonal block? */
176:               for (PetscInt j=0; j<A->rmap->n && flg; ++j) if (PetscUnlikely(it[j] != A->rmap->rstart+j)) flg = PETSC_FALSE;
177:               if (flg) { /* complete local diagonal block in IS? */
178:                 /* fast extraction when the local diagonal block is part of the submatrix, e.g., for PCASM or PCHPDDM
179:                  *      [   B   C   E   ]
180:                  *  A = [   B   D   E   ]
181:                  *      [   B   F   E   ]
182:                  */
183:                 m = std::distance(idxr,it); /* shift of the coefficient (0,0) of block D from above */
184:                 MatGetDiagonalBlock_Htool(A,&D);
185:                 MatDenseGetArrayRead(D,&copy);
186:                 for (PetscInt k=0; k<A->rmap->n; ++k) {
187:                   PetscArraycpy(ptr+(m+k)*nrow+m,copy+k*A->rmap->n,A->rmap->n); /* block D from above */
188:                 }
189:                 MatDenseRestoreArrayRead(D,&copy);
190:                 if (m) {
191:                   a->wrapper->copy_submatrix(nrow,m,idxr,idxc,ptr); /* vertical block B from above */
192:                   /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
193:                   if (A->symmetric || A->hermitian) {
194:                     MatCreateDense(PETSC_COMM_SELF,A->rmap->n,m,A->rmap->n,m,ptr+m,&B);
195:                     MatDenseSetLDA(B,nrow);
196:                     MatCreateDense(PETSC_COMM_SELF,m,A->rmap->n,m,A->rmap->n,ptr+m*nrow,&BT);
197:                     MatDenseSetLDA(BT,nrow);
198:                     if (A->hermitian && PetscDefined(USE_COMPLEX)) {
199:                       MatHermitianTranspose(B,MAT_REUSE_MATRIX,&BT);
200:                     } else {
201:                       MatTranspose(B,MAT_REUSE_MATRIX,&BT);
202:                     }
203:                     MatDestroy(&B);
204:                     MatDestroy(&BT);
205:                   } else {
206:                     for (PetscInt k=0; k<A->rmap->n; ++k) { /* block C from above */
207:                       a->wrapper->copy_submatrix(m,1,idxr,idxc+m+k,ptr+(m+k)*nrow);
208:                     }
209:                   }
210:                 }
211:                 if (m+A->rmap->n != nrow) {
212:                   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 */
213:                   /* entry-wise assembly may be costly, so transpose already-computed entries when possible */
214:                   if (A->symmetric || A->hermitian) {
215:                     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);
216:                     MatDenseSetLDA(B,nrow);
217:                     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);
218:                     MatDenseSetLDA(BT,nrow);
219:                     if (A->hermitian && PetscDefined(USE_COMPLEX)) {
220:                       MatHermitianTranspose(B,MAT_REUSE_MATRIX,&BT);
221:                     } else {
222:                       MatTranspose(B,MAT_REUSE_MATRIX,&BT);
223:                     }
224:                     MatDestroy(&B);
225:                     MatDestroy(&BT);
226:                   } else {
227:                     for (PetscInt k=0; k<A->rmap->n; ++k) { /* block F from above */
228:                       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);
229:                     }
230:                   }
231:                 }
232:               } /* complete local diagonal block not in IS */
233:             } else flg = PETSC_FALSE; /* IS not long enough to store the local diagonal block */
234:           } else flg = PETSC_FALSE; /* rmap->rstart not in IS */
235:         } /* unsorted IS */
236:       }
237:     } else flg = PETSC_FALSE; /* different row and column IS */
238:     if (!flg) a->wrapper->copy_submatrix(nrow,m,idxr,idxc,ptr); /* reassemble everything */
239:     ISRestoreIndices(irow[i],&idxr);
240:     ISRestoreIndices(icol[i],&idxc);
241:     MatDenseRestoreArrayWrite((*submat)[i],&ptr);
242:     MatScale((*submat)[i],a->s);
243:   }
244:   return 0;
245: }

247: static PetscErrorCode MatDestroy_Htool(Mat A)
248: {
249:   Mat_Htool               *a = (Mat_Htool*)A->data;
250:   PetscContainer          container;
251:   MatHtoolKernelTranspose *kernelt;

253:   PetscObjectChangeTypeName((PetscObject)A,NULL);
254:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_seqdense_C",NULL);
255:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_mpidense_C",NULL);
256:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_seqdense_C",NULL);
257:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_mpidense_C",NULL);
258:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetHierarchicalMat_C",NULL);
259:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolSetKernel_C",NULL);
260:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationSource_C",NULL);
261:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationTarget_C",NULL);
262:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolUsePermutation_C",NULL);
263:   PetscObjectQuery((PetscObject)A,"KernelTranspose",(PetscObject*)&container);
264:   if (container) { /* created in MatTranspose_Htool() */
265:     PetscContainerGetPointer(container,(void**)&kernelt);
266:     MatDestroy(&kernelt->A);
267:     PetscFree(kernelt);
268:     PetscContainerDestroy(&container);
269:     PetscObjectCompose((PetscObject)A,"KernelTranspose",NULL);
270:   }
271:   if (a->gcoords_source != a->gcoords_target) {
272:     PetscFree(a->gcoords_source);
273:   }
274:   PetscFree(a->gcoords_target);
275:   PetscFree2(a->work_source,a->work_target);
276:   delete a->wrapper;
277:   delete a->hmatrix;
278:   PetscFree(A->data);
279:   return 0;
280: }

282: static PetscErrorCode MatView_Htool(Mat A,PetscViewer pv)
283: {
284:   Mat_Htool      *a = (Mat_Htool*)A->data;
285:   PetscBool      flg;

287:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&flg);
288:   if (flg) {
289:     PetscViewerASCIIPrintf(pv,"symmetry: %c\n",a->hmatrix->get_symmetry_type());
290:     if (PetscAbsScalar(a->s-1.0) > PETSC_MACHINE_EPSILON) {
291: #if defined(PETSC_USE_COMPLEX)
292:       PetscViewerASCIIPrintf(pv,"scaling: %g+%gi\n",(double)PetscRealPart(a->s),(double)PetscImaginaryPart(a->s));
293: #else
294:       PetscViewerASCIIPrintf(pv,"scaling: %g\n",(double)a->s);
295: #endif
296:     }
297:     PetscViewerASCIIPrintf(pv,"minimum cluster size: %" PetscInt_FMT "\n",a->bs[0]);
298:     PetscViewerASCIIPrintf(pv,"maximum block size: %" PetscInt_FMT "\n",a->bs[1]);
299:     PetscViewerASCIIPrintf(pv,"epsilon: %g\n",(double)a->epsilon);
300:     PetscViewerASCIIPrintf(pv,"eta: %g\n",(double)a->eta);
301:     PetscViewerASCIIPrintf(pv,"minimum target depth: %" PetscInt_FMT "\n",a->depth[0]);
302:     PetscViewerASCIIPrintf(pv,"minimum source depth: %" PetscInt_FMT "\n",a->depth[1]);
303:     PetscViewerASCIIPrintf(pv,"compressor: %s\n",MatHtoolCompressorTypes[a->compressor]);
304:     PetscViewerASCIIPrintf(pv,"clustering: %s\n",MatHtoolClusteringTypes[a->clustering]);
305:     PetscViewerASCIIPrintf(pv,"compression ratio: %s\n",a->hmatrix->get_infos("Compression_ratio").c_str());
306:     PetscViewerASCIIPrintf(pv,"space saving: %s\n",a->hmatrix->get_infos("Space_saving").c_str());
307:     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());
308:     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(),a->hmatrix->get_infos("Dense_block_size_max").c_str());
309:     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(),a->hmatrix->get_infos("Low_rank_block_size_max").c_str());
310:     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());
311:   }
312:   return 0;
313: }

315: static PetscErrorCode MatScale_Htool(Mat A,PetscScalar s)
316: {
317:   Mat_Htool *a = (Mat_Htool*)A->data;

319:   a->s *= s;
320:   return 0;
321: }

323: /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */
324: static PetscErrorCode MatGetRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
325: {
326:   Mat_Htool      *a = (Mat_Htool*)A->data;
327:   PetscInt       *idxc;
328:   PetscBLASInt   one = 1,bn;

330:   if (nz) *nz = A->cmap->N;
331:   if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */
332:     PetscMalloc1(A->cmap->N,&idxc);
333:     for (PetscInt i=0; i<A->cmap->N; ++i) idxc[i] = i;
334:   }
335:   if (idx) *idx = idxc;
336:   if (v) {
337:     PetscMalloc1(A->cmap->N,v);
338:     if (a->wrapper) a->wrapper->copy_submatrix(1,A->cmap->N,&row,idxc,*v);
339:     else reinterpret_cast<htool::VirtualGenerator<PetscScalar>*>(a->kernelctx)->copy_submatrix(1,A->cmap->N,&row,idxc,*v);
340:     PetscBLASIntCast(A->cmap->N,&bn);
341:     PetscStackCallBLAS("BLASscal",BLASscal_(&bn,&a->s,*v,&one));
342:   }
343:   if (!idx) {
344:     PetscFree(idxc);
345:   }
346:   return 0;
347: }

349: static PetscErrorCode MatRestoreRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
350: {
351:   if (nz) *nz = 0;
352:   if (idx) {
353:     PetscFree(*idx);
354:   }
355:   if (v) {
356:     PetscFree(*v);
357:   }
358:   return 0;
359: }

361: static PetscErrorCode MatSetFromOptions_Htool(PetscOptionItems *PetscOptionsObject,Mat A)
362: {
363:   Mat_Htool      *a = (Mat_Htool*)A->data;
364:   PetscInt       n;
365:   PetscBool      flg;

367:   PetscOptionsHeadBegin(PetscOptionsObject,"Htool options");
368:   PetscOptionsInt("-mat_htool_min_cluster_size","Minimal leaf size in cluster tree",NULL,a->bs[0],a->bs,NULL);
369:   PetscOptionsInt("-mat_htool_max_block_size","Maximal number of coefficients in a dense block",NULL,a->bs[1],a->bs + 1,NULL);
370:   PetscOptionsReal("-mat_htool_epsilon","Relative error in Frobenius norm when approximating a block",NULL,a->epsilon,&a->epsilon,NULL);
371:   PetscOptionsReal("-mat_htool_eta","Admissibility condition tolerance",NULL,a->eta,&a->eta,NULL);
372:   PetscOptionsInt("-mat_htool_min_target_depth","Minimal cluster tree depth associated with the rows",NULL,a->depth[0],a->depth,NULL);
373:   PetscOptionsInt("-mat_htool_min_source_depth","Minimal cluster tree depth associated with the columns",NULL,a->depth[1],a->depth + 1,NULL);
374:   n = 0;
375:   PetscOptionsEList("-mat_htool_compressor","Type of compression","MatHtoolCompressorType",MatHtoolCompressorTypes,PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes),MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA],&n,&flg);
376:   if (flg) a->compressor = MatHtoolCompressorType(n);
377:   n = 0;
378:   PetscOptionsEList("-mat_htool_clustering","Type of clustering","MatHtoolClusteringType",MatHtoolClusteringTypes,PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes),MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR],&n,&flg);
379:   if (flg) a->clustering = MatHtoolClusteringType(n);
380:   PetscOptionsHeadEnd();
381:   return 0;
382: }

384: static PetscErrorCode MatAssemblyEnd_Htool(Mat A,MatAssemblyType type)
385: {
386:   Mat_Htool                                                    *a = (Mat_Htool*)A->data;
387:   const PetscInt                                               *ranges;
388:   PetscInt                                                     *offset;
389:   PetscMPIInt                                                  size;
390:   char                                                         S = PetscDefined(USE_COMPLEX) && A->hermitian ? 'H' : (A->symmetric ? 'S' : 'N'),uplo = S == 'N' ? 'N' : 'U';
391:   htool::VirtualGenerator<PetscScalar>                         *generator = nullptr;
392:   std::shared_ptr<htool::VirtualCluster>                       t,s = nullptr;
393:   std::shared_ptr<htool::VirtualLowRankGenerator<PetscScalar>> compressor = nullptr;

395:   PetscCitationsRegister(HtoolCitation,&HtoolCite);
396:   delete a->wrapper;
397:   delete a->hmatrix;
398:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
399:   PetscMalloc1(2*size,&offset);
400:   MatGetOwnershipRanges(A,&ranges);
401:   for (PetscInt i=0; i<size; ++i) {
402:     offset[2*i] = ranges[i];
403:     offset[2*i+1] = ranges[i+1] - ranges[i];
404:   }
405:   switch (a->clustering) {
406:   case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
407:     t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
408:     break;
409:   case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
410:     t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
411:     break;
412:   case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
413:     t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
414:     break;
415:   default:
416:     t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
417:   }
418:   t->set_minclustersize(a->bs[0]);
419:   t->build(A->rmap->N,a->gcoords_target,offset);
420:   if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N,A->cmap->N,a->dim,a->kernel,a->kernelctx);
421:   else {
422:     a->wrapper = NULL;
423:     generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar>*>(a->kernelctx);
424:   }
425:   if (a->gcoords_target != a->gcoords_source) {
426:     MatGetOwnershipRangesColumn(A,&ranges);
427:     for (PetscInt i=0; i<size; ++i) {
428:       offset[2*i] = ranges[i];
429:       offset[2*i+1] = ranges[i+1] - ranges[i];
430:     }
431:     switch (a->clustering) {
432:     case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC:
433:       s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
434:       break;
435:     case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC:
436:       s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim);
437:       break;
438:     case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR:
439:       s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim);
440:       break;
441:     default:
442:       s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim);
443:     }
444:     s->set_minclustersize(a->bs[0]);
445:     s->build(A->cmap->N,a->gcoords_source,offset);
446:     S = uplo = 'N';
447:   }
448:   PetscFree(offset);
449:   switch (a->compressor) {
450:   case MAT_HTOOL_COMPRESSOR_FULL_ACA:
451:     compressor = std::make_shared<htool::fullACA<PetscScalar>>();
452:     break;
453:   case MAT_HTOOL_COMPRESSOR_SVD:
454:     compressor = std::make_shared<htool::SVD<PetscScalar>>();
455:     break;
456:   default:
457:     compressor = std::make_shared<htool::sympartialACA<PetscScalar>>();
458:   }
459:   a->hmatrix = dynamic_cast<htool::VirtualHMatrix<PetscScalar>*>(new htool::HMatrix<PetscScalar>(t,s ? s : t,a->epsilon,a->eta,S,uplo));
460:   a->hmatrix->set_compression(compressor);
461:   a->hmatrix->set_maxblocksize(a->bs[1]);
462:   a->hmatrix->set_mintargetdepth(a->depth[0]);
463:   a->hmatrix->set_minsourcedepth(a->depth[1]);
464:   if (s) a->hmatrix->build(a->wrapper ? *a->wrapper : *generator,a->gcoords_target,a->gcoords_source);
465:   else   a->hmatrix->build(a->wrapper ? *a->wrapper : *generator,a->gcoords_target);
466:   return 0;
467: }

469: static PetscErrorCode MatProductNumeric_Htool(Mat C)
470: {
471:   Mat_Product       *product = C->product;
472:   Mat_Htool         *a = (Mat_Htool*)product->A->data;
473:   const PetscScalar *in;
474:   PetscScalar       *out;
475:   PetscInt          N,lda;

477:   MatCheckProduct(C,1);
478:   MatGetSize(C,NULL,&N);
479:   MatDenseGetLDA(C,&lda);
481:   MatDenseGetArrayRead(product->B,&in);
482:   MatDenseGetArrayWrite(C,&out);
483:   switch (product->type) {
484:   case MATPRODUCT_AB:
485:     a->hmatrix->mvprod_local_to_local(in,out,N);
486:     break;
487:   case MATPRODUCT_AtB:
488:     a->hmatrix->mvprod_transp_local_to_local(in,out,N);
489:     break;
490:   default:
491:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductType %s is not supported",MatProductTypes[product->type]);
492:   }
493:   MatDenseRestoreArrayWrite(C,&out);
494:   MatDenseRestoreArrayRead(product->B,&in);
495:   MatScale(C,a->s);
496:   return 0;
497: }

499: static PetscErrorCode MatProductSymbolic_Htool(Mat C)
500: {
501:   Mat_Product    *product = C->product;
502:   Mat            A,B;
503:   PetscBool      flg;

505:   MatCheckProduct(C,1);
506:   A = product->A;
507:   B = product->B;
508:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,"");
510:   switch (product->type) {
511:   case MATPRODUCT_AB:
512:     if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
513:       MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
514:     }
515:     break;
516:   case MATPRODUCT_AtB:
517:     if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) {
518:       MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
519:     }
520:     break;
521:   default:
522:     SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]);
523:   }
524:   MatSetType(C,MATDENSE);
525:   MatSetUp(C);
526:   MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
527:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
528:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
529:   C->ops->productsymbolic = NULL;
530:   C->ops->productnumeric = MatProductNumeric_Htool;
531:   return 0;
532: }

534: static PetscErrorCode MatProductSetFromOptions_Htool(Mat C)
535: {
536:   MatCheckProduct(C,1);
537:   if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool;
538:   return 0;
539: }

541: static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix)
542: {
543:   Mat_Htool *a = (Mat_Htool*)A->data;

545:   *hmatrix = a->hmatrix;
546:   return 0;
547: }

549: /*@C
550:      MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a MATHTOOL.

552:    Input Parameter:
553: .     A - hierarchical matrix

555:    Output Parameter:
556: .     hmatrix - opaque pointer to a Htool virtual matrix

558:    Level: advanced

560: .seealso: `MATHTOOL`
561: @*/
562: PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix)
563: {
566:   PetscTryMethod(A,"MatHtoolGetHierarchicalMat_C",(Mat,const htool::VirtualHMatrix<PetscScalar>**),(A,hmatrix));
567:   return 0;
568: }

570: static PetscErrorCode MatHtoolSetKernel_Htool(Mat A,MatHtoolKernel kernel,void *kernelctx)
571: {
572:   Mat_Htool *a = (Mat_Htool*)A->data;

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:   return 0;
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: `MATHTOOL`, `MatCreateHtoolFromKernel()`
592: @*/
593: PETSC_EXTERN PetscErrorCode MatHtoolSetKernel(Mat A,MatHtoolKernel kernel,void *kernelctx)
594: {
598:   PetscTryMethod(A,"MatHtoolSetKernel_C",(Mat,MatHtoolKernel,void*),(A,kernel,kernelctx));
599:   return 0;
600: }

602: static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A,IS* is)
603: {
604:   Mat_Htool             *a = (Mat_Htool*)A->data;
605:   std::vector<PetscInt> source;

607:   source = a->hmatrix->get_source_cluster()->get_local_perm();
608:   ISCreateGeneral(PetscObjectComm((PetscObject)A),source.size(),source.data(),PETSC_COPY_VALUES,is);
609:   ISSetPermutation(*is);
610:   return 0;
611: }

613: /*@C
614:      MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster.

616:    Input Parameter:
617: .     A - hierarchical matrix

619:    Output Parameter:
620: .     is - permutation

622:    Level: advanced

624: .seealso: `MATHTOOL`, `MatHtoolGetPermutationTarget()`, `MatHtoolUsePermutation()`
625: @*/
626: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationSource(Mat A,IS* is)
627: {
630:   PetscTryMethod(A,"MatHtoolGetPermutationSource_C",(Mat,IS*),(A,is));
631:   return 0;
632: }

634: static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A,IS* is)
635: {
636:   Mat_Htool             *a = (Mat_Htool*)A->data;
637:   std::vector<PetscInt> target;

639:   target = a->hmatrix->get_target_cluster()->get_local_perm();
640:   ISCreateGeneral(PetscObjectComm((PetscObject)A),target.size(),target.data(),PETSC_COPY_VALUES,is);
641:   ISSetPermutation(*is);
642:   return 0;
643: }

645: /*@C
646:      MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster.

648:    Input Parameter:
649: .     A - hierarchical matrix

651:    Output Parameter:
652: .     is - permutation

654:    Level: advanced

656: .seealso: `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolUsePermutation()`
657: @*/
658: PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationTarget(Mat A,IS* is)
659: {
662:   PetscTryMethod(A,"MatHtoolGetPermutationTarget_C",(Mat,IS*),(A,is));
663:   return 0;
664: }

666: static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A,PetscBool use)
667: {
668:   Mat_Htool *a = (Mat_Htool*)A->data;

670:   a->hmatrix->set_use_permutation(use);
671:   return 0;
672: }

674: /*@C
675:      MatHtoolUsePermutation - Sets whether MATHTOOL should permute input (resp. output) vectors following its internal source (resp. target) permutation.

677:    Input Parameters:
678: +     A - hierarchical matrix
679: -     use - Boolean value

681:    Level: advanced

683: .seealso: `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolGetPermutationTarget()`
684: @*/
685: PETSC_EXTERN PetscErrorCode MatHtoolUsePermutation(Mat A,PetscBool use)
686: {
689:   PetscTryMethod(A,"MatHtoolUsePermutation_C",(Mat,PetscBool),(A,use));
690:   return 0;
691: }

693: static PetscErrorCode MatConvert_Htool_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B)
694: {
695:   Mat            C;
696:   Mat_Htool      *a = (Mat_Htool*)A->data;
697:   PetscInt       lda;
698:   PetscScalar    *array;

700:   if (reuse == MAT_REUSE_MATRIX) {
701:     C = *B;
703:     MatDenseGetLDA(C,&lda);
705:   } else {
706:     MatCreate(PetscObjectComm((PetscObject)A),&C);
707:     MatSetSizes(C,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
708:     MatSetType(C,MATDENSE);
709:     MatSetUp(C);
710:     MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
711:   }
712:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
713:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
714:   MatDenseGetArrayWrite(C,&array);
715:   a->hmatrix->copy_local_dense_perm(array);
716:   MatDenseRestoreArrayWrite(C,&array);
717:   MatScale(C,a->s);
718:   if (reuse == MAT_INPLACE_MATRIX) {
719:     MatHeaderReplace(A,&C);
720:   } else *B = C;
721:   return 0;
722: }

724: static PetscErrorCode GenEntriesTranspose(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *rows,const PetscInt *cols,PetscScalar *ptr,void *ctx)
725: {
726:   MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose*)ctx;
727:   PetscScalar             *tmp;

729:   generator->kernel(sdim,N,M,cols,rows,ptr,generator->kernelctx);
730:   PetscMalloc1(M*N,&tmp);
731:   PetscArraycpy(tmp,ptr,M*N);
732:   for (PetscInt i=0; i<M; ++i) {
733:     for (PetscInt j=0; j<N; ++j) ptr[i+j*M] = tmp[j+i*N];
734:   }
735:   PetscFree(tmp);
736:   return 0;
737: }

739: /* naive implementation which keeps a reference to the original Mat */
740: static PetscErrorCode MatTranspose_Htool(Mat A,MatReuse reuse,Mat *B)
741: {
742:   Mat                     C;
743:   Mat_Htool               *a = (Mat_Htool*)A->data,*c;
744:   PetscInt                M = A->rmap->N,N = A->cmap->N,m = A->rmap->n,n = A->cmap->n;
745:   PetscContainer          container;
746:   MatHtoolKernelTranspose *kernelt;

749:   if (reuse == MAT_INITIAL_MATRIX) {
750:     MatCreate(PetscObjectComm((PetscObject)A),&C);
751:     MatSetSizes(C,n,m,N,M);
752:     MatSetType(C,((PetscObject)A)->type_name);
753:     MatSetUp(C);
754:     PetscContainerCreate(PetscObjectComm((PetscObject)C),&container);
755:     PetscNew(&kernelt);
756:     PetscContainerSetPointer(container,kernelt);
757:     PetscObjectCompose((PetscObject)C,"KernelTranspose",(PetscObject)container);
758:   } else {
759:     C = *B;
760:     PetscObjectQuery((PetscObject)C,"KernelTranspose",(PetscObject*)&container);
762:     PetscContainerGetPointer(container,(void**)&kernelt);
763:   }
764:   c                  = (Mat_Htool*)C->data;
765:   c->dim             = a->dim;
766:   c->s               = a->s;
767:   c->kernel          = GenEntriesTranspose;
768:   if (kernelt->A != A) {
769:     MatDestroy(&kernelt->A);
770:     kernelt->A       = A;
771:     PetscObjectReference((PetscObject)A);
772:   }
773:   kernelt->kernel    = a->kernel;
774:   kernelt->kernelctx = a->kernelctx;
775:   c->kernelctx       = kernelt;
776:   if (reuse == MAT_INITIAL_MATRIX) {
777:     PetscMalloc1(N*c->dim,&c->gcoords_target);
778:     PetscArraycpy(c->gcoords_target,a->gcoords_source,N*c->dim);
779:     if (a->gcoords_target != a->gcoords_source) {
780:       PetscMalloc1(M*c->dim,&c->gcoords_source);
781:       PetscArraycpy(c->gcoords_source,a->gcoords_target,M*c->dim);
782:     } else c->gcoords_source = c->gcoords_target;
783:     PetscCalloc2(M,&c->work_source,N,&c->work_target);
784:   }
785:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
786:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
787:   if (reuse == MAT_INITIAL_MATRIX) *B = C;
788:   return 0;
789: }

791: /*@C
792:      MatCreateHtoolFromKernel - Creates a MATHTOOL from a user-supplied kernel.

794:    Input Parameters:
795: +     comm - MPI communicator
796: .     m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
797: .     n - number of local columns (or PETSC_DECIDE to have calculated if N is given)
798: .     M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
799: .     N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
800: .     spacedim - dimension of the space coordinates
801: .     coords_target - coordinates of the target
802: .     coords_source - coordinates of the source
803: .     kernel - computational kernel (or NULL)
804: -     kernelctx - kernel context (if kernel is NULL, the pointer must be of type htool::VirtualGenerator<PetscScalar>*)

806:    Output Parameter:
807: .     B - matrix

809:    Options Database Keys:
810: +     -mat_htool_min_cluster_size <PetscInt> - minimal leaf size in cluster tree
811: .     -mat_htool_max_block_size <PetscInt> - maximal number of coefficients in a dense block
812: .     -mat_htool_epsilon <PetscReal> - relative error in Frobenius norm when approximating a block
813: .     -mat_htool_eta <PetscReal> - admissibility condition tolerance
814: .     -mat_htool_min_target_depth <PetscInt> - minimal cluster tree depth associated with the rows
815: .     -mat_htool_min_source_depth <PetscInt> - minimal cluster tree depth associated with the columns
816: .     -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression
817: -     -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering

819:    Level: intermediate

821: .seealso: `MatCreate()`, `MATHTOOL`, `PCSetCoordinates()`, `MatHtoolSetKernel()`, `MatHtoolCompressorType`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
822: @*/
823: 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)
824: {
825:   Mat            A;
826:   Mat_Htool      *a;

828:   MatCreate(comm,&A);
834:   MatSetSizes(A,m,n,M,N);
835:   MatSetType(A,MATHTOOL);
836:   MatSetUp(A);
837:   a            = (Mat_Htool*)A->data;
838:   a->dim       = spacedim;
839:   a->s         = 1.0;
840:   a->kernel    = kernel;
841:   a->kernelctx = kernelctx;
842:   PetscCalloc1(A->rmap->N*spacedim,&a->gcoords_target);
843:   PetscArraycpy(a->gcoords_target+A->rmap->rstart*spacedim,coords_target,A->rmap->n*spacedim);
844:   MPIU_Allreduce(MPI_IN_PLACE,a->gcoords_target,A->rmap->N*spacedim,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)A)); /* global target coordinates */
845:   if (coords_target != coords_source) {
846:     PetscCalloc1(A->cmap->N*spacedim,&a->gcoords_source);
847:     PetscArraycpy(a->gcoords_source+A->cmap->rstart*spacedim,coords_source,A->cmap->n*spacedim);
848:     MPIU_Allreduce(MPI_IN_PLACE,a->gcoords_source,A->cmap->N*spacedim,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)A)); /* global source coordinates */
849:   } else a->gcoords_source = a->gcoords_target;
850:   PetscCalloc2(A->cmap->N,&a->work_source,A->rmap->N,&a->work_target);
851:   *B = A;
852:   return 0;
853: }

855: /*MC
856:      MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package.

858:   Use ./configure --download-htool to install PETSc to use Htool.

860:    Options Database Keys:
861: .     -mat_type htool - matrix type to "htool" during a call to MatSetFromOptions()

863:    Level: beginner

865: .seealso: `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()`
866: M*/
867: PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A)
868: {
869:   Mat_Htool *a;

871:   PetscNewLog(A,&a);
872:   A->data = (void*)a;
873:   PetscObjectChangeTypeName((PetscObject)A,MATHTOOL);
874:   PetscMemzero(A->ops,sizeof(struct _MatOps));
875:   A->ops->getdiagonal       = MatGetDiagonal_Htool;
876:   A->ops->getdiagonalblock  = MatGetDiagonalBlock_Htool;
877:   A->ops->mult              = MatMult_Htool;
878:   A->ops->multadd           = MatMultAdd_Htool;
879:   A->ops->multtranspose     = MatMultTranspose_Htool;
880:   if (!PetscDefined(USE_COMPLEX)) A->ops->multhermitiantranspose = MatMultTranspose_Htool;
881:   A->ops->increaseoverlap   = MatIncreaseOverlap_Htool;
882:   A->ops->createsubmatrices = MatCreateSubMatrices_Htool;
883:   A->ops->transpose         = MatTranspose_Htool;
884:   A->ops->destroy           = MatDestroy_Htool;
885:   A->ops->view              = MatView_Htool;
886:   A->ops->setfromoptions    = MatSetFromOptions_Htool;
887:   A->ops->scale             = MatScale_Htool;
888:   A->ops->getrow            = MatGetRow_Htool;
889:   A->ops->restorerow        = MatRestoreRow_Htool;
890:   A->ops->assemblyend       = MatAssemblyEnd_Htool;
891:   a->dim                    = 0;
892:   a->gcoords_target         = NULL;
893:   a->gcoords_source         = NULL;
894:   a->s                      = 1.0;
895:   a->bs[0]                  = 10;
896:   a->bs[1]                  = 1000000;
897:   a->epsilon                = PetscSqrtReal(PETSC_SMALL);
898:   a->eta                    = 10.0;
899:   a->depth[0]               = 0;
900:   a->depth[1]               = 0;
901:   a->compressor             = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA;
902:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_seqdense_C",MatProductSetFromOptions_Htool);
903:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_mpidense_C",MatProductSetFromOptions_Htool);
904:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_seqdense_C",MatConvert_Htool_Dense);
905:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_mpidense_C",MatConvert_Htool_Dense);
906:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetHierarchicalMat_C",MatHtoolGetHierarchicalMat_Htool);
907:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolSetKernel_C",MatHtoolSetKernel_Htool);
908:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationSource_C",MatHtoolGetPermutationSource_Htool);
909:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationTarget_C",MatHtoolGetPermutationTarget_Htool);
910:   PetscObjectComposeFunction((PetscObject)A,"MatHtoolUsePermutation_C",MatHtoolUsePermutation_Htool);
911:   return 0;
912: }