Actual source code: math2opus.cu

  1: #include <h2opusconf.h>
  2: /* skip compilation of this .cu file if H2OPUS is CPU only while PETSc has GPU support */
  4:   #include <h2opus.h>
  5:   #if defined(H2OPUS_USE_MPI)
  6:     #include <h2opus/distributed/distributed_h2opus_handle.h>
  7:     #include <h2opus/distributed/distributed_geometric_construction.h>
  8:     #include <h2opus/distributed/distributed_hgemv.h>
  9:     #include <h2opus/distributed/distributed_horthog.h>
 10:     #include <h2opus/distributed/distributed_hcompress.h>
 11:   #endif
 12:   #include <h2opus/util/boxentrygen.h>
 13: #include <petsc/private/matimpl.h>
 14: #include <petsc/private/vecimpl.h>
 15: #include <petsc/private/deviceimpl.h>
 16: #include <petscsf.h>

 18: /* math2opusutils */
 19: PETSC_INTERN PetscErrorCode MatDenseGetH2OpusStridedSF(Mat, PetscSF, PetscSF *);
 20: PETSC_INTERN PetscErrorCode VecSign(Vec, Vec);
 21: PETSC_INTERN PetscErrorCode VecSetDelta(Vec, PetscInt);
 22: PETSC_INTERN PetscErrorCode MatApproximateNorm_Private(Mat, NormType, PetscInt, PetscReal *);

 24:   #define MatH2OpusGetThrustPointer(v) thrust::raw_pointer_cast((v).data())

 26:   /* Use GPU only if H2OPUS is configured for GPU */
 27:   #if defined(PETSC_HAVE_CUDA) && defined(H2OPUS_USE_GPU)
 28:     #define PETSC_H2OPUS_USE_GPU
 29:   #endif
 30:   #if defined(PETSC_H2OPUS_USE_GPU)
 31:     #define MatH2OpusUpdateIfNeeded(A, B) MatBindToCPU(A, (PetscBool)((A)->boundtocpu || (B)))
 32:   #else
 33:     #define MatH2OpusUpdateIfNeeded(A, B) PETSC_SUCCESS
 34:   #endif

 36: // TODO H2OPUS:
 37: // DistributedHMatrix
 38: //   unsymmetric ?
 39: //   transpose for distributed_hgemv?
 40: //   clearData()
 41: // Unify interface for sequential and parallel?
 42: // Reuse geometric construction (almost possible, only the unsymmetric case is explicitly handled)
 43: //
 44: template <class T>
 45: class PetscPointCloud : public H2OpusDataSet<T> {
 46: private:
 47:   int            dimension;
 48:   size_t         num_points;
 49:   std::vector<T> pts;

 51: public:
 52:   PetscPointCloud(int dim, size_t num_pts, const T coords[])
 53:   {
 54:     dim              = dim > 0 ? dim : 1;
 55:     this->dimension  = dim;
 56:     this->num_points = num_pts;

 58:     pts.resize(num_pts * dim);
 59:     if (coords) {
 60:       for (size_t n = 0; n < num_pts; n++)
 61:         for (int i = 0; i < dim; i++) pts[n * dim + i] = coords[n * dim + i];
 62:     } else {
 63:       PetscReal h = 1.0; //num_pts > 1 ? 1./(num_pts - 1) : 0.0;
 64:       for (size_t n = 0; n < num_pts; n++) {
 65:         pts[n * dim] = n * h;
 66:         for (int i = 1; i < dim; i++) pts[n * dim + i] = 0.0;
 67:       }
 68:     }
 69:   }

 71:   PetscPointCloud(const PetscPointCloud<T> &other)
 72:   {
 73:     size_t N         = other.dimension * other.num_points;
 74:     this->dimension  = other.dimension;
 75:     this->num_points = other.num_points;
 76:     this->pts.resize(N);
 77:     for (size_t i = 0; i < N; i++) this->pts[i] = other.pts[i];
 78:   }

 80:   int getDimension() const { return dimension; }

 82:   size_t getDataSetSize() const { return num_points; }

 84:   T getDataPoint(size_t idx, int dim) const
 85:   {
 86:     assert(dim < dimension && idx < num_points);
 87:     return pts[idx * dimension + dim];
 88:   }

 90:   void Print(std::ostream &out = std::cout)
 91:   {
 92:     out << "Dimension: " << dimension << std::endl;
 93:     out << "NumPoints: " << num_points << std::endl;
 94:     for (size_t n = 0; n < num_points; n++) {
 95:       for (int d = 0; d < dimension; d++) out << pts[n * dimension + d] << " ";
 96:       out << std::endl;
 97:     }
 98:   }
 99: };

101: template <class T>
102: class PetscFunctionGenerator {
103: private:
104:   MatH2OpusKernelFn *k;
105:   int                dim;
106:   void              *ctx;

108: public:
109:   PetscFunctionGenerator(MatH2OpusKernelFn *k, int dim, void *ctx)
110:   {
111:     this->k   = k;
112:     this->dim = dim;
113:     this->ctx = ctx;
114:   }
115:   PetscFunctionGenerator(PetscFunctionGenerator &other)
116:   {
117:     this->k   = other.k;
118:     this->dim = other.dim;
119:     this->ctx = other.ctx;
120:   }
121:   T operator()(PetscReal *pt1, PetscReal *pt2) { return (T)((*this->k)(this->dim, pt1, pt2, this->ctx)); }
122: };

124: #include <../src/mat/impls/h2opus/math2opussampler.hpp>

126:   /* just to not clutter the code */
127:   #if !defined(H2OPUS_USE_GPU)
128: typedef HMatrix HMatrix_GPU;
129:     #if defined(H2OPUS_USE_MPI)
130: typedef DistributedHMatrix DistributedHMatrix_GPU;
131:     #endif
132:   #endif

134: typedef struct {
135:   #if defined(H2OPUS_USE_MPI)
136:   distributedH2OpusHandle_t handle;
137:   #else
138:   h2opusHandle_t handle;
139:   #endif
140:   /* Sequential and parallel matrices are two different classes at the moment */
141:   HMatrix *hmatrix;
142:   #if defined(H2OPUS_USE_MPI)
143:   DistributedHMatrix *dist_hmatrix;
144:   #else
145:   HMatrix *dist_hmatrix; /* just to not clutter the code */
146:   #endif
147:   /* May use permutations */
148:   PetscSF                           sf;
149:   PetscLayout                       h2opus_rmap, h2opus_cmap;
150:   IS                                h2opus_indexmap;
151:   thrust::host_vector<PetscScalar> *xx, *yy;
152:   PetscInt                          xxs, yys;
153:   PetscBool                         multsetup;

155:   /* GPU */
156:   HMatrix_GPU *hmatrix_gpu;
157:   #if defined(H2OPUS_USE_MPI)
158:   DistributedHMatrix_GPU *dist_hmatrix_gpu;
159:   #else
160:   HMatrix_GPU *dist_hmatrix_gpu; /* just to not clutter the code */
161:   #endif
162:   #if defined(PETSC_H2OPUS_USE_GPU)
163:   thrust::device_vector<PetscScalar> *xx_gpu, *yy_gpu;
164:   PetscInt                            xxs_gpu, yys_gpu;
165:   #endif

167:   /* construction from matvecs */
168:   PetscMatrixSampler *sampler;
169:   PetscBool           nativemult;

171:   /* Admissibility */
172:   PetscReal eta;
173:   PetscInt  leafsize;

175:   /* for dof reordering */
176:   PetscPointCloud<PetscReal> *ptcloud;

178:   /* kernel for generating matrix entries */
179:   PetscFunctionGenerator<PetscScalar> *kernel;

181:   /* basis orthogonalized? */
182:   PetscBool orthogonal;

184:   /* customization */
185:   PetscInt  basisord;
186:   PetscInt  max_rank;
187:   PetscInt  bs;
188:   PetscReal rtol;
189:   PetscInt  norm_max_samples;
190:   PetscBool check_construction;
191:   PetscBool hara_verbose;
192:   PetscBool resize;

194:   /* keeps track of MatScale values */
195:   PetscScalar s;
196: } Mat_H2OPUS;

198: static PetscErrorCode MatDestroy_H2OPUS(Mat A)
199: {
200:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

202:   PetscFunctionBegin;
203:   #if defined(H2OPUS_USE_MPI)
204:   h2opusDestroyDistributedHandle(a->handle);
205:   #else
206:   h2opusDestroyHandle(a->handle);
207:   #endif
208:   delete a->dist_hmatrix;
209:   delete a->hmatrix;
210:   PetscCall(PetscSFDestroy(&a->sf));
211:   PetscCall(PetscLayoutDestroy(&a->h2opus_rmap));
212:   PetscCall(PetscLayoutDestroy(&a->h2opus_cmap));
213:   PetscCall(ISDestroy(&a->h2opus_indexmap));
214:   delete a->xx;
215:   delete a->yy;
216:   delete a->hmatrix_gpu;
217:   delete a->dist_hmatrix_gpu;
218:   #if defined(PETSC_H2OPUS_USE_GPU)
219:   delete a->xx_gpu;
220:   delete a->yy_gpu;
221:   #endif
222:   delete a->sampler;
223:   delete a->ptcloud;
224:   delete a->kernel;
225:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", NULL));
226:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", NULL));
227:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", NULL));
228:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", NULL));
229:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
230:   PetscCall(PetscFree(A->data));
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: PetscErrorCode MatH2OpusSetNativeMult(Mat A, PetscBool nm)
235: {
236:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
237:   PetscBool   ish2opus;

239:   PetscFunctionBegin;
242:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
243:   if (ish2opus) {
244:     if (a->h2opus_rmap) { /* need to swap layouts for vector creation */
245:       if ((!a->nativemult && nm) || (a->nativemult && !nm)) {
246:         PetscLayout t;
247:         t              = A->rmap;
248:         A->rmap        = a->h2opus_rmap;
249:         a->h2opus_rmap = t;
250:         t              = A->cmap;
251:         A->cmap        = a->h2opus_cmap;
252:         a->h2opus_cmap = t;
253:       }
254:     }
255:     a->nativemult = nm;
256:   }
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: PetscErrorCode MatH2OpusGetNativeMult(Mat A, PetscBool *nm)
261: {
262:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
263:   PetscBool   ish2opus;

265:   PetscFunctionBegin;
267:   PetscAssertPointer(nm, 2);
268:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
269:   PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
270:   *nm = a->nativemult;
271:   PetscFunctionReturn(PETSC_SUCCESS);
272: }

274: PETSC_EXTERN PetscErrorCode MatNorm_H2OPUS(Mat A, NormType normtype, PetscReal *n)
275: {
276:   PetscBool   ish2opus;
277:   PetscInt    nmax = PETSC_DECIDE;
278:   Mat_H2OPUS *a    = NULL;
279:   PetscBool   mult = PETSC_FALSE;

281:   PetscFunctionBegin;
282:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
283:   if (ish2opus) { /* set userdefine number of samples and fastpath for mult (norms are order independent) */
284:     a = (Mat_H2OPUS *)A->data;

286:     nmax = a->norm_max_samples;
287:     mult = a->nativemult;
288:     PetscCall(MatH2OpusSetNativeMult(A, PETSC_TRUE));
289:   } else {
290:     PetscCall(PetscOptionsGetInt(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_approximate_norm_samples", &nmax, NULL));
291:   }
292:   PetscCall(MatApproximateNorm_Private(A, normtype, nmax, n));
293:   if (a) PetscCall(MatH2OpusSetNativeMult(A, mult));
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: static PetscErrorCode MatH2OpusResizeBuffers_Private(Mat A, PetscInt xN, PetscInt yN)
298: {
299:   Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
300:   PetscInt    n;
301:   PetscBool   boundtocpu = PETSC_TRUE;

303:   PetscFunctionBegin;
304:   #if defined(PETSC_H2OPUS_USE_GPU)
305:   boundtocpu = A->boundtocpu;
306:   #endif
307:   PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
308:   if (boundtocpu) {
309:     if (h2opus->xxs < xN) {
310:       h2opus->xx->resize(n * xN);
311:       h2opus->xxs = xN;
312:     }
313:     if (h2opus->yys < yN) {
314:       h2opus->yy->resize(n * yN);
315:       h2opus->yys = yN;
316:     }
317:   }
318:   #if defined(PETSC_H2OPUS_USE_GPU)
319:   if (!boundtocpu) {
320:     if (h2opus->xxs_gpu < xN) {
321:       h2opus->xx_gpu->resize(n * xN);
322:       h2opus->xxs_gpu = xN;
323:     }
324:     if (h2opus->yys_gpu < yN) {
325:       h2opus->yy_gpu->resize(n * yN);
326:       h2opus->yys_gpu = yN;
327:     }
328:   }
329:   #endif
330:   PetscFunctionReturn(PETSC_SUCCESS);
331: }

333: static PetscErrorCode MatMultNKernel_H2OPUS(Mat A, PetscBool transA, Mat B, Mat C)
334: {
335:   Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
336:   #if defined(H2OPUS_USE_MPI)
337:   h2opusHandle_t handle = h2opus->handle->handle;
338:   #else
339:   h2opusHandle_t handle = h2opus->handle;
340:   #endif
341:   PetscBool    boundtocpu = PETSC_TRUE;
342:   PetscScalar *xx, *yy, *uxx, *uyy;
343:   PetscInt     blda, clda;
344:   PetscMPIInt  size;
345:   PetscSF      bsf, csf;
346:   PetscBool    usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);

348:   PetscFunctionBegin;
349:   HLibProfile::clear();
350:   #if defined(PETSC_H2OPUS_USE_GPU)
351:   boundtocpu = A->boundtocpu;
352:   #endif
353:   PetscCall(MatDenseGetLDA(B, &blda));
354:   PetscCall(MatDenseGetLDA(C, &clda));
355:   if (usesf) {
356:     PetscInt n;

358:     PetscCall(MatDenseGetH2OpusStridedSF(B, h2opus->sf, &bsf));
359:     PetscCall(MatDenseGetH2OpusStridedSF(C, h2opus->sf, &csf));

361:     PetscCall(MatH2OpusResizeBuffers_Private(A, B->cmap->N, C->cmap->N));
362:     PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
363:     blda = n;
364:     clda = n;
365:   }
366:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
367:   if (boundtocpu) {
368:     PetscCall(MatDenseGetArrayRead(B, (const PetscScalar **)&xx));
369:     PetscCall(MatDenseGetArrayWrite(C, &yy));
370:     if (usesf) {
371:       uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
372:       uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
373:       PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
374:       PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
375:     } else {
376:       uxx = xx;
377:       uyy = yy;
378:     }
379:     if (size > 1) {
380:       PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
381:       PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
382:   #if defined(H2OPUS_USE_MPI)
383:       distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
384:   #endif
385:     } else {
386:       PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
387:       hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
388:     }
389:     PetscCall(MatDenseRestoreArrayRead(B, (const PetscScalar **)&xx));
390:     if (usesf) {
391:       PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
392:       PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
393:     }
394:     PetscCall(MatDenseRestoreArrayWrite(C, &yy));
395:   #if defined(PETSC_H2OPUS_USE_GPU)
396:   } else {
397:     PetscBool ciscuda, biscuda;

399:     /* If not of type seqdensecuda, convert on the fly (i.e. allocate GPU memory) */
400:     PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &biscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
401:     if (!biscuda) PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
402:     PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &ciscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
403:     if (!ciscuda) {
404:       C->assembled = PETSC_TRUE;
405:       PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
406:     }
407:     PetscCall(MatDenseCUDAGetArrayRead(B, (const PetscScalar **)&xx));
408:     PetscCall(MatDenseCUDAGetArrayWrite(C, &yy));
409:     if (usesf) {
410:       uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
411:       uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
412:       PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
413:       PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
414:     } else {
415:       uxx = xx;
416:       uyy = yy;
417:     }
418:     PetscCall(PetscLogGpuTimeBegin());
419:     if (size > 1) {
420:       PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
421:       PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
422:     #if defined(H2OPUS_USE_MPI)
423:       distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
424:     #endif
425:     } else {
426:       PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
427:       hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
428:     }
429:     PetscCall(PetscLogGpuTimeEnd());
430:     PetscCall(MatDenseCUDARestoreArrayRead(B, (const PetscScalar **)&xx));
431:     if (usesf) {
432:       PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
433:       PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
434:     }
435:     PetscCall(MatDenseCUDARestoreArrayWrite(C, &yy));
436:     if (!biscuda) PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
437:     if (!ciscuda) PetscCall(MatConvert(C, MATDENSE, MAT_INPLACE_MATRIX, &C));
438:   #endif
439:   }
440:   { /* log flops */
441:     double gops, time, perf, dev;
442:     HLibProfile::getHgemvPerf(gops, time, perf, dev);
443:   #if defined(PETSC_H2OPUS_USE_GPU)
444:     if (boundtocpu) {
445:       PetscCall(PetscLogFlops(1e9 * gops));
446:     } else {
447:       PetscCall(PetscLogGpuFlops(1e9 * gops));
448:     }
449:   #else
450:     PetscCall(PetscLogFlops(1e9 * gops));
451:   #endif
452:   }
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: static PetscErrorCode MatProductNumeric_H2OPUS(Mat C)
457: {
458:   Mat_Product *product = C->product;

460:   PetscFunctionBegin;
461:   MatCheckProduct(C, 1);
462:   switch (product->type) {
463:   case MATPRODUCT_AB:
464:     PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_FALSE, product->B, C));
465:     break;
466:   case MATPRODUCT_AtB:
467:     PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_TRUE, product->B, C));
468:     break;
469:   default:
470:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
471:   }
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: static PetscErrorCode MatProductSymbolic_H2OPUS(Mat C)
476: {
477:   Mat_Product *product = C->product;
478:   PetscBool    cisdense;
479:   Mat          A, B;

481:   PetscFunctionBegin;
482:   MatCheckProduct(C, 1);
483:   A = product->A;
484:   B = product->B;
485:   switch (product->type) {
486:   case MATPRODUCT_AB:
487:     PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
488:     PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
489:     PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
490:     if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
491:     PetscCall(MatSetUp(C));
492:     break;
493:   case MATPRODUCT_AtB:
494:     PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
495:     PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
496:     PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
497:     if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
498:     PetscCall(MatSetUp(C));
499:     break;
500:   default:
501:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
502:   }
503:   C->ops->productsymbolic = NULL;
504:   C->ops->productnumeric  = MatProductNumeric_H2OPUS;
505:   PetscFunctionReturn(PETSC_SUCCESS);
506: }

508: static PetscErrorCode MatProductSetFromOptions_H2OPUS(Mat C)
509: {
510:   PetscFunctionBegin;
511:   MatCheckProduct(C, 1);
512:   if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_H2OPUS;
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: static PetscErrorCode MatMultKernel_H2OPUS(Mat A, Vec x, PetscScalar sy, Vec y, PetscBool trans)
517: {
518:   Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
519:   #if defined(H2OPUS_USE_MPI)
520:   h2opusHandle_t handle = h2opus->handle->handle;
521:   #else
522:   h2opusHandle_t handle = h2opus->handle;
523:   #endif
524:   PetscBool    boundtocpu = PETSC_TRUE;
525:   PetscInt     n;
526:   PetscScalar *xx, *yy, *uxx, *uyy;
527:   PetscMPIInt  size;
528:   PetscBool    usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);

530:   PetscFunctionBegin;
531:   HLibProfile::clear();
532:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
533:   #if defined(PETSC_H2OPUS_USE_GPU)
534:   boundtocpu = A->boundtocpu;
535:   #endif
536:   if (usesf) {
537:     PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
538:   } else n = A->rmap->n;
539:   if (boundtocpu) {
540:     PetscCall(VecGetArrayRead(x, (const PetscScalar **)&xx));
541:     if (sy == 0.0) {
542:       PetscCall(VecGetArrayWrite(y, &yy));
543:     } else {
544:       PetscCall(VecGetArray(y, &yy));
545:     }
546:     if (usesf) {
547:       uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
548:       uyy = MatH2OpusGetThrustPointer(*h2opus->yy);

550:       PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
551:       PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
552:       if (sy != 0.0) {
553:         PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
554:         PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
555:       }
556:     } else {
557:       uxx = xx;
558:       uyy = yy;
559:     }
560:     if (size > 1) {
561:       PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
562:       PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
563:   #if defined(H2OPUS_USE_MPI)
564:       distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, n, sy, uyy, n, 1, h2opus->handle);
565:   #endif
566:     } else {
567:       PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
568:       hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, n, sy, uyy, n, 1, handle);
569:     }
570:     PetscCall(VecRestoreArrayRead(x, (const PetscScalar **)&xx));
571:     if (usesf) {
572:       PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
573:       PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
574:     }
575:     if (sy == 0.0) {
576:       PetscCall(VecRestoreArrayWrite(y, &yy));
577:     } else {
578:       PetscCall(VecRestoreArray(y, &yy));
579:     }
580:   #if defined(PETSC_H2OPUS_USE_GPU)
581:   } else {
582:     PetscCall(VecCUDAGetArrayRead(x, (const PetscScalar **)&xx));
583:     if (sy == 0.0) {
584:       PetscCall(VecCUDAGetArrayWrite(y, &yy));
585:     } else {
586:       PetscCall(VecCUDAGetArray(y, &yy));
587:     }
588:     if (usesf) {
589:       uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
590:       uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);

592:       PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
593:       PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
594:       if (sy != 0.0) {
595:         PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
596:         PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
597:       }
598:     } else {
599:       uxx = xx;
600:       uyy = yy;
601:     }
602:     PetscCall(PetscLogGpuTimeBegin());
603:     if (size > 1) {
604:       PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
605:       PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
606:     #if defined(H2OPUS_USE_MPI)
607:       distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, n, sy, uyy, n, 1, h2opus->handle);
608:     #endif
609:     } else {
610:       PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
611:       hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, n, sy, uyy, n, 1, handle);
612:     }
613:     PetscCall(PetscLogGpuTimeEnd());
614:     PetscCall(VecCUDARestoreArrayRead(x, (const PetscScalar **)&xx));
615:     if (usesf) {
616:       PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
617:       PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
618:     }
619:     if (sy == 0.0) {
620:       PetscCall(VecCUDARestoreArrayWrite(y, &yy));
621:     } else {
622:       PetscCall(VecCUDARestoreArray(y, &yy));
623:     }
624:   #endif
625:   }
626:   { /* log flops */
627:     double gops, time, perf, dev;
628:     HLibProfile::getHgemvPerf(gops, time, perf, dev);
629:   #if defined(PETSC_H2OPUS_USE_GPU)
630:     if (boundtocpu) {
631:       PetscCall(PetscLogFlops(1e9 * gops));
632:     } else {
633:       PetscCall(PetscLogGpuFlops(1e9 * gops));
634:     }
635:   #else
636:     PetscCall(PetscLogFlops(1e9 * gops));
637:   #endif
638:   }
639:   PetscFunctionReturn(PETSC_SUCCESS);
640: }

642: static PetscErrorCode MatMultTranspose_H2OPUS(Mat A, Vec x, Vec y)
643: {
644:   PetscBool xiscuda, yiscuda;

646:   PetscFunctionBegin;
647:   PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
648:   PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
649:   PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
650:   PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_TRUE));
651:   PetscFunctionReturn(PETSC_SUCCESS);
652: }

654: static PetscErrorCode MatMult_H2OPUS(Mat A, Vec x, Vec y)
655: {
656:   PetscBool xiscuda, yiscuda;

658:   PetscFunctionBegin;
659:   PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
660:   PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
661:   PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
662:   PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_FALSE));
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: static PetscErrorCode MatMultTransposeAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
667: {
668:   PetscBool xiscuda, ziscuda;

670:   PetscFunctionBegin;
671:   PetscCall(VecCopy(y, z));
672:   PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
673:   PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
674:   PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
675:   PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_TRUE));
676:   PetscFunctionReturn(PETSC_SUCCESS);
677: }

679: static PetscErrorCode MatMultAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
680: {
681:   PetscBool xiscuda, ziscuda;

683:   PetscFunctionBegin;
684:   PetscCall(VecCopy(y, z));
685:   PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
686:   PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
687:   PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
688:   PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_FALSE));
689:   PetscFunctionReturn(PETSC_SUCCESS);
690: }

692: static PetscErrorCode MatScale_H2OPUS(Mat A, PetscScalar s)
693: {
694:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

696:   PetscFunctionBegin;
697:   a->s *= s;
698:   PetscFunctionReturn(PETSC_SUCCESS);
699: }

701: static PetscErrorCode MatSetFromOptions_H2OPUS(Mat A, PetscOptionItems *PetscOptionsObject)
702: {
703:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

705:   PetscFunctionBegin;
706:   PetscOptionsHeadBegin(PetscOptionsObject, "H2OPUS options");
707:   PetscCall(PetscOptionsInt("-mat_h2opus_leafsize", "Leaf size of cluster tree", NULL, a->leafsize, &a->leafsize, NULL));
708:   PetscCall(PetscOptionsReal("-mat_h2opus_eta", "Admissibility condition tolerance", NULL, a->eta, &a->eta, NULL));
709:   PetscCall(PetscOptionsInt("-mat_h2opus_order", "Basis order for off-diagonal sampling when constructed from kernel", NULL, a->basisord, &a->basisord, NULL));
710:   PetscCall(PetscOptionsInt("-mat_h2opus_maxrank", "Maximum rank when constructed from matvecs", NULL, a->max_rank, &a->max_rank, NULL));
711:   PetscCall(PetscOptionsInt("-mat_h2opus_samples", "Maximum number of samples to be taken concurrently when constructing from matvecs", NULL, a->bs, &a->bs, NULL));
712:   PetscCall(PetscOptionsInt("-mat_h2opus_normsamples", "Maximum number of samples to be when estimating norms", NULL, a->norm_max_samples, &a->norm_max_samples, NULL));
713:   PetscCall(PetscOptionsReal("-mat_h2opus_rtol", "Relative tolerance for construction from sampling", NULL, a->rtol, &a->rtol, NULL));
714:   PetscCall(PetscOptionsBool("-mat_h2opus_check", "Check error when constructing from sampling during MatAssemblyEnd()", NULL, a->check_construction, &a->check_construction, NULL));
715:   PetscCall(PetscOptionsBool("-mat_h2opus_hara_verbose", "Verbose output from hara construction", NULL, a->hara_verbose, &a->hara_verbose, NULL));
716:   PetscCall(PetscOptionsBool("-mat_h2opus_resize", "Resize after compression", NULL, a->resize, &a->resize, NULL));
717:   PetscOptionsHeadEnd();
718:   PetscFunctionReturn(PETSC_SUCCESS);
719: }

721: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat, PetscInt, const PetscReal[], PetscBool, MatH2OpusKernelFn *, void *);

723: static PetscErrorCode MatH2OpusInferCoordinates_Private(Mat A)
724: {
725:   Mat_H2OPUS        *a = (Mat_H2OPUS *)A->data;
726:   Vec                c;
727:   PetscInt           spacedim;
728:   const PetscScalar *coords;

730:   PetscFunctionBegin;
731:   if (a->ptcloud) PetscFunctionReturn(PETSC_SUCCESS);
732:   PetscCall(PetscObjectQuery((PetscObject)A, "__math2opus_coords", (PetscObject *)&c));
733:   if (!c && a->sampler) {
734:     Mat S = a->sampler->GetSamplingMat();

736:     PetscCall(PetscObjectQuery((PetscObject)S, "__math2opus_coords", (PetscObject *)&c));
737:   }
738:   if (!c) {
739:     PetscCall(MatH2OpusSetCoords_H2OPUS(A, -1, NULL, PETSC_FALSE, NULL, NULL));
740:   } else {
741:     PetscCall(VecGetArrayRead(c, &coords));
742:     PetscCall(VecGetBlockSize(c, &spacedim));
743:     PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, PETSC_FALSE, NULL, NULL));
744:     PetscCall(VecRestoreArrayRead(c, &coords));
745:   }
746:   PetscFunctionReturn(PETSC_SUCCESS);
747: }

749: static PetscErrorCode MatSetUpMultiply_H2OPUS(Mat A)
750: {
751:   MPI_Comm      comm;
752:   PetscMPIInt   size;
753:   Mat_H2OPUS   *a = (Mat_H2OPUS *)A->data;
754:   PetscInt      n = 0, *idx = NULL;
755:   int          *iidx = NULL;
756:   PetscCopyMode own;
757:   PetscBool     rid;

759:   PetscFunctionBegin;
760:   if (a->multsetup) PetscFunctionReturn(PETSC_SUCCESS);
761:   if (a->sf) { /* MatDuplicate_H2OPUS takes reference to the SF */
762:     PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
763:   #if defined(PETSC_H2OPUS_USE_GPU)
764:     a->xx_gpu  = new thrust::device_vector<PetscScalar>(n);
765:     a->yy_gpu  = new thrust::device_vector<PetscScalar>(n);
766:     a->xxs_gpu = 1;
767:     a->yys_gpu = 1;
768:   #endif
769:     a->xx  = new thrust::host_vector<PetscScalar>(n);
770:     a->yy  = new thrust::host_vector<PetscScalar>(n);
771:     a->xxs = 1;
772:     a->yys = 1;
773:   } else {
774:     IS is;
775:     PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
776:     PetscCallMPI(MPI_Comm_size(comm, &size));
777:     if (!a->h2opus_indexmap) {
778:       if (size > 1) {
779:         PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
780:   #if defined(H2OPUS_USE_MPI)
781:         iidx = MatH2OpusGetThrustPointer(a->dist_hmatrix->basis_tree.basis_branch.index_map);
782:         n    = a->dist_hmatrix->basis_tree.basis_branch.index_map.size();
783:   #endif
784:       } else {
785:         iidx = MatH2OpusGetThrustPointer(a->hmatrix->u_basis_tree.index_map);
786:         n    = a->hmatrix->u_basis_tree.index_map.size();
787:       }

789:       if (PetscDefined(USE_64BIT_INDICES)) {
790:         PetscInt i;

792:         own = PETSC_OWN_POINTER;
793:         PetscCall(PetscMalloc1(n, &idx));
794:         for (i = 0; i < n; i++) idx[i] = iidx[i];
795:       } else {
796:         own = PETSC_COPY_VALUES;
797:         idx = (PetscInt *)iidx;
798:       }
799:       PetscCall(ISCreateGeneral(comm, n, idx, own, &is));
800:       PetscCall(ISSetPermutation(is));
801:       PetscCall(ISViewFromOptions(is, (PetscObject)A, "-mat_h2opus_indexmap_view"));
802:       a->h2opus_indexmap = is;
803:     }
804:     PetscCall(ISGetLocalSize(a->h2opus_indexmap, &n));
805:     PetscCall(ISGetIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
806:     rid = (PetscBool)(n == A->rmap->n);
807:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &rid, 1, MPIU_BOOL, MPI_LAND, comm));
808:     if (rid) PetscCall(ISIdentity(a->h2opus_indexmap, &rid));
809:     if (!rid) {
810:       if (size > 1) { /* Parallel distribution may be different, save it here for fast path in MatMult (see MatH2OpusSetNativeMult) */
811:         PetscCall(PetscLayoutCreate(comm, &a->h2opus_rmap));
812:         PetscCall(PetscLayoutSetLocalSize(a->h2opus_rmap, n));
813:         PetscCall(PetscLayoutSetUp(a->h2opus_rmap));
814:         PetscCall(PetscLayoutReference(a->h2opus_rmap, &a->h2opus_cmap));
815:       }
816:       PetscCall(PetscSFCreate(comm, &a->sf));
817:       PetscCall(PetscSFSetGraphLayout(a->sf, A->rmap, n, NULL, PETSC_OWN_POINTER, idx));
818:       PetscCall(PetscSFSetUp(a->sf));
819:       PetscCall(PetscSFViewFromOptions(a->sf, (PetscObject)A, "-mat_h2opus_sf_view"));
820:   #if defined(PETSC_H2OPUS_USE_GPU)
821:       a->xx_gpu  = new thrust::device_vector<PetscScalar>(n);
822:       a->yy_gpu  = new thrust::device_vector<PetscScalar>(n);
823:       a->xxs_gpu = 1;
824:       a->yys_gpu = 1;
825:   #endif
826:       a->xx  = new thrust::host_vector<PetscScalar>(n);
827:       a->yy  = new thrust::host_vector<PetscScalar>(n);
828:       a->xxs = 1;
829:       a->yys = 1;
830:     }
831:     PetscCall(ISRestoreIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
832:   }
833:   a->multsetup = PETSC_TRUE;
834:   PetscFunctionReturn(PETSC_SUCCESS);
835: }

837: static PetscErrorCode MatAssemblyEnd_H2OPUS(Mat A, MatAssemblyType assemblytype)
838: {
839:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
840:   #if defined(H2OPUS_USE_MPI)
841:   h2opusHandle_t handle = a->handle->handle;
842:   #else
843:   h2opusHandle_t handle = a->handle;
844:   #endif
845:   PetscBool   kernel       = PETSC_FALSE;
846:   PetscBool   boundtocpu   = PETSC_TRUE;
847:   PetscBool   samplingdone = PETSC_FALSE;
848:   MPI_Comm    comm;
849:   PetscMPIInt size;

851:   PetscFunctionBegin;
852:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
853:   PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
854:   PetscCheck(A->rmap->N == A->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");

856:   /* XXX */
857:   a->leafsize = PetscMin(a->leafsize, PetscMin(A->rmap->N, A->cmap->N));

859:   PetscCallMPI(MPI_Comm_size(comm, &size));
860:   /* TODO REUSABILITY of geometric construction */
861:   delete a->hmatrix;
862:   delete a->dist_hmatrix;
863:   #if defined(PETSC_H2OPUS_USE_GPU)
864:   delete a->hmatrix_gpu;
865:   delete a->dist_hmatrix_gpu;
866:   #endif
867:   a->orthogonal = PETSC_FALSE;

869:   /* TODO: other? */
870:   H2OpusBoxCenterAdmissibility adm(a->eta);

872:   PetscCall(PetscLogEventBegin(MAT_H2Opus_Build, A, 0, 0, 0));
873:   if (size > 1) {
874:   #if defined(H2OPUS_USE_MPI)
875:     a->dist_hmatrix = new DistributedHMatrix(A->rmap->n /* ,A->symmetric */);
876:   #else
877:     a->dist_hmatrix = NULL;
878:   #endif
879:   } else a->hmatrix = new HMatrix(A->rmap->n, A->symmetric == PETSC_BOOL3_TRUE);
880:   PetscCall(MatH2OpusInferCoordinates_Private(A));
881:   PetscCheck(a->ptcloud, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing pointcloud");
882:   if (a->kernel) {
883:     BoxEntryGen<PetscScalar, H2OPUS_HWTYPE_CPU, PetscFunctionGenerator<PetscScalar>> entry_gen(*a->kernel);
884:     if (size > 1) {
885:       PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
886:   #if defined(H2OPUS_USE_MPI)
887:       buildDistributedHMatrix(*a->dist_hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord, a->handle);
888:   #endif
889:     } else {
890:       buildHMatrix(*a->hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord);
891:     }
892:     kernel = PETSC_TRUE;
893:   } else {
894:     PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Construction from sampling not supported in parallel");
895:     buildHMatrixStructure(*a->hmatrix, a->ptcloud, a->leafsize, adm);
896:   }
897:   PetscCall(MatSetUpMultiply_H2OPUS(A));

899:   #if defined(PETSC_H2OPUS_USE_GPU)
900:   boundtocpu = A->boundtocpu;
901:   if (!boundtocpu) {
902:     if (size > 1) {
903:       PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
904:     #if defined(H2OPUS_USE_MPI)
905:       a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
906:     #endif
907:     } else {
908:       a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
909:     }
910:   }
911:   #endif
912:   if (size == 1) {
913:     if (!kernel && a->sampler && a->sampler->GetSamplingMat()) {
914:       PetscReal Anorm;
915:       bool      verbose;

917:       PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_hara_verbose", &a->hara_verbose, NULL));
918:       verbose = a->hara_verbose;
919:       PetscCall(MatApproximateNorm_Private(a->sampler->GetSamplingMat(), NORM_2, a->norm_max_samples, &Anorm));
920:       if (a->hara_verbose) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Sampling uses max rank %d, tol %g (%g*%g), %s samples %d\n", a->max_rank, a->rtol * Anorm, a->rtol, Anorm, boundtocpu ? "CPU" : "GPU", a->bs));
921:       if (a->sf && !a->nativemult) a->sampler->SetIndexMap(a->hmatrix->u_basis_tree.index_map.size(), a->hmatrix->u_basis_tree.index_map.data());
922:       a->sampler->SetStream(handle->getMainStream());
923:       if (boundtocpu) {
924:         a->sampler->SetGPUSampling(false);
925:         hara(a->sampler, *a->hmatrix, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
926:   #if defined(PETSC_H2OPUS_USE_GPU)
927:       } else {
928:         a->sampler->SetGPUSampling(true);
929:         hara(a->sampler, *a->hmatrix_gpu, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
930:   #endif
931:       }
932:       samplingdone = PETSC_TRUE;
933:     }
934:   }
935:   #if defined(PETSC_H2OPUS_USE_GPU)
936:   if (!boundtocpu) {
937:     delete a->hmatrix;
938:     delete a->dist_hmatrix;
939:     a->hmatrix      = NULL;
940:     a->dist_hmatrix = NULL;
941:   }
942:   A->offloadmask = boundtocpu ? PETSC_OFFLOAD_CPU : PETSC_OFFLOAD_GPU;
943:   #endif
944:   PetscCall(PetscLogEventEnd(MAT_H2Opus_Build, A, 0, 0, 0));

946:   if (!a->s) a->s = 1.0;
947:   A->assembled = PETSC_TRUE;

949:   if (samplingdone) {
950:     PetscBool check  = a->check_construction;
951:     PetscBool checke = PETSC_FALSE;

953:     PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check", &check, NULL));
954:     PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check_explicit", &checke, NULL));
955:     if (check) {
956:       Mat       E, Ae;
957:       PetscReal n1, ni, n2;
958:       PetscReal n1A, niA, n2A;
959:       void (*normfunc)(void);

961:       Ae = a->sampler->GetSamplingMat();
962:       PetscCall(MatConvert(A, MATSHELL, MAT_INITIAL_MATRIX, &E));
963:       PetscCall(MatShellSetOperation(E, MATOP_NORM, (void (*)(void))MatNorm_H2OPUS));
964:       PetscCall(MatAXPY(E, -1.0, Ae, DIFFERENT_NONZERO_PATTERN));
965:       PetscCall(MatNorm(E, NORM_1, &n1));
966:       PetscCall(MatNorm(E, NORM_INFINITY, &ni));
967:       PetscCall(MatNorm(E, NORM_2, &n2));
968:       if (checke) {
969:         Mat eA, eE, eAe;

971:         PetscCall(MatComputeOperator(A, MATAIJ, &eA));
972:         PetscCall(MatComputeOperator(E, MATAIJ, &eE));
973:         PetscCall(MatComputeOperator(Ae, MATAIJ, &eAe));
974:         PetscCall(MatFilter(eA, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
975:         PetscCall(MatFilter(eE, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
976:         PetscCall(MatFilter(eAe, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
977:         PetscCall(PetscObjectSetName((PetscObject)eA, "H2Mat"));
978:         PetscCall(MatView(eA, NULL));
979:         PetscCall(PetscObjectSetName((PetscObject)eAe, "S"));
980:         PetscCall(MatView(eAe, NULL));
981:         PetscCall(PetscObjectSetName((PetscObject)eE, "H2Mat - S"));
982:         PetscCall(MatView(eE, NULL));
983:         PetscCall(MatDestroy(&eA));
984:         PetscCall(MatDestroy(&eE));
985:         PetscCall(MatDestroy(&eAe));
986:       }

988:       PetscCall(MatGetOperation(Ae, MATOP_NORM, &normfunc));
989:       PetscCall(MatSetOperation(Ae, MATOP_NORM, (void (*)(void))MatNorm_H2OPUS));
990:       PetscCall(MatNorm(Ae, NORM_1, &n1A));
991:       PetscCall(MatNorm(Ae, NORM_INFINITY, &niA));
992:       PetscCall(MatNorm(Ae, NORM_2, &n2A));
993:       n1A = PetscMax(n1A, PETSC_SMALL);
994:       n2A = PetscMax(n2A, PETSC_SMALL);
995:       niA = PetscMax(niA, PETSC_SMALL);
996:       PetscCall(MatSetOperation(Ae, MATOP_NORM, normfunc));
997:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)A), "MATH2OPUS construction errors: NORM_1 %g, NORM_INFINITY %g, NORM_2 %g (%g %g %g)\n", (double)n1, (double)ni, (double)n2, (double)(n1 / n1A), (double)(ni / niA), (double)(n2 / n2A)));
998:       PetscCall(MatDestroy(&E));
999:     }
1000:     a->sampler->SetSamplingMat(NULL);
1001:   }
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

1005: static PetscErrorCode MatZeroEntries_H2OPUS(Mat A)
1006: {
1007:   PetscMPIInt size;
1008:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

1010:   PetscFunctionBegin;
1011:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1012:   PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet supported");
1013:   a->hmatrix->clearData();
1014:   #if defined(PETSC_H2OPUS_USE_GPU)
1015:   if (a->hmatrix_gpu) a->hmatrix_gpu->clearData();
1016:   #endif
1017:   PetscFunctionReturn(PETSC_SUCCESS);
1018: }

1020: static PetscErrorCode MatDuplicate_H2OPUS(Mat B, MatDuplicateOption op, Mat *nA)
1021: {
1022:   Mat         A;
1023:   Mat_H2OPUS *a, *b = (Mat_H2OPUS *)B->data;
1024:   #if defined(PETSC_H2OPUS_USE_GPU)
1025:   PetscBool iscpu = PETSC_FALSE;
1026:   #else
1027:   PetscBool iscpu = PETSC_TRUE;
1028:   #endif
1029:   MPI_Comm comm;

1031:   PetscFunctionBegin;
1032:   PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1033:   PetscCall(MatCreate(comm, &A));
1034:   PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1035:   PetscCall(MatSetType(A, MATH2OPUS));
1036:   PetscCall(MatPropagateSymmetryOptions(B, A));
1037:   a = (Mat_H2OPUS *)A->data;

1039:   a->eta              = b->eta;
1040:   a->leafsize         = b->leafsize;
1041:   a->basisord         = b->basisord;
1042:   a->max_rank         = b->max_rank;
1043:   a->bs               = b->bs;
1044:   a->rtol             = b->rtol;
1045:   a->norm_max_samples = b->norm_max_samples;
1046:   if (op == MAT_COPY_VALUES) a->s = b->s;

1048:   a->ptcloud = new PetscPointCloud<PetscReal>(*b->ptcloud);
1049:   if (op == MAT_COPY_VALUES && b->kernel) a->kernel = new PetscFunctionGenerator<PetscScalar>(*b->kernel);

1051:   #if defined(H2OPUS_USE_MPI)
1052:   if (b->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*b->dist_hmatrix);
1053:     #if defined(PETSC_H2OPUS_USE_GPU)
1054:   if (b->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*b->dist_hmatrix_gpu);
1055:     #endif
1056:   #endif
1057:   if (b->hmatrix) {
1058:     a->hmatrix = new HMatrix(*b->hmatrix);
1059:     if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix->clearData();
1060:   }
1061:   #if defined(PETSC_H2OPUS_USE_GPU)
1062:   if (b->hmatrix_gpu) {
1063:     a->hmatrix_gpu = new HMatrix_GPU(*b->hmatrix_gpu);
1064:     if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix_gpu->clearData();
1065:   }
1066:   #endif
1067:   if (b->sf) {
1068:     PetscCall(PetscObjectReference((PetscObject)b->sf));
1069:     a->sf = b->sf;
1070:   }
1071:   if (b->h2opus_indexmap) {
1072:     PetscCall(PetscObjectReference((PetscObject)b->h2opus_indexmap));
1073:     a->h2opus_indexmap = b->h2opus_indexmap;
1074:   }

1076:   PetscCall(MatSetUp(A));
1077:   PetscCall(MatSetUpMultiply_H2OPUS(A));
1078:   if (op == MAT_COPY_VALUES) {
1079:     A->assembled  = PETSC_TRUE;
1080:     a->orthogonal = b->orthogonal;
1081:   #if defined(PETSC_H2OPUS_USE_GPU)
1082:     A->offloadmask = B->offloadmask;
1083:   #endif
1084:   }
1085:   #if defined(PETSC_H2OPUS_USE_GPU)
1086:   iscpu = B->boundtocpu;
1087:   #endif
1088:   PetscCall(MatBindToCPU(A, iscpu));

1090:   *nA = A;
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

1094: static PetscErrorCode MatView_H2OPUS(Mat A, PetscViewer view)
1095: {
1096:   Mat_H2OPUS       *h2opus = (Mat_H2OPUS *)A->data;
1097:   PetscBool         isascii, vieweps;
1098:   PetscMPIInt       size;
1099:   PetscViewerFormat format;

1101:   PetscFunctionBegin;
1102:   PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
1103:   PetscCall(PetscViewerGetFormat(view, &format));
1104:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1105:   if (isascii) {
1106:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1107:       if (size == 1) {
1108:         FILE *fp;
1109:         PetscCall(PetscViewerASCIIGetPointer(view, &fp));
1110:         dumpHMatrix(*h2opus->hmatrix, 6, fp);
1111:       }
1112:     } else {
1113:       PetscCall(PetscViewerASCIIPrintf(view, "  H-Matrix constructed from %s\n", h2opus->kernel ? "Kernel" : "Mat"));
1114:       PetscCall(PetscViewerASCIIPrintf(view, "  PointCloud dim %" PetscInt_FMT "\n", h2opus->ptcloud ? h2opus->ptcloud->getDimension() : 0));
1115:       PetscCall(PetscViewerASCIIPrintf(view, "  Admissibility parameters: leaf size %" PetscInt_FMT ", eta %g\n", h2opus->leafsize, (double)h2opus->eta));
1116:       if (!h2opus->kernel) {
1117:         PetscCall(PetscViewerASCIIPrintf(view, "  Sampling parameters: max_rank %" PetscInt_FMT ", samples %" PetscInt_FMT ", tolerance %g\n", h2opus->max_rank, h2opus->bs, (double)h2opus->rtol));
1118:       } else {
1119:         PetscCall(PetscViewerASCIIPrintf(view, "  Off-diagonal blocks approximation order %" PetscInt_FMT "\n", h2opus->basisord));
1120:       }
1121:       PetscCall(PetscViewerASCIIPrintf(view, "  Number of samples for norms %" PetscInt_FMT "\n", h2opus->norm_max_samples));
1122:       if (size == 1) {
1123:         double dense_mem_cpu = h2opus->hmatrix ? h2opus->hmatrix->getDenseMemoryUsage() : 0;
1124:         double low_rank_cpu  = h2opus->hmatrix ? h2opus->hmatrix->getLowRankMemoryUsage() : 0;
1125:   #if defined(PETSC_HAVE_CUDA)
1126:         double dense_mem_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getDenseMemoryUsage() : 0;
1127:         double low_rank_gpu  = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getLowRankMemoryUsage() : 0;
1128:   #endif
1129:         PetscCall(PetscViewerASCIIPrintf(view, "  Memory consumption GB (CPU): %g (dense) %g (low rank) %g (total)\n", dense_mem_cpu, low_rank_cpu, low_rank_cpu + dense_mem_cpu));
1130:   #if defined(PETSC_HAVE_CUDA)
1131:         PetscCall(PetscViewerASCIIPrintf(view, "  Memory consumption GB (GPU): %g (dense) %g (low rank) %g (total)\n", dense_mem_gpu, low_rank_gpu, low_rank_gpu + dense_mem_gpu));
1132:   #endif
1133:       } else {
1134:   #if defined(PETSC_HAVE_CUDA)
1135:         double      matrix_mem[4] = {0., 0., 0., 0.};
1136:         PetscMPIInt rsize         = 4;
1137:   #else
1138:         double      matrix_mem[2] = {0., 0.};
1139:         PetscMPIInt rsize         = 2;
1140:   #endif
1141:   #if defined(H2OPUS_USE_MPI)
1142:         matrix_mem[0] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalDenseMemoryUsage() : 0;
1143:         matrix_mem[1] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalLowRankMemoryUsage() : 0;
1144:     #if defined(PETSC_HAVE_CUDA)
1145:         matrix_mem[2] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalDenseMemoryUsage() : 0;
1146:         matrix_mem[3] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalLowRankMemoryUsage() : 0;
1147:     #endif
1148:   #endif
1149:         PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, matrix_mem, rsize, MPI_DOUBLE_PRECISION, MPI_SUM, PetscObjectComm((PetscObject)A)));
1150:         PetscCall(PetscViewerASCIIPrintf(view, "  Memory consumption GB (CPU): %g (dense) %g (low rank) %g (total)\n", matrix_mem[0], matrix_mem[1], matrix_mem[0] + matrix_mem[1]));
1151:   #if defined(PETSC_HAVE_CUDA)
1152:         PetscCall(PetscViewerASCIIPrintf(view, "  Memory consumption GB (GPU): %g (dense) %g (low rank) %g (total)\n", matrix_mem[2], matrix_mem[3], matrix_mem[2] + matrix_mem[3]));
1153:   #endif
1154:       }
1155:     }
1156:   }
1157:   vieweps = PETSC_FALSE;
1158:   PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps", &vieweps, NULL));
1159:   if (vieweps) {
1160:     char        filename[256];
1161:     const char *name;

1163:     PetscCall(PetscObjectGetName((PetscObject)A, &name));
1164:     PetscCall(PetscSNPrintf(filename, sizeof(filename), "%s_structure.eps", name));
1165:     PetscCall(PetscOptionsGetString(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps_filename", filename, sizeof(filename), NULL));
1166:     outputEps(*h2opus->hmatrix, filename);
1167:   }
1168:   PetscFunctionReturn(PETSC_SUCCESS);
1169: }

1171: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat A, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, void *kernelctx)
1172: {
1173:   Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1174:   PetscReal  *gcoords;
1175:   PetscInt    N;
1176:   MPI_Comm    comm;
1177:   PetscMPIInt size;
1178:   PetscBool   cong;

1180:   PetscFunctionBegin;
1181:   PetscCall(PetscLayoutSetUp(A->rmap));
1182:   PetscCall(PetscLayoutSetUp(A->cmap));
1183:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1184:   PetscCall(MatHasCongruentLayouts(A, &cong));
1185:   PetscCheck(cong, comm, PETSC_ERR_SUP, "Only for square matrices with congruent layouts");
1186:   N = A->rmap->N;
1187:   PetscCallMPI(MPI_Comm_size(comm, &size));
1188:   if (spacedim > 0 && size > 1 && cdist) {
1189:     PetscSF      sf;
1190:     MPI_Datatype dtype;

1192:     PetscCallMPI(MPI_Type_contiguous(spacedim, MPIU_REAL, &dtype));
1193:     PetscCallMPI(MPI_Type_commit(&dtype));

1195:     PetscCall(PetscSFCreate(comm, &sf));
1196:     PetscCall(PetscSFSetGraphWithPattern(sf, A->rmap, PETSCSF_PATTERN_ALLGATHER));
1197:     PetscCall(PetscMalloc1(spacedim * N, &gcoords));
1198:     PetscCall(PetscSFBcastBegin(sf, dtype, coords, gcoords, MPI_REPLACE));
1199:     PetscCall(PetscSFBcastEnd(sf, dtype, coords, gcoords, MPI_REPLACE));
1200:     PetscCall(PetscSFDestroy(&sf));
1201:     PetscCallMPI(MPI_Type_free(&dtype));
1202:   } else gcoords = (PetscReal *)coords;

1204:   delete h2opus->ptcloud;
1205:   delete h2opus->kernel;
1206:   h2opus->ptcloud = new PetscPointCloud<PetscReal>(spacedim, N, gcoords);
1207:   if (kernel) h2opus->kernel = new PetscFunctionGenerator<PetscScalar>(kernel, spacedim, kernelctx);
1208:   if (gcoords != coords) PetscCall(PetscFree(gcoords));
1209:   A->preallocated = PETSC_TRUE;
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213:   #if defined(PETSC_H2OPUS_USE_GPU)
1214: static PetscErrorCode MatBindToCPU_H2OPUS(Mat A, PetscBool flg)
1215: {
1216:   PetscMPIInt size;
1217:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

1219:   PetscFunctionBegin;
1220:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1221:   if (flg && A->offloadmask == PETSC_OFFLOAD_GPU) {
1222:     if (size > 1) {
1223:       PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1224:     #if defined(H2OPUS_USE_MPI)
1225:       if (!a->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix_gpu);
1226:       else *a->dist_hmatrix = *a->dist_hmatrix_gpu;
1227:     #endif
1228:     } else {
1229:       PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1230:       if (!a->hmatrix) a->hmatrix = new HMatrix(*a->hmatrix_gpu);
1231:       else *a->hmatrix = *a->hmatrix_gpu;
1232:     }
1233:     delete a->hmatrix_gpu;
1234:     delete a->dist_hmatrix_gpu;
1235:     a->hmatrix_gpu      = NULL;
1236:     a->dist_hmatrix_gpu = NULL;
1237:     A->offloadmask      = PETSC_OFFLOAD_CPU;
1238:   } else if (!flg && A->offloadmask == PETSC_OFFLOAD_CPU) {
1239:     if (size > 1) {
1240:       PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1241:     #if defined(H2OPUS_USE_MPI)
1242:       if (!a->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
1243:       else *a->dist_hmatrix_gpu = *a->dist_hmatrix;
1244:     #endif
1245:     } else {
1246:       PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1247:       if (!a->hmatrix_gpu) a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
1248:       else *a->hmatrix_gpu = *a->hmatrix;
1249:     }
1250:     delete a->hmatrix;
1251:     delete a->dist_hmatrix;
1252:     a->hmatrix      = NULL;
1253:     a->dist_hmatrix = NULL;
1254:     A->offloadmask  = PETSC_OFFLOAD_GPU;
1255:   }
1256:   PetscCall(PetscFree(A->defaultvectype));
1257:   if (!flg) {
1258:     PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1259:   } else {
1260:     PetscCall(PetscStrallocpy(VECSTANDARD, &A->defaultvectype));
1261:   }
1262:   A->boundtocpu = flg;
1263:   PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1265:   #endif

1267: /*MC
1268:    MATH2OPUS = "h2opus" - A matrix type for hierarchical matrices using the H2Opus package {cite}`zampinibouakaramturkiyyahkniokeyes2022`.

1270:    Options Database Key:
1271: .  -mat_type h2opus - matrix type to "h2opus"

1273:    Level: beginner

1275:    Notes:
1276:    H2Opus implements hierarchical matrices in the $H^2$ flavour. It supports CPU or NVIDIA GPUs.

1278:    For CPU only builds, use `./configure --download-h2opus --download-thrust` to install PETSc to use H2Opus.
1279:    In order to run on NVIDIA GPUs, use `./configure --download-h2opus --download-magma --download-kblas`.

1281: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATHTOOL`, `MATDENSE`, `MatCreateH2OpusFromKernel()`, `MatCreateH2OpusFromMat()`
1282: M*/
1283: PETSC_EXTERN PetscErrorCode MatCreate_H2OPUS(Mat A)
1284: {
1285:   Mat_H2OPUS *a;
1286:   PetscMPIInt size;

1288:   PetscFunctionBegin;
1289:   #if defined(PETSC_H2OPUS_USE_GPU)
1290:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
1291:   #endif
1292:   PetscCall(PetscNew(&a));
1293:   A->data = (void *)a;

1295:   a->eta              = 0.9;
1296:   a->leafsize         = 32;
1297:   a->basisord         = 4;
1298:   a->max_rank         = 64;
1299:   a->bs               = 32;
1300:   a->rtol             = 1.e-4;
1301:   a->s                = 1.0;
1302:   a->norm_max_samples = 10;
1303:   a->resize           = PETSC_TRUE; /* reallocate after compression */
1304:   #if defined(H2OPUS_USE_MPI)
1305:   h2opusCreateDistributedHandleComm(&a->handle, PetscObjectComm((PetscObject)A));
1306:   #else
1307:   h2opusCreateHandle(&a->handle);
1308:   #endif
1309:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1310:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATH2OPUS));
1311:   PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));

1313:   A->ops->destroy          = MatDestroy_H2OPUS;
1314:   A->ops->view             = MatView_H2OPUS;
1315:   A->ops->assemblyend      = MatAssemblyEnd_H2OPUS;
1316:   A->ops->mult             = MatMult_H2OPUS;
1317:   A->ops->multtranspose    = MatMultTranspose_H2OPUS;
1318:   A->ops->multadd          = MatMultAdd_H2OPUS;
1319:   A->ops->multtransposeadd = MatMultTransposeAdd_H2OPUS;
1320:   A->ops->scale            = MatScale_H2OPUS;
1321:   A->ops->duplicate        = MatDuplicate_H2OPUS;
1322:   A->ops->setfromoptions   = MatSetFromOptions_H2OPUS;
1323:   A->ops->norm             = MatNorm_H2OPUS;
1324:   A->ops->zeroentries      = MatZeroEntries_H2OPUS;
1325:   #if defined(PETSC_H2OPUS_USE_GPU)
1326:   A->ops->bindtocpu = MatBindToCPU_H2OPUS;
1327:   #endif

1329:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", MatProductSetFromOptions_H2OPUS));
1330:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", MatProductSetFromOptions_H2OPUS));
1331:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", MatProductSetFromOptions_H2OPUS));
1332:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", MatProductSetFromOptions_H2OPUS));
1333:   #if defined(PETSC_H2OPUS_USE_GPU)
1334:   PetscCall(PetscFree(A->defaultvectype));
1335:   PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1336:   #endif
1337:   PetscFunctionReturn(PETSC_SUCCESS);
1338: }

1340: /*@
1341:   MatH2OpusOrthogonalize - Orthogonalize the basis tree of a hierarchical matrix.

1343:   Input Parameter:
1344: . A - the matrix

1346:   Level: intermediate

1348: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`
1349: @*/
1350: PetscErrorCode MatH2OpusOrthogonalize(Mat A)
1351: {
1352:   PetscBool   ish2opus;
1353:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1354:   PetscMPIInt size;
1355:   PetscBool   boundtocpu = PETSC_TRUE;

1357:   PetscFunctionBegin;
1360:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1361:   if (!ish2opus) PetscFunctionReturn(PETSC_SUCCESS);
1362:   if (a->orthogonal) PetscFunctionReturn(PETSC_SUCCESS);
1363:   HLibProfile::clear();
1364:   PetscCall(PetscLogEventBegin(MAT_H2Opus_Orthog, A, 0, 0, 0));
1365:   #if defined(PETSC_H2OPUS_USE_GPU)
1366:   boundtocpu = A->boundtocpu;
1367:   #endif
1368:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1369:   if (size > 1) {
1370:     if (boundtocpu) {
1371:       PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1372:   #if defined(H2OPUS_USE_MPI)
1373:       distributed_horthog(*a->dist_hmatrix, a->handle);
1374:   #endif
1375:   #if defined(PETSC_H2OPUS_USE_GPU)
1376:       A->offloadmask = PETSC_OFFLOAD_CPU;
1377:     } else {
1378:       PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1379:       PetscCall(PetscLogGpuTimeBegin());
1380:     #if defined(H2OPUS_USE_MPI)
1381:       distributed_horthog(*a->dist_hmatrix_gpu, a->handle);
1382:     #endif
1383:       PetscCall(PetscLogGpuTimeEnd());
1384:   #endif
1385:     }
1386:   } else {
1387:   #if defined(H2OPUS_USE_MPI)
1388:     h2opusHandle_t handle = a->handle->handle;
1389:   #else
1390:     h2opusHandle_t handle = a->handle;
1391:   #endif
1392:     if (boundtocpu) {
1393:       PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1394:       horthog(*a->hmatrix, handle);
1395:   #if defined(PETSC_H2OPUS_USE_GPU)
1396:       A->offloadmask = PETSC_OFFLOAD_CPU;
1397:     } else {
1398:       PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1399:       PetscCall(PetscLogGpuTimeBegin());
1400:       horthog(*a->hmatrix_gpu, handle);
1401:       PetscCall(PetscLogGpuTimeEnd());
1402:   #endif
1403:     }
1404:   }
1405:   a->orthogonal = PETSC_TRUE;
1406:   { /* log flops */
1407:     double gops, time, perf, dev;
1408:     HLibProfile::getHorthogPerf(gops, time, perf, dev);
1409:   #if defined(PETSC_H2OPUS_USE_GPU)
1410:     if (boundtocpu) {
1411:       PetscCall(PetscLogFlops(1e9 * gops));
1412:     } else {
1413:       PetscCall(PetscLogGpuFlops(1e9 * gops));
1414:     }
1415:   #else
1416:     PetscCall(PetscLogFlops(1e9 * gops));
1417:   #endif
1418:   }
1419:   PetscCall(PetscLogEventEnd(MAT_H2Opus_Orthog, A, 0, 0, 0));
1420:   PetscFunctionReturn(PETSC_SUCCESS);
1421: }

1423: /*@
1424:   MatH2OpusCompress - Compress a hierarchical matrix.

1426:   Input Parameters:
1427: + A   - the matrix
1428: - tol - the absolute truncation threshold

1430:   Level: intermediate

1432: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusOrthogonalize()`
1433: @*/
1434: PetscErrorCode MatH2OpusCompress(Mat A, PetscReal tol)
1435: {
1436:   PetscBool   ish2opus;
1437:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1438:   PetscMPIInt size;
1439:   PetscBool   boundtocpu = PETSC_TRUE;

1441:   PetscFunctionBegin;
1445:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1446:   if (!ish2opus || tol <= 0.0) PetscFunctionReturn(PETSC_SUCCESS);
1447:   PetscCall(MatH2OpusOrthogonalize(A));
1448:   HLibProfile::clear();
1449:   PetscCall(PetscLogEventBegin(MAT_H2Opus_Compress, A, 0, 0, 0));
1450:   #if defined(PETSC_H2OPUS_USE_GPU)
1451:   boundtocpu = A->boundtocpu;
1452:   #endif
1453:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1454:   if (size > 1) {
1455:     if (boundtocpu) {
1456:       PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1457:   #if defined(H2OPUS_USE_MPI)
1458:       distributed_hcompress(*a->dist_hmatrix, tol, a->handle);
1459:       if (a->resize) {
1460:         DistributedHMatrix *dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix);
1461:         delete a->dist_hmatrix;
1462:         a->dist_hmatrix = dist_hmatrix;
1463:       }
1464:   #endif
1465:   #if defined(PETSC_H2OPUS_USE_GPU)
1466:       A->offloadmask = PETSC_OFFLOAD_CPU;
1467:     } else {
1468:       PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1469:       PetscCall(PetscLogGpuTimeBegin());
1470:     #if defined(H2OPUS_USE_MPI)
1471:       distributed_hcompress(*a->dist_hmatrix_gpu, tol, a->handle);

1473:       if (a->resize) {
1474:         DistributedHMatrix_GPU *dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix_gpu);
1475:         delete a->dist_hmatrix_gpu;
1476:         a->dist_hmatrix_gpu = dist_hmatrix_gpu;
1477:       }
1478:     #endif
1479:       PetscCall(PetscLogGpuTimeEnd());
1480:   #endif
1481:     }
1482:   } else {
1483:   #if defined(H2OPUS_USE_MPI)
1484:     h2opusHandle_t handle = a->handle->handle;
1485:   #else
1486:     h2opusHandle_t handle = a->handle;
1487:   #endif
1488:     if (boundtocpu) {
1489:       PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1490:       hcompress(*a->hmatrix, tol, handle);

1492:       if (a->resize) {
1493:         HMatrix *hmatrix = new HMatrix(*a->hmatrix);
1494:         delete a->hmatrix;
1495:         a->hmatrix = hmatrix;
1496:       }
1497:   #if defined(PETSC_H2OPUS_USE_GPU)
1498:       A->offloadmask = PETSC_OFFLOAD_CPU;
1499:     } else {
1500:       PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1501:       PetscCall(PetscLogGpuTimeBegin());
1502:       hcompress(*a->hmatrix_gpu, tol, handle);
1503:       PetscCall(PetscLogGpuTimeEnd());

1505:       if (a->resize) {
1506:         HMatrix_GPU *hmatrix_gpu = new HMatrix_GPU(*a->hmatrix_gpu);
1507:         delete a->hmatrix_gpu;
1508:         a->hmatrix_gpu = hmatrix_gpu;
1509:       }
1510:   #endif
1511:     }
1512:   }
1513:   { /* log flops */
1514:     double gops, time, perf, dev;
1515:     HLibProfile::getHcompressPerf(gops, time, perf, dev);
1516:   #if defined(PETSC_H2OPUS_USE_GPU)
1517:     if (boundtocpu) {
1518:       PetscCall(PetscLogFlops(1e9 * gops));
1519:     } else {
1520:       PetscCall(PetscLogGpuFlops(1e9 * gops));
1521:     }
1522:   #else
1523:     PetscCall(PetscLogFlops(1e9 * gops));
1524:   #endif
1525:   }
1526:   PetscCall(PetscLogEventEnd(MAT_H2Opus_Compress, A, 0, 0, 0));
1527:   PetscFunctionReturn(PETSC_SUCCESS);
1528: }

1530: /*@
1531:   MatH2OpusSetSamplingMat - Set a matrix to be sampled from matrix-vector products on another matrix to construct a hierarchical matrix.

1533:   Input Parameters:
1534: + A   - the hierarchical matrix
1535: . B   - the matrix to be sampled
1536: . bs  - maximum number of samples to be taken concurrently
1537: - tol - relative tolerance for construction

1539:   Level: intermediate

1541:   Notes:
1542:   You need to call `MatAssemblyBegin()` and `MatAssemblyEnd()` to update the hierarchical matrix.

1544: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`
1545: @*/
1546: PetscErrorCode MatH2OpusSetSamplingMat(Mat A, Mat B, PetscInt bs, PetscReal tol)
1547: {
1548:   PetscBool ish2opus;

1550:   PetscFunctionBegin;
1556:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1557:   if (ish2opus) {
1558:     Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

1560:     if (!a->sampler) a->sampler = new PetscMatrixSampler();
1561:     a->sampler->SetSamplingMat(B);
1562:     if (bs > 0) a->bs = bs;
1563:     if (tol > 0.) a->rtol = tol;
1564:     delete a->kernel;
1565:   }
1566:   PetscFunctionReturn(PETSC_SUCCESS);
1567: }

1569: /*@C
1570:   MatCreateH2OpusFromKernel - Creates a `MATH2OPUS` from a user-supplied kernel.

1572:   Input Parameters:
1573: + comm      - MPI communicator
1574: . m         - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1575: . n         - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1576: . M         - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1577: . N         - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1578: . spacedim  - dimension of the space coordinates
1579: . coords    - coordinates of the points
1580: . cdist     - whether or not coordinates are distributed
1581: . kernel    - computational kernel (or `NULL`)
1582: . kernelctx - kernel context
1583: . eta       - admissibility condition tolerance
1584: . leafsize  - leaf size in cluster tree
1585: - basisord  - approximation order for Chebychev interpolation of low-rank blocks

1587:   Output Parameter:
1588: . nA - matrix

1590:   Options Database Keys:
1591: + -mat_h2opus_leafsize <`PetscInt`>    - Leaf size of cluster tree
1592: . -mat_h2opus_eta <`PetscReal`>        - Admissibility condition tolerance
1593: . -mat_h2opus_order <`PetscInt`>       - Chebychev approximation order
1594: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be used when estimating norms

1596:   Level: intermediate

1598: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`
1599: @*/
1600: PetscErrorCode MatCreateH2OpusFromKernel(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, void *kernelctx, PetscReal eta, PetscInt leafsize, PetscInt basisord, Mat *nA)
1601: {
1602:   Mat         A;
1603:   Mat_H2OPUS *h2opus;
1604:   #if defined(PETSC_H2OPUS_USE_GPU)
1605:   PetscBool iscpu = PETSC_FALSE;
1606:   #else
1607:   PetscBool iscpu = PETSC_TRUE;
1608:   #endif

1610:   PetscFunctionBegin;
1611:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1612:   PetscCall(MatCreate(comm, &A));
1613:   PetscCall(MatSetSizes(A, m, n, M, N));
1614:   PetscCheck(M == N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1615:   PetscCall(MatSetType(A, MATH2OPUS));
1616:   PetscCall(MatBindToCPU(A, iscpu));
1617:   PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, kernel, kernelctx));

1619:   h2opus = (Mat_H2OPUS *)A->data;
1620:   if (eta > 0.) h2opus->eta = eta;
1621:   if (leafsize > 0) h2opus->leafsize = leafsize;
1622:   if (basisord > 0) h2opus->basisord = basisord;

1624:   *nA = A;
1625:   PetscFunctionReturn(PETSC_SUCCESS);
1626: }

1628: /*@
1629:   MatCreateH2OpusFromMat - Creates a `MATH2OPUS` sampling from a user-supplied operator.

1631:   Input Parameters:
1632: + B        - the matrix to be sampled
1633: . spacedim - dimension of the space coordinates
1634: . coords   - coordinates of the points
1635: . cdist    - whether or not coordinates are distributed
1636: . eta      - admissibility condition tolerance
1637: . leafsize - leaf size in cluster tree
1638: . maxrank  - maximum rank allowed
1639: . bs       - maximum number of samples to be taken concurrently
1640: - rtol     - relative tolerance for construction

1642:   Output Parameter:
1643: . nA - matrix

1645:   Options Database Keys:
1646: + -mat_h2opus_leafsize <`PetscInt`>      - Leaf size of cluster tree
1647: . -mat_h2opus_eta <`PetscReal`>          - Admissibility condition tolerance
1648: . -mat_h2opus_maxrank <`PetscInt`>       - Maximum rank when constructed from matvecs
1649: . -mat_h2opus_samples <`PetscInt`>       - Maximum number of samples to be taken concurrently when constructing from matvecs
1650: . -mat_h2opus_rtol <`PetscReal`>         - Relative tolerance for construction from sampling
1651: . -mat_h2opus_check <`PetscBool`>        - Check error when constructing from sampling during MatAssemblyEnd()
1652: . -mat_h2opus_hara_verbose <`PetscBool`> - Verbose output from hara construction
1653: - -mat_h2opus_normsamples <`PetscInt`>   - Maximum number of samples to be when estimating norms

1655:   Level: intermediate

1657:   Note:
1658:   Not available in parallel

1660: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
1661: @*/
1662: PetscErrorCode MatCreateH2OpusFromMat(Mat B, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, PetscReal eta, PetscInt leafsize, PetscInt maxrank, PetscInt bs, PetscReal rtol, Mat *nA)
1663: {
1664:   Mat         A;
1665:   Mat_H2OPUS *h2opus;
1666:   MPI_Comm    comm;
1667:   PetscBool   boundtocpu = PETSC_TRUE;

1669:   PetscFunctionBegin;
1678:   PetscAssertPointer(nA, 10);
1679:   PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1680:   PetscCheck(B->rmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1681:   PetscCheck(B->rmap->N == B->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1682:   PetscCall(MatCreate(comm, &A));
1683:   PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1684:   #if defined(PETSC_H2OPUS_USE_GPU)
1685:   {
1686:     VecType   vtype;
1687:     PetscBool isstd, iscuda, iskok;

1689:     PetscCall(MatGetVecType(B, &vtype));
1690:     PetscCall(PetscStrcmpAny(vtype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
1691:     PetscCall(PetscStrcmpAny(vtype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
1692:     PetscCall(PetscStrcmpAny(vtype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
1693:     PetscCheck(isstd || iscuda || iskok, comm, PETSC_ERR_SUP, "Not for type %s", vtype);
1694:     if (iscuda && !B->boundtocpu) boundtocpu = PETSC_FALSE;
1695:     if (iskok && PetscDefined(HAVE_MACRO_KOKKOS_ENABLE_CUDA)) boundtocpu = PETSC_FALSE;
1696:   }
1697:   #endif
1698:   PetscCall(MatSetType(A, MATH2OPUS));
1699:   PetscCall(MatBindToCPU(A, boundtocpu));
1700:   if (spacedim) PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, NULL, NULL));
1701:   PetscCall(MatPropagateSymmetryOptions(B, A));
1702:   /* PetscCheck(A->symmetric,comm,PETSC_ERR_SUP,"Unsymmetric sampling does not work"); */

1704:   h2opus          = (Mat_H2OPUS *)A->data;
1705:   h2opus->sampler = new PetscMatrixSampler(B);
1706:   if (eta > 0.) h2opus->eta = eta;
1707:   if (leafsize > 0) h2opus->leafsize = leafsize;
1708:   if (maxrank > 0) h2opus->max_rank = maxrank;
1709:   if (bs > 0) h2opus->bs = bs;
1710:   if (rtol > 0.) h2opus->rtol = rtol;
1711:   *nA             = A;
1712:   A->preallocated = PETSC_TRUE;
1713:   PetscFunctionReturn(PETSC_SUCCESS);
1714: }

1716: /*@
1717:   MatH2OpusGetIndexMap - Access reordering index set.

1719:   Input Parameter:
1720: . A - the matrix

1722:   Output Parameter:
1723: . indexmap - the index set for the reordering

1725:   Level: intermediate

1727: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1728: @*/
1729: PetscErrorCode MatH2OpusGetIndexMap(Mat A, IS *indexmap)
1730: {
1731:   PetscBool   ish2opus;
1732:   Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;

1734:   PetscFunctionBegin;
1737:   PetscAssertPointer(indexmap, 2);
1738:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1739:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1740:   PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1741:   *indexmap = a->h2opus_indexmap;
1742:   PetscFunctionReturn(PETSC_SUCCESS);
1743: }

1745: /*@
1746:   MatH2OpusMapVec - Maps a vector between PETSc and H2Opus ordering

1748:   Input Parameters:
1749: + A             - the matrix
1750: . nativetopetsc - if true, maps from H2Opus ordering to PETSc ordering. If false, applies the reverse map
1751: - in            - the vector to be mapped

1753:   Output Parameter:
1754: . out - the newly created mapped vector

1756:   Level: intermediate

1758: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1759: @*/
1760: PetscErrorCode MatH2OpusMapVec(Mat A, PetscBool nativetopetsc, Vec in, Vec *out)
1761: {
1762:   PetscBool    ish2opus;
1763:   Mat_H2OPUS  *a = (Mat_H2OPUS *)A->data;
1764:   PetscScalar *xin, *xout;
1765:   PetscBool    nm;

1767:   PetscFunctionBegin;
1772:   PetscAssertPointer(out, 4);
1773:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1774:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1775:   PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1776:   nm = a->nativemult;
1777:   PetscCall(MatH2OpusSetNativeMult(A, (PetscBool)!nativetopetsc));
1778:   PetscCall(MatCreateVecs(A, out, NULL));
1779:   PetscCall(MatH2OpusSetNativeMult(A, nm));
1780:   if (!a->sf) { /* same ordering */
1781:     PetscCall(VecCopy(in, *out));
1782:     PetscFunctionReturn(PETSC_SUCCESS);
1783:   }
1784:   PetscCall(VecGetArrayRead(in, (const PetscScalar **)&xin));
1785:   PetscCall(VecGetArrayWrite(*out, &xout));
1786:   if (nativetopetsc) {
1787:     PetscCall(PetscSFReduceBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1788:     PetscCall(PetscSFReduceEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1789:   } else {
1790:     PetscCall(PetscSFBcastBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1791:     PetscCall(PetscSFBcastEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1792:   }
1793:   PetscCall(VecRestoreArrayRead(in, (const PetscScalar **)&xin));
1794:   PetscCall(VecRestoreArrayWrite(*out, &xout));
1795:   PetscFunctionReturn(PETSC_SUCCESS);
1796: }

1798: /*@
1799:   MatH2OpusLowRankUpdate - Perform a low-rank update of the form $ A = A + s * U * V^T $

1801:   Input Parameters:
1802: + A - the hierarchical `MATH2OPUS` matrix
1803: . s - the scaling factor
1804: . U - the dense low-rank update matrix
1805: - V - (optional) the dense low-rank update matrix (if `NULL`, then `V` = `U` is assumed)

1807:   Note:
1808:   The `U` and `V` matrices must be in `MATDENSE` dense format

1810:   Level: intermediate

1812: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`, `MATDENSE`
1813: @*/
1814: PetscErrorCode MatH2OpusLowRankUpdate(Mat A, Mat U, Mat V, PetscScalar s)
1815: {
1816:   PetscBool flg;

1818:   PetscFunctionBegin;
1821:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1823:   PetscCheckSameComm(A, 1, U, 2);
1824:   if (V) {
1826:     PetscCheckSameComm(A, 1, V, 3);
1827:   }

1830:   if (!V) V = U;
1831:   PetscCheck(U->cmap->N == V->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Non matching rank update %" PetscInt_FMT " != %" PetscInt_FMT, U->cmap->N, V->cmap->N);
1832:   if (!U->cmap->N) PetscFunctionReturn(PETSC_SUCCESS);
1833:   PetscCall(PetscLayoutCompare(U->rmap, A->rmap, &flg));
1834:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A and U must have the same row layout");
1835:   PetscCall(PetscLayoutCompare(V->rmap, A->cmap, &flg));
1836:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A column layout must match V row column layout");
1837:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &flg));
1838:   if (flg) {
1839:     Mat_H2OPUS        *a = (Mat_H2OPUS *)A->data;
1840:     const PetscScalar *u, *v, *uu, *vv;
1841:     PetscInt           ldu, ldv;
1842:     PetscMPIInt        size;
1843:   #if defined(H2OPUS_USE_MPI)
1844:     h2opusHandle_t handle = a->handle->handle;
1845:   #else
1846:     h2opusHandle_t handle = a->handle;
1847:   #endif
1848:     PetscBool usesf = (PetscBool)(a->sf && !a->nativemult);
1849:     PetscSF   usf, vsf;

1851:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1852:     PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet implemented in parallel");
1853:     PetscCall(PetscLogEventBegin(MAT_H2Opus_LR, A, 0, 0, 0));
1854:     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)U, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1855:     PetscCheck(flg, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "Not for U of type %s", ((PetscObject)U)->type_name);
1856:     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)V, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1857:     PetscCheck(flg, PetscObjectComm((PetscObject)V), PETSC_ERR_SUP, "Not for V of type %s", ((PetscObject)V)->type_name);
1858:     PetscCall(MatDenseGetLDA(U, &ldu));
1859:     PetscCall(MatDenseGetLDA(V, &ldv));
1860:     PetscCall(MatBoundToCPU(A, &flg));
1861:     if (usesf) {
1862:       PetscInt n;

1864:       PetscCall(MatDenseGetH2OpusStridedSF(U, a->sf, &usf));
1865:       PetscCall(MatDenseGetH2OpusStridedSF(V, a->sf, &vsf));
1866:       PetscCall(MatH2OpusResizeBuffers_Private(A, U->cmap->N, V->cmap->N));
1867:       PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
1868:       ldu = n;
1869:       ldv = n;
1870:     }
1871:     if (flg) {
1872:       PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1873:       PetscCall(MatDenseGetArrayRead(U, &u));
1874:       PetscCall(MatDenseGetArrayRead(V, &v));
1875:       if (usesf) {
1876:         vv = MatH2OpusGetThrustPointer(*a->yy);
1877:         PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1878:         PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1879:         if (U != V) {
1880:           uu = MatH2OpusGetThrustPointer(*a->xx);
1881:           PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1882:           PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1883:         } else uu = vv;
1884:       } else {
1885:         uu = u;
1886:         vv = v;
1887:       }
1888:       hlru_global(*a->hmatrix, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1889:       PetscCall(MatDenseRestoreArrayRead(U, &u));
1890:       PetscCall(MatDenseRestoreArrayRead(V, &v));
1891:     } else {
1892:   #if defined(PETSC_H2OPUS_USE_GPU)
1893:       PetscBool flgU, flgV;

1895:       PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1896:       PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &flgU, MATSEQDENSE, MATMPIDENSE, ""));
1897:       if (flgU) PetscCall(MatConvert(U, MATDENSECUDA, MAT_INPLACE_MATRIX, &U));
1898:       PetscCall(PetscObjectTypeCompareAny((PetscObject)V, &flgV, MATSEQDENSE, MATMPIDENSE, ""));
1899:       if (flgV) PetscCall(MatConvert(V, MATDENSECUDA, MAT_INPLACE_MATRIX, &V));
1900:       PetscCall(MatDenseCUDAGetArrayRead(U, &u));
1901:       PetscCall(MatDenseCUDAGetArrayRead(V, &v));
1902:       if (usesf) {
1903:         vv = MatH2OpusGetThrustPointer(*a->yy_gpu);
1904:         PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1905:         PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1906:         if (U != V) {
1907:           uu = MatH2OpusGetThrustPointer(*a->xx_gpu);
1908:           PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1909:           PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1910:         } else uu = vv;
1911:       } else {
1912:         uu = u;
1913:         vv = v;
1914:       }
1915:   #else
1916:       SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
1917:   #endif
1918:       hlru_global(*a->hmatrix_gpu, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1919:   #if defined(PETSC_H2OPUS_USE_GPU)
1920:       PetscCall(MatDenseCUDARestoreArrayRead(U, &u));
1921:       PetscCall(MatDenseCUDARestoreArrayRead(V, &v));
1922:       if (flgU) PetscCall(MatConvert(U, MATDENSE, MAT_INPLACE_MATRIX, &U));
1923:       if (flgV) PetscCall(MatConvert(V, MATDENSE, MAT_INPLACE_MATRIX, &V));
1924:   #endif
1925:     }
1926:     PetscCall(PetscLogEventEnd(MAT_H2Opus_LR, A, 0, 0, 0));
1927:     a->orthogonal = PETSC_FALSE;
1928:   }
1929:   PetscFunctionReturn(PETSC_SUCCESS);
1930: }
1931: #endif