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 *);
21: #define MatH2OpusGetThrustPointer(v) thrust::raw_pointer_cast((v).data())
23: /* Use GPU only if H2OPUS is configured for GPU */
24: #if PetscDefined(HAVE_CUDA) && defined(H2OPUS_USE_GPU)
25: #define PETSC_H2OPUS_USE_GPU
26: #endif
27: #if PetscDefined(H2OPUS_USE_GPU)
28: #define MatH2OpusUpdateIfNeeded(A, B) MatBindToCPU(A, (PetscBool)((A)->boundtocpu || (B)))
29: #else
30: #define MatH2OpusUpdateIfNeeded(A, B) PETSC_SUCCESS
31: #endif
33: // TODO H2OPUS:
34: // DistributedHMatrix
35: // unsymmetric ?
36: // transpose for distributed_hgemv?
37: // clearData()
38: // Unify interface for sequential and parallel?
39: // Reuse geometric construction (almost possible, only the unsymmetric case is explicitly handled)
40: //
41: template <class T>
42: class PetscPointCloud : public H2OpusDataSet<T> {
43: private:
44: int dimension;
45: size_t num_points;
46: std::vector<T> pts;
48: public:
49: PetscPointCloud(int dim, size_t num_pts, const T coords[])
50: {
51: dim = dim > 0 ? dim : 1;
52: this->dimension = dim;
53: this->num_points = num_pts;
55: pts.resize(num_pts * dim);
56: if (coords) {
57: for (size_t n = 0; n < num_pts; n++)
58: for (int i = 0; i < dim; i++) pts[n * dim + i] = coords[n * dim + i];
59: } else {
60: PetscReal h = 1.0; //num_pts > 1 ? 1./(num_pts - 1) : 0.0;
61: for (size_t n = 0; n < num_pts; n++) {
62: pts[n * dim] = n * h;
63: for (int i = 1; i < dim; i++) pts[n * dim + i] = 0.0;
64: }
65: }
66: }
68: PetscPointCloud(const PetscPointCloud<T> &other)
69: {
70: size_t N = other.dimension * other.num_points;
71: this->dimension = other.dimension;
72: this->num_points = other.num_points;
73: this->pts.resize(N);
74: for (size_t i = 0; i < N; i++) this->pts[i] = other.pts[i];
75: }
77: int getDimension() const { return dimension; }
79: size_t getDataSetSize() const { return num_points; }
81: T getDataPoint(size_t idx, int dim) const
82: {
83: assert(dim < dimension && idx < num_points);
84: return pts[idx * dimension + dim];
85: }
87: void Print(std::ostream &out = std::cout)
88: {
89: out << "Dimension: " << dimension << std::endl;
90: out << "NumPoints: " << num_points << std::endl;
91: for (size_t n = 0; n < num_points; n++) {
92: for (int d = 0; d < dimension; d++) out << pts[n * dimension + d] << " ";
93: out << std::endl;
94: }
95: }
96: };
98: template <class T>
99: class PetscFunctionGenerator {
100: private:
101: MatH2OpusKernelFn *k;
102: int dim;
103: void *ctx;
105: public:
106: PetscFunctionGenerator(MatH2OpusKernelFn *k, int dim, PetscCtx ctx)
107: {
108: this->k = k;
109: this->dim = dim;
110: this->ctx = ctx;
111: }
112: PetscFunctionGenerator(PetscFunctionGenerator &other)
113: {
114: this->k = other.k;
115: this->dim = other.dim;
116: this->ctx = other.ctx;
117: }
118: T operator()(PetscReal *pt1, PetscReal *pt2) { return (T)((*this->k)(this->dim, pt1, pt2, this->ctx)); }
119: };
121: #include <../src/mat/impls/h2opus/math2opussampler.hpp>
123: /* just to not clutter the code */
124: #if !defined(H2OPUS_USE_GPU)
125: typedef HMatrix HMatrix_GPU;
126: #if defined(H2OPUS_USE_MPI)
127: typedef DistributedHMatrix DistributedHMatrix_GPU;
128: #endif
129: #endif
131: typedef struct {
132: #if defined(H2OPUS_USE_MPI)
133: distributedH2OpusHandle_t handle;
134: #else
135: h2opusHandle_t handle;
136: #endif
137: /* Sequential and parallel matrices are two different classes at the moment */
138: HMatrix *hmatrix;
139: #if defined(H2OPUS_USE_MPI)
140: DistributedHMatrix *dist_hmatrix;
141: #else
142: HMatrix *dist_hmatrix; /* just to not clutter the code */
143: #endif
144: /* May use permutations */
145: PetscSF sf;
146: PetscLayout h2opus_rmap, h2opus_cmap;
147: IS h2opus_indexmap;
148: thrust::host_vector<PetscScalar> *xx, *yy;
149: PetscInt xxs, yys;
150: PetscBool multsetup;
152: /* GPU */
153: HMatrix_GPU *hmatrix_gpu;
154: #if defined(H2OPUS_USE_MPI)
155: DistributedHMatrix_GPU *dist_hmatrix_gpu;
156: #else
157: HMatrix_GPU *dist_hmatrix_gpu; /* just to not clutter the code */
158: #endif
159: #if PetscDefined(H2OPUS_USE_GPU)
160: thrust::device_vector<PetscScalar> *xx_gpu, *yy_gpu;
161: PetscInt xxs_gpu, yys_gpu;
162: #endif
164: /* construction from matvecs */
165: PetscMatrixSampler *sampler;
166: PetscBool nativemult;
168: /* Admissibility */
169: PetscReal eta;
170: PetscInt leafsize;
172: /* for dof reordering */
173: PetscPointCloud<PetscReal> *ptcloud;
175: /* kernel for generating matrix entries */
176: PetscFunctionGenerator<PetscScalar> *kernel;
178: /* basis orthogonalized? */
179: PetscBool orthogonal;
181: /* customization */
182: PetscInt basisord;
183: PetscInt max_rank;
184: PetscInt bs;
185: PetscReal rtol;
186: PetscInt norm_max_samples;
187: PetscBool check_construction;
188: PetscBool hara_verbose;
189: PetscBool resize;
191: /* keeps track of MatScale values */
192: PetscScalar s;
193: } Mat_H2OPUS;
195: static PetscErrorCode MatDestroy_H2OPUS(Mat A)
196: {
197: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
199: PetscFunctionBegin;
200: #if defined(H2OPUS_USE_MPI)
201: h2opusDestroyDistributedHandle(a->handle);
202: #else
203: h2opusDestroyHandle(a->handle);
204: #endif
205: delete a->dist_hmatrix;
206: delete a->hmatrix;
207: PetscCall(PetscSFDestroy(&a->sf));
208: PetscCall(PetscLayoutDestroy(&a->h2opus_rmap));
209: PetscCall(PetscLayoutDestroy(&a->h2opus_cmap));
210: PetscCall(ISDestroy(&a->h2opus_indexmap));
211: delete a->xx;
212: delete a->yy;
213: delete a->hmatrix_gpu;
214: delete a->dist_hmatrix_gpu;
215: #if PetscDefined(H2OPUS_USE_GPU)
216: delete a->xx_gpu;
217: delete a->yy_gpu;
218: #endif
219: delete a->sampler;
220: delete a->ptcloud;
221: delete a->kernel;
222: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", NULL));
223: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", NULL));
224: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", NULL));
225: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", NULL));
226: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
227: PetscCall(PetscFree(A->data));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: MatH2OpusSetNativeMult - Enable or disable the native H2Opus matrix-vector multiplication path for a `MATH2OPUS`.
234: Logically Collective
236: Input Parameters:
237: + A - the `MATH2OPUS` matrix
238: - nm - `PETSC_TRUE` to enable the native H2Opus multiply layout, `PETSC_FALSE` to use the PETSc layout
240: Level: advanced
242: Note:
243: Switching this flag swaps the row and column `PetscLayout`s of `A` with those needed by H2Opus so that
244: vectors created by `MatCreateVecs()` are compatible with the currently selected multiplication path.
246: .seealso: `Mat`, `MATH2OPUS`, `MatH2OpusGetNativeMult()`
247: @*/
248: PetscErrorCode MatH2OpusSetNativeMult(Mat A, PetscBool nm)
249: {
250: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
251: PetscBool ish2opus;
253: PetscFunctionBegin;
256: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
257: if (ish2opus) {
258: if (a->h2opus_rmap) { /* need to swap layouts for vector creation */
259: if ((!a->nativemult && nm) || (a->nativemult && !nm)) {
260: PetscLayout t;
261: t = A->rmap;
262: A->rmap = a->h2opus_rmap;
263: a->h2opus_rmap = t;
264: t = A->cmap;
265: A->cmap = a->h2opus_cmap;
266: a->h2opus_cmap = t;
267: }
268: }
269: a->nativemult = nm;
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: MatH2OpusGetNativeMult - Query whether the native H2Opus matrix-vector multiplication path is enabled for a `MATH2OPUS`.
277: Not Collective
279: Input Parameter:
280: . A - the `MATH2OPUS` matrix
282: Output Parameter:
283: . nm - `PETSC_TRUE` if the native H2Opus multiply is enabled, `PETSC_FALSE` otherwise
285: Level: advanced
287: .seealso: `Mat`, `MATH2OPUS`, `MatH2OpusSetNativeMult()`
288: @*/
289: PetscErrorCode MatH2OpusGetNativeMult(Mat A, PetscBool *nm)
290: {
291: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
292: PetscBool ish2opus;
294: PetscFunctionBegin;
296: PetscAssertPointer(nm, 2);
297: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
298: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
299: *nm = a->nativemult;
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: PETSC_EXTERN PetscErrorCode MatNorm_H2OPUS(Mat A, NormType normtype, PetscReal *n)
304: {
305: PetscBool ish2opus;
306: PetscInt nmax = PETSC_DECIDE;
307: Mat_H2OPUS *a = NULL;
308: PetscBool mult = PETSC_FALSE;
310: PetscFunctionBegin;
311: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
312: if (ish2opus) { /* set userdefine number of samples and fastpath for mult (norms are order independent) */
313: a = (Mat_H2OPUS *)A->data;
315: nmax = a->norm_max_samples;
316: mult = a->nativemult;
317: PetscCall(MatH2OpusSetNativeMult(A, PETSC_TRUE));
318: } else {
319: PetscCall(PetscOptionsGetInt(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_approximate_norm_samples", &nmax, NULL));
320: }
321: PetscCall(MatNormApproximate(A, normtype, nmax, n));
322: if (a) PetscCall(MatH2OpusSetNativeMult(A, mult));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode MatH2OpusResizeBuffers_Private(Mat A, PetscInt xN, PetscInt yN)
327: {
328: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
329: PetscInt n;
330: PetscBool boundtocpu = PETSC_TRUE;
332: PetscFunctionBegin;
333: #if PetscDefined(H2OPUS_USE_GPU)
334: boundtocpu = A->boundtocpu;
335: #endif
336: PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
337: if (boundtocpu) {
338: if (h2opus->xxs < xN) {
339: h2opus->xx->resize(n * xN);
340: h2opus->xxs = xN;
341: }
342: if (h2opus->yys < yN) {
343: h2opus->yy->resize(n * yN);
344: h2opus->yys = yN;
345: }
346: }
347: #if PetscDefined(H2OPUS_USE_GPU)
348: if (!boundtocpu) {
349: if (h2opus->xxs_gpu < xN) {
350: h2opus->xx_gpu->resize(n * xN);
351: h2opus->xxs_gpu = xN;
352: }
353: if (h2opus->yys_gpu < yN) {
354: h2opus->yy_gpu->resize(n * yN);
355: h2opus->yys_gpu = yN;
356: }
357: }
358: #endif
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode MatMultNKernel_H2OPUS(Mat A, PetscBool transA, Mat B, Mat C)
363: {
364: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
365: #if defined(H2OPUS_USE_MPI)
366: h2opusHandle_t handle = h2opus->handle->handle;
367: #else
368: h2opusHandle_t handle = h2opus->handle;
369: #endif
370: PetscBool boundtocpu = PETSC_TRUE;
371: PetscScalar *xx, *yy, *uxx, *uyy;
372: PetscInt blda, clda;
373: PetscMPIInt size;
374: PetscSF bsf, csf;
375: PetscBool usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);
377: PetscFunctionBegin;
378: HLibProfile::clear();
379: #if PetscDefined(H2OPUS_USE_GPU)
380: boundtocpu = A->boundtocpu;
381: #endif
382: PetscCall(MatDenseGetLDA(B, &blda));
383: PetscCall(MatDenseGetLDA(C, &clda));
384: if (usesf) {
385: PetscInt n;
387: PetscCall(MatDenseGetH2OpusStridedSF(B, h2opus->sf, &bsf));
388: PetscCall(MatDenseGetH2OpusStridedSF(C, h2opus->sf, &csf));
390: PetscCall(MatH2OpusResizeBuffers_Private(A, B->cmap->N, C->cmap->N));
391: PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
392: blda = n;
393: clda = n;
394: }
395: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
396: if (boundtocpu) {
397: PetscCall(MatDenseGetArrayRead(B, (const PetscScalar **)&xx));
398: PetscCall(MatDenseGetArrayWrite(C, &yy));
399: if (usesf) {
400: uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
401: uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
402: PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
403: PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
404: } else {
405: uxx = xx;
406: uyy = yy;
407: }
408: if (size > 1) {
409: PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
410: PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
411: #if defined(H2OPUS_USE_MPI)
412: distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
413: #endif
414: } else {
415: PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
416: hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
417: }
418: PetscCall(MatDenseRestoreArrayRead(B, (const PetscScalar **)&xx));
419: if (usesf) {
420: PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
421: PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
422: }
423: PetscCall(MatDenseRestoreArrayWrite(C, &yy));
424: #if PetscDefined(H2OPUS_USE_GPU)
425: } else {
426: PetscBool ciscuda, biscuda;
428: /* If not of type seqdensecuda, convert on the fly (i.e. allocate GPU memory) */
429: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &biscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
430: if (!biscuda) PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
431: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &ciscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
432: if (!ciscuda) {
433: C->assembled = PETSC_TRUE;
434: PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
435: }
436: PetscCall(MatDenseCUDAGetArrayRead(B, (const PetscScalar **)&xx));
437: PetscCall(MatDenseCUDAGetArrayWrite(C, &yy));
438: if (usesf) {
439: uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
440: uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
441: PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
442: PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
443: } else {
444: uxx = xx;
445: uyy = yy;
446: }
447: PetscCall(PetscLogGpuTimeBegin());
448: if (size > 1) {
449: PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
450: PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
451: #if defined(H2OPUS_USE_MPI)
452: distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
453: #endif
454: } else {
455: PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
456: hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
457: }
458: PetscCall(PetscLogGpuTimeEnd());
459: PetscCall(MatDenseCUDARestoreArrayRead(B, (const PetscScalar **)&xx));
460: if (usesf) {
461: PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
462: PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
463: }
464: PetscCall(MatDenseCUDARestoreArrayWrite(C, &yy));
465: if (!biscuda) PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
466: if (!ciscuda) PetscCall(MatConvert(C, MATDENSE, MAT_INPLACE_MATRIX, &C));
467: #endif
468: }
469: { /* log flops */
470: double gops, time, perf, dev;
471: HLibProfile::getHgemvPerf(gops, time, perf, dev);
472: #if PetscDefined(H2OPUS_USE_GPU)
473: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
474: else PetscCall(PetscLogGpuFlops(1e9 * gops));
475: #else
476: PetscCall(PetscLogFlops(1e9 * gops));
477: #endif
478: }
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: static PetscErrorCode MatProductNumeric_H2OPUS(Mat C)
483: {
484: Mat_Product *product = C->product;
486: PetscFunctionBegin;
487: MatCheckProduct(C, 1);
488: switch (product->type) {
489: case MATPRODUCT_AB:
490: PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_FALSE, product->B, C));
491: break;
492: case MATPRODUCT_AtB:
493: PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_TRUE, product->B, C));
494: break;
495: default:
496: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
497: }
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: static PetscErrorCode MatProductSymbolic_H2OPUS(Mat C)
502: {
503: Mat_Product *product = C->product;
504: PetscBool cisdense;
505: Mat A, B;
507: PetscFunctionBegin;
508: MatCheckProduct(C, 1);
509: A = product->A;
510: B = product->B;
511: switch (product->type) {
512: case MATPRODUCT_AB:
513: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
514: PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
515: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
516: if (!cisdense) {
517: PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
518: PetscCall(MatSetVecType(C, product->B->defaultvectype));
519: }
520: PetscCall(MatSetUp(C));
521: break;
522: case MATPRODUCT_AtB:
523: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
524: PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
525: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
526: if (!cisdense) {
527: PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
528: PetscCall(MatSetVecType(C, product->B->defaultvectype));
529: }
530: PetscCall(MatSetUp(C));
531: break;
532: default:
533: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
534: }
535: C->ops->productsymbolic = NULL;
536: C->ops->productnumeric = MatProductNumeric_H2OPUS;
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode MatProductSetFromOptions_H2OPUS(Mat C)
541: {
542: PetscFunctionBegin;
543: MatCheckProduct(C, 1);
544: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_H2OPUS;
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode MatMultKernel_H2OPUS(Mat A, Vec x, PetscScalar sy, Vec y, PetscBool trans)
549: {
550: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
551: #if defined(H2OPUS_USE_MPI)
552: h2opusHandle_t handle = h2opus->handle->handle;
553: #else
554: h2opusHandle_t handle = h2opus->handle;
555: #endif
556: PetscBool boundtocpu = PETSC_TRUE;
557: PetscInt n;
558: PetscScalar *xx, *yy, *uxx, *uyy;
559: PetscMPIInt size;
560: PetscBool usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);
562: PetscFunctionBegin;
563: HLibProfile::clear();
564: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
565: #if PetscDefined(H2OPUS_USE_GPU)
566: boundtocpu = A->boundtocpu;
567: #endif
568: if (usesf) PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
569: else n = A->rmap->n;
570: if (boundtocpu) {
571: PetscCall(VecGetArrayRead(x, (const PetscScalar **)&xx));
572: if (sy == 0.0) {
573: PetscCall(VecGetArrayWrite(y, &yy));
574: } else {
575: PetscCall(VecGetArray(y, &yy));
576: }
577: if (usesf) {
578: uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
579: uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
581: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
582: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
583: if (sy != 0.0) {
584: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
585: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
586: }
587: } else {
588: uxx = xx;
589: uyy = yy;
590: }
591: if (size > 1) {
592: PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
593: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
594: #if defined(H2OPUS_USE_MPI)
595: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, n, sy, uyy, n, 1, h2opus->handle);
596: #endif
597: } else {
598: PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
599: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, n, sy, uyy, n, 1, handle);
600: }
601: PetscCall(VecRestoreArrayRead(x, (const PetscScalar **)&xx));
602: if (usesf) {
603: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
604: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
605: }
606: if (sy == 0.0) {
607: PetscCall(VecRestoreArrayWrite(y, &yy));
608: } else {
609: PetscCall(VecRestoreArray(y, &yy));
610: }
611: #if PetscDefined(H2OPUS_USE_GPU)
612: } else {
613: PetscCall(VecCUDAGetArrayRead(x, (const PetscScalar **)&xx));
614: if (sy == 0.0) {
615: PetscCall(VecCUDAGetArrayWrite(y, &yy));
616: } else {
617: PetscCall(VecCUDAGetArray(y, &yy));
618: }
619: if (usesf) {
620: uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
621: uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
623: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
624: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
625: if (sy != 0.0) {
626: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
627: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
628: }
629: } else {
630: uxx = xx;
631: uyy = yy;
632: }
633: PetscCall(PetscLogGpuTimeBegin());
634: if (size > 1) {
635: PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
636: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
637: #if defined(H2OPUS_USE_MPI)
638: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, n, sy, uyy, n, 1, h2opus->handle);
639: #endif
640: } else {
641: PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
642: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, n, sy, uyy, n, 1, handle);
643: }
644: PetscCall(PetscLogGpuTimeEnd());
645: PetscCall(VecCUDARestoreArrayRead(x, (const PetscScalar **)&xx));
646: if (usesf) {
647: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
648: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
649: }
650: if (sy == 0.0) {
651: PetscCall(VecCUDARestoreArrayWrite(y, &yy));
652: } else {
653: PetscCall(VecCUDARestoreArray(y, &yy));
654: }
655: #endif
656: }
657: { /* log flops */
658: double gops, time, perf, dev;
659: HLibProfile::getHgemvPerf(gops, time, perf, dev);
660: #if PetscDefined(H2OPUS_USE_GPU)
661: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
662: else PetscCall(PetscLogGpuFlops(1e9 * gops));
663: #else
664: PetscCall(PetscLogFlops(1e9 * gops));
665: #endif
666: }
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: static PetscErrorCode MatMultTranspose_H2OPUS(Mat A, Vec x, Vec y)
671: {
672: PetscBool xiscuda, yiscuda;
674: PetscFunctionBegin;
675: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
676: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
677: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
678: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_TRUE));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode MatMult_H2OPUS(Mat A, Vec x, Vec y)
683: {
684: PetscBool xiscuda, yiscuda;
686: PetscFunctionBegin;
687: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
688: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
689: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
690: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_FALSE));
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: static PetscErrorCode MatMultTransposeAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
695: {
696: PetscBool xiscuda, ziscuda;
698: PetscFunctionBegin;
699: PetscCall(VecCopy(y, z));
700: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
701: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
702: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
703: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_TRUE));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: static PetscErrorCode MatMultAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
708: {
709: PetscBool xiscuda, ziscuda;
711: PetscFunctionBegin;
712: PetscCall(VecCopy(y, z));
713: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
714: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
715: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
716: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_FALSE));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode MatScale_H2OPUS(Mat A, PetscScalar s)
721: {
722: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
724: PetscFunctionBegin;
725: a->s *= s;
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: static PetscErrorCode MatSetFromOptions_H2OPUS(Mat A, PetscOptionItems PetscOptionsObject)
730: {
731: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
733: PetscFunctionBegin;
734: PetscOptionsHeadBegin(PetscOptionsObject, "H2OPUS options");
735: PetscCall(PetscOptionsInt("-mat_h2opus_leafsize", "Leaf size of cluster tree", NULL, a->leafsize, &a->leafsize, NULL));
736: PetscCall(PetscOptionsReal("-mat_h2opus_eta", "Admissibility condition tolerance", NULL, a->eta, &a->eta, NULL));
737: PetscCall(PetscOptionsInt("-mat_h2opus_order", "Basis order for off-diagonal sampling when constructed from kernel", NULL, a->basisord, &a->basisord, NULL));
738: PetscCall(PetscOptionsInt("-mat_h2opus_maxrank", "Maximum rank when constructed from matvecs", NULL, a->max_rank, &a->max_rank, NULL));
739: PetscCall(PetscOptionsInt("-mat_h2opus_samples", "Maximum number of samples to be taken concurrently when constructing from matvecs", NULL, a->bs, &a->bs, NULL));
740: PetscCall(PetscOptionsInt("-mat_h2opus_normsamples", "Maximum number of samples to be when estimating norms", NULL, a->norm_max_samples, &a->norm_max_samples, NULL));
741: PetscCall(PetscOptionsReal("-mat_h2opus_rtol", "Relative tolerance for construction from sampling", NULL, a->rtol, &a->rtol, NULL));
742: PetscCall(PetscOptionsBool("-mat_h2opus_check", "Check error when constructing from sampling during MatAssemblyEnd()", NULL, a->check_construction, &a->check_construction, NULL));
743: PetscCall(PetscOptionsBool("-mat_h2opus_hara_verbose", "Verbose output from hara construction", NULL, a->hara_verbose, &a->hara_verbose, NULL));
744: PetscCall(PetscOptionsBool("-mat_h2opus_resize", "Resize after compression", NULL, a->resize, &a->resize, NULL));
745: PetscOptionsHeadEnd();
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat, PetscInt, const PetscReal[], PetscBool, MatH2OpusKernelFn *, void *);
751: static PetscErrorCode MatH2OpusInferCoordinates_Private(Mat A)
752: {
753: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
754: Vec c;
755: PetscInt spacedim;
756: const PetscScalar *coords;
758: PetscFunctionBegin;
759: if (a->ptcloud) PetscFunctionReturn(PETSC_SUCCESS);
760: PetscCall(PetscObjectQuery((PetscObject)A, "__math2opus_coords", (PetscObject *)&c));
761: if (!c && a->sampler) {
762: Mat S = a->sampler->GetSamplingMat();
764: PetscCall(PetscObjectQuery((PetscObject)S, "__math2opus_coords", (PetscObject *)&c));
765: }
766: if (!c) {
767: PetscCall(MatH2OpusSetCoords_H2OPUS(A, -1, NULL, PETSC_FALSE, NULL, NULL));
768: } else {
769: PetscCall(VecGetArrayRead(c, &coords));
770: PetscCall(VecGetBlockSize(c, &spacedim));
771: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, PETSC_FALSE, NULL, NULL));
772: PetscCall(VecRestoreArrayRead(c, &coords));
773: }
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: static PetscErrorCode MatSetUpMultiply_H2OPUS(Mat A)
778: {
779: MPI_Comm comm;
780: PetscMPIInt size;
781: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
782: PetscInt n = 0, *idx = NULL;
783: int *iidx = NULL;
784: PetscCopyMode own;
785: PetscBool rid;
787: PetscFunctionBegin;
788: if (a->multsetup) PetscFunctionReturn(PETSC_SUCCESS);
789: if (a->sf) { /* MatDuplicate_H2OPUS takes reference to the SF */
790: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
791: #if PetscDefined(H2OPUS_USE_GPU)
792: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
793: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
794: a->xxs_gpu = 1;
795: a->yys_gpu = 1;
796: #endif
797: a->xx = new thrust::host_vector<PetscScalar>(n);
798: a->yy = new thrust::host_vector<PetscScalar>(n);
799: a->xxs = 1;
800: a->yys = 1;
801: } else {
802: IS is;
803: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
804: PetscCallMPI(MPI_Comm_size(comm, &size));
805: if (!a->h2opus_indexmap) {
806: if (size > 1) {
807: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
808: #if defined(H2OPUS_USE_MPI)
809: iidx = MatH2OpusGetThrustPointer(a->dist_hmatrix->basis_tree.basis_branch.index_map);
810: n = a->dist_hmatrix->basis_tree.basis_branch.index_map.size();
811: #endif
812: } else {
813: iidx = MatH2OpusGetThrustPointer(a->hmatrix->u_basis_tree.index_map);
814: n = a->hmatrix->u_basis_tree.index_map.size();
815: }
817: if (PetscDefined(USE_64BIT_INDICES)) {
818: PetscInt i;
820: own = PETSC_OWN_POINTER;
821: PetscCall(PetscMalloc1(n, &idx));
822: for (i = 0; i < n; i++) idx[i] = iidx[i];
823: } else {
824: own = PETSC_COPY_VALUES;
825: idx = (PetscInt *)iidx;
826: }
827: PetscCall(ISCreateGeneral(comm, n, idx, own, &is));
828: PetscCall(ISSetPermutation(is));
829: PetscCall(ISViewFromOptions(is, (PetscObject)A, "-mat_h2opus_indexmap_view"));
830: a->h2opus_indexmap = is;
831: }
832: PetscCall(ISGetLocalSize(a->h2opus_indexmap, &n));
833: PetscCall(ISGetIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
834: rid = (PetscBool)(n == A->rmap->n);
835: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &rid, 1, MPI_C_BOOL, MPI_LAND, comm));
836: if (rid) PetscCall(ISIdentity(a->h2opus_indexmap, &rid));
837: if (!rid) {
838: if (size > 1) { /* Parallel distribution may be different, save it here for fast path in MatMult (see MatH2OpusSetNativeMult) */
839: PetscCall(PetscLayoutCreate(comm, &a->h2opus_rmap));
840: PetscCall(PetscLayoutSetLocalSize(a->h2opus_rmap, n));
841: PetscCall(PetscLayoutSetUp(a->h2opus_rmap));
842: PetscCall(PetscLayoutReference(a->h2opus_rmap, &a->h2opus_cmap));
843: }
844: PetscCall(PetscSFCreate(comm, &a->sf));
845: PetscCall(PetscSFSetGraphLayout(a->sf, A->rmap, n, NULL, PETSC_OWN_POINTER, idx));
846: PetscCall(PetscSFSetUp(a->sf));
847: PetscCall(PetscSFViewFromOptions(a->sf, (PetscObject)A, "-mat_h2opus_sf_view"));
848: #if PetscDefined(H2OPUS_USE_GPU)
849: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
850: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
851: a->xxs_gpu = 1;
852: a->yys_gpu = 1;
853: #endif
854: a->xx = new thrust::host_vector<PetscScalar>(n);
855: a->yy = new thrust::host_vector<PetscScalar>(n);
856: a->xxs = 1;
857: a->yys = 1;
858: }
859: PetscCall(ISRestoreIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
860: }
861: a->multsetup = PETSC_TRUE;
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
865: static PetscErrorCode MatAssemblyEnd_H2OPUS(Mat A, MatAssemblyType assemblytype)
866: {
867: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
868: #if defined(H2OPUS_USE_MPI)
869: h2opusHandle_t handle = a->handle->handle;
870: #else
871: h2opusHandle_t handle = a->handle;
872: #endif
873: PetscBool kernel = PETSC_FALSE;
874: PetscBool boundtocpu = PETSC_TRUE;
875: PetscBool samplingdone = PETSC_FALSE;
876: MPI_Comm comm;
877: PetscMPIInt size;
879: PetscFunctionBegin;
880: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
881: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
882: PetscCheck(A->rmap->N == A->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
884: /* XXX */
885: a->leafsize = PetscMin(a->leafsize, PetscMin(A->rmap->N, A->cmap->N));
887: PetscCallMPI(MPI_Comm_size(comm, &size));
888: /* TODO REUSABILITY of geometric construction */
889: delete a->hmatrix;
890: delete a->dist_hmatrix;
891: #if PetscDefined(H2OPUS_USE_GPU)
892: delete a->hmatrix_gpu;
893: delete a->dist_hmatrix_gpu;
894: #endif
895: a->orthogonal = PETSC_FALSE;
897: /* TODO: other? */
898: H2OpusBoxCenterAdmissibility adm(a->eta);
900: PetscCall(PetscLogEventBegin(MAT_H2Opus_Build, A, 0, 0, 0));
901: if (size > 1) {
902: #if defined(H2OPUS_USE_MPI)
903: a->dist_hmatrix = new DistributedHMatrix(A->rmap->n /* ,A->symmetric */);
904: #else
905: a->dist_hmatrix = NULL;
906: #endif
907: } else a->hmatrix = new HMatrix(A->rmap->n, A->symmetric == PETSC_BOOL3_TRUE);
908: PetscCall(MatH2OpusInferCoordinates_Private(A));
909: PetscCheck(a->ptcloud, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing pointcloud");
910: if (a->kernel) {
911: BoxEntryGen<PetscScalar, H2OPUS_HWTYPE_CPU, PetscFunctionGenerator<PetscScalar>> entry_gen(*a->kernel);
912: if (size > 1) {
913: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
914: #if defined(H2OPUS_USE_MPI)
915: buildDistributedHMatrix(*a->dist_hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord, a->handle);
916: #endif
917: } else {
918: buildHMatrix(*a->hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord);
919: }
920: kernel = PETSC_TRUE;
921: } else {
922: PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Construction from sampling not supported in parallel");
923: buildHMatrixStructure(*a->hmatrix, a->ptcloud, a->leafsize, adm);
924: }
925: PetscCall(MatSetUpMultiply_H2OPUS(A));
927: #if PetscDefined(H2OPUS_USE_GPU)
928: boundtocpu = A->boundtocpu;
929: if (!boundtocpu) {
930: if (size > 1) {
931: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
932: #if defined(H2OPUS_USE_MPI)
933: a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
934: #endif
935: } else {
936: a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
937: }
938: }
939: #endif
940: if (size == 1) {
941: if (!kernel && a->sampler && a->sampler->GetSamplingMat()) {
942: PetscReal Anorm;
943: bool verbose;
945: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_hara_verbose", &a->hara_verbose, NULL));
946: verbose = a->hara_verbose;
947: PetscCall(MatNormApproximate(a->sampler->GetSamplingMat(), NORM_2, a->norm_max_samples, &Anorm));
948: 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));
949: if (a->sf && !a->nativemult) a->sampler->SetIndexMap(a->hmatrix->u_basis_tree.index_map.size(), a->hmatrix->u_basis_tree.index_map.data());
950: a->sampler->SetStream(handle->getMainStream());
951: if (boundtocpu) {
952: a->sampler->SetGPUSampling(false);
953: hara(a->sampler, *a->hmatrix, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
954: #if PetscDefined(H2OPUS_USE_GPU)
955: } else {
956: a->sampler->SetGPUSampling(true);
957: hara(a->sampler, *a->hmatrix_gpu, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
958: #endif
959: }
960: samplingdone = PETSC_TRUE;
961: }
962: }
963: #if PetscDefined(H2OPUS_USE_GPU)
964: if (!boundtocpu) {
965: delete a->hmatrix;
966: delete a->dist_hmatrix;
967: a->hmatrix = NULL;
968: a->dist_hmatrix = NULL;
969: }
970: A->offloadmask = boundtocpu ? PETSC_OFFLOAD_CPU : PETSC_OFFLOAD_GPU;
971: #endif
972: PetscCall(PetscLogEventEnd(MAT_H2Opus_Build, A, 0, 0, 0));
974: if (!a->s) a->s = 1.0;
975: A->assembled = PETSC_TRUE;
977: if (samplingdone) {
978: PetscBool check = a->check_construction;
979: PetscBool checke = PETSC_FALSE;
981: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check", &check, NULL));
982: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check_explicit", &checke, NULL));
983: if (check) {
984: Mat E, Ae;
985: PetscReal n1, ni, n2;
986: PetscReal n1A, niA, n2A;
987: PetscErrorCodeFn *normfunc;
989: Ae = a->sampler->GetSamplingMat();
990: PetscCall(MatConvert(A, MATSHELL, MAT_INITIAL_MATRIX, &E));
991: PetscCall(MatShellSetOperation(E, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
992: PetscCall(MatAXPY(E, -1.0, Ae, DIFFERENT_NONZERO_PATTERN));
993: PetscCall(MatNorm(E, NORM_1, &n1));
994: PetscCall(MatNorm(E, NORM_INFINITY, &ni));
995: PetscCall(MatNorm(E, NORM_2, &n2));
996: if (checke) {
997: Mat eA, eE, eAe;
999: PetscCall(MatComputeOperator(A, MATAIJ, &eA));
1000: PetscCall(MatComputeOperator(E, MATAIJ, &eE));
1001: PetscCall(MatComputeOperator(Ae, MATAIJ, &eAe));
1002: PetscCall(MatFilter(eA, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
1003: PetscCall(MatFilter(eE, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
1004: PetscCall(MatFilter(eAe, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
1005: PetscCall(PetscObjectSetName((PetscObject)eA, "H2Mat"));
1006: PetscCall(MatView(eA, NULL));
1007: PetscCall(PetscObjectSetName((PetscObject)eAe, "S"));
1008: PetscCall(MatView(eAe, NULL));
1009: PetscCall(PetscObjectSetName((PetscObject)eE, "H2Mat - S"));
1010: PetscCall(MatView(eE, NULL));
1011: PetscCall(MatDestroy(&eA));
1012: PetscCall(MatDestroy(&eE));
1013: PetscCall(MatDestroy(&eAe));
1014: }
1016: PetscCall(MatGetOperation(Ae, MATOP_NORM, &normfunc));
1017: PetscCall(MatSetOperation(Ae, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
1018: PetscCall(MatNorm(Ae, NORM_1, &n1A));
1019: PetscCall(MatNorm(Ae, NORM_INFINITY, &niA));
1020: PetscCall(MatNorm(Ae, NORM_2, &n2A));
1021: n1A = PetscMax(n1A, PETSC_SMALL);
1022: n2A = PetscMax(n2A, PETSC_SMALL);
1023: niA = PetscMax(niA, PETSC_SMALL);
1024: PetscCall(MatSetOperation(Ae, MATOP_NORM, normfunc));
1025: 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)));
1026: PetscCall(MatDestroy(&E));
1027: }
1028: a->sampler->SetSamplingMat(NULL);
1029: }
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: static PetscErrorCode MatZeroEntries_H2OPUS(Mat A)
1034: {
1035: PetscMPIInt size;
1036: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1038: PetscFunctionBegin;
1039: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1040: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet supported");
1041: a->hmatrix->clearData();
1042: #if PetscDefined(H2OPUS_USE_GPU)
1043: if (a->hmatrix_gpu) a->hmatrix_gpu->clearData();
1044: #endif
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: static PetscErrorCode MatDuplicate_H2OPUS(Mat B, MatDuplicateOption op, Mat *nA)
1049: {
1050: Mat A;
1051: Mat_H2OPUS *a, *b = (Mat_H2OPUS *)B->data;
1052: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1053: MPI_Comm comm;
1055: PetscFunctionBegin;
1056: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1057: PetscCall(MatCreate(comm, &A));
1058: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1059: PetscCall(MatSetType(A, MATH2OPUS));
1060: PetscCall(MatPropagateSymmetryOptions(B, A));
1061: a = (Mat_H2OPUS *)A->data;
1063: a->eta = b->eta;
1064: a->leafsize = b->leafsize;
1065: a->basisord = b->basisord;
1066: a->max_rank = b->max_rank;
1067: a->bs = b->bs;
1068: a->rtol = b->rtol;
1069: a->norm_max_samples = b->norm_max_samples;
1070: if (op == MAT_COPY_VALUES) a->s = b->s;
1072: a->ptcloud = new PetscPointCloud<PetscReal>(*b->ptcloud);
1073: if (op == MAT_COPY_VALUES && b->kernel) a->kernel = new PetscFunctionGenerator<PetscScalar>(*b->kernel);
1075: #if defined(H2OPUS_USE_MPI)
1076: if (b->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*b->dist_hmatrix);
1077: #if PetscDefined(H2OPUS_USE_GPU)
1078: if (b->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*b->dist_hmatrix_gpu);
1079: #endif
1080: #endif
1081: if (b->hmatrix) {
1082: a->hmatrix = new HMatrix(*b->hmatrix);
1083: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix->clearData();
1084: }
1085: #if PetscDefined(H2OPUS_USE_GPU)
1086: if (b->hmatrix_gpu) {
1087: a->hmatrix_gpu = new HMatrix_GPU(*b->hmatrix_gpu);
1088: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix_gpu->clearData();
1089: }
1090: #endif
1091: if (b->sf) {
1092: PetscCall(PetscObjectReference((PetscObject)b->sf));
1093: a->sf = b->sf;
1094: }
1095: if (b->h2opus_indexmap) {
1096: PetscCall(PetscObjectReference((PetscObject)b->h2opus_indexmap));
1097: a->h2opus_indexmap = b->h2opus_indexmap;
1098: }
1100: PetscCall(MatSetUp(A));
1101: PetscCall(MatSetUpMultiply_H2OPUS(A));
1102: if (op == MAT_COPY_VALUES) {
1103: A->assembled = PETSC_TRUE;
1104: a->orthogonal = b->orthogonal;
1105: #if PetscDefined(H2OPUS_USE_GPU)
1106: A->offloadmask = B->offloadmask;
1107: #endif
1108: }
1109: #if PetscDefined(H2OPUS_USE_GPU)
1110: iscpu = B->boundtocpu;
1111: #endif
1112: PetscCall(MatBindToCPU(A, iscpu));
1114: *nA = A;
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: static PetscErrorCode MatView_H2OPUS(Mat A, PetscViewer view)
1119: {
1120: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1121: PetscBool isascii, vieweps;
1122: PetscMPIInt size;
1123: PetscViewerFormat format;
1125: PetscFunctionBegin;
1126: PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
1127: PetscCall(PetscViewerGetFormat(view, &format));
1128: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1129: if (isascii) {
1130: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1131: if (size == 1) {
1132: FILE *fp;
1133: PetscCall(PetscViewerASCIIGetPointer(view, &fp));
1134: dumpHMatrix(*h2opus->hmatrix, 6, fp);
1135: }
1136: } else {
1137: PetscCall(PetscViewerASCIIPrintf(view, " H-Matrix constructed from %s\n", h2opus->kernel ? "Kernel" : "Mat"));
1138: PetscCall(PetscViewerASCIIPrintf(view, " PointCloud dim %" PetscInt_FMT "\n", h2opus->ptcloud ? h2opus->ptcloud->getDimension() : 0));
1139: PetscCall(PetscViewerASCIIPrintf(view, " Admissibility parameters: leaf size %" PetscInt_FMT ", eta %g\n", h2opus->leafsize, (double)h2opus->eta));
1140: if (!h2opus->kernel) {
1141: PetscCall(PetscViewerASCIIPrintf(view, " Sampling parameters: max_rank %" PetscInt_FMT ", samples %" PetscInt_FMT ", tolerance %g\n", h2opus->max_rank, h2opus->bs, (double)h2opus->rtol));
1142: } else {
1143: PetscCall(PetscViewerASCIIPrintf(view, " Off-diagonal blocks approximation order %" PetscInt_FMT "\n", h2opus->basisord));
1144: }
1145: PetscCall(PetscViewerASCIIPrintf(view, " Number of samples for norms %" PetscInt_FMT "\n", h2opus->norm_max_samples));
1146: if (size == 1) {
1147: double dense_mem_cpu = h2opus->hmatrix ? h2opus->hmatrix->getDenseMemoryUsage() : 0;
1148: double low_rank_cpu = h2opus->hmatrix ? h2opus->hmatrix->getLowRankMemoryUsage() : 0;
1149: #if PetscDefined(HAVE_CUDA)
1150: double dense_mem_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getDenseMemoryUsage() : 0;
1151: double low_rank_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getLowRankMemoryUsage() : 0;
1152: #endif
1153: 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));
1154: #if PetscDefined(HAVE_CUDA)
1155: 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));
1156: #endif
1157: } else {
1158: #if PetscDefined(HAVE_CUDA)
1159: double matrix_mem[4] = {0., 0., 0., 0.};
1160: PetscMPIInt rsize = 4;
1161: #else
1162: double matrix_mem[2] = {0., 0.};
1163: PetscMPIInt rsize = 2;
1164: #endif
1165: #if defined(H2OPUS_USE_MPI)
1166: matrix_mem[0] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalDenseMemoryUsage() : 0;
1167: matrix_mem[1] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalLowRankMemoryUsage() : 0;
1168: #if PetscDefined(HAVE_CUDA)
1169: matrix_mem[2] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalDenseMemoryUsage() : 0;
1170: matrix_mem[3] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalLowRankMemoryUsage() : 0;
1171: #endif
1172: #endif
1173: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, matrix_mem, rsize, MPI_DOUBLE_PRECISION, MPI_SUM, PetscObjectComm((PetscObject)A)));
1174: 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]));
1175: #if PetscDefined(HAVE_CUDA)
1176: 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]));
1177: #endif
1178: }
1179: }
1180: }
1181: vieweps = PETSC_FALSE;
1182: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps", &vieweps, NULL));
1183: if (vieweps) {
1184: char filename[256];
1185: const char *name;
1187: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1188: PetscCall(PetscSNPrintf(filename, sizeof(filename), "%s_structure.eps", name));
1189: PetscCall(PetscOptionsGetString(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps_filename", filename, sizeof(filename), NULL));
1190: outputEps(*h2opus->hmatrix, filename);
1191: }
1192: PetscFunctionReturn(PETSC_SUCCESS);
1193: }
1195: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat A, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, void *kernelctx)
1196: {
1197: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1198: PetscReal *gcoords;
1199: PetscInt N;
1200: MPI_Comm comm;
1201: PetscMPIInt size;
1202: PetscBool cong;
1204: PetscFunctionBegin;
1205: PetscCall(PetscLayoutSetUp(A->rmap));
1206: PetscCall(PetscLayoutSetUp(A->cmap));
1207: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1208: PetscCall(MatHasCongruentLayouts(A, &cong));
1209: PetscCheck(cong, comm, PETSC_ERR_SUP, "Only for square matrices with congruent layouts");
1210: N = A->rmap->N;
1211: PetscCallMPI(MPI_Comm_size(comm, &size));
1212: if (spacedim > 0 && size > 1 && cdist) {
1213: PetscSF sf;
1214: MPI_Datatype dtype;
1216: PetscCallMPI(MPI_Type_contiguous(spacedim, MPIU_REAL, &dtype));
1217: PetscCallMPI(MPI_Type_commit(&dtype));
1219: PetscCall(PetscSFCreate(comm, &sf));
1220: PetscCall(PetscSFSetGraphWithPattern(sf, A->rmap, PETSCSF_PATTERN_ALLGATHER));
1221: PetscCall(PetscMalloc1(spacedim * N, &gcoords));
1222: PetscCall(PetscSFBcastBegin(sf, dtype, coords, gcoords, MPI_REPLACE));
1223: PetscCall(PetscSFBcastEnd(sf, dtype, coords, gcoords, MPI_REPLACE));
1224: PetscCall(PetscSFDestroy(&sf));
1225: PetscCallMPI(MPI_Type_free(&dtype));
1226: } else gcoords = (PetscReal *)coords;
1228: delete h2opus->ptcloud;
1229: delete h2opus->kernel;
1230: h2opus->ptcloud = new PetscPointCloud<PetscReal>(spacedim, N, gcoords);
1231: if (kernel) h2opus->kernel = new PetscFunctionGenerator<PetscScalar>(kernel, spacedim, kernelctx);
1232: if (gcoords != coords) PetscCall(PetscFree(gcoords));
1233: A->preallocated = PETSC_TRUE;
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: #if PetscDefined(H2OPUS_USE_GPU)
1238: static PetscErrorCode MatBindToCPU_H2OPUS(Mat A, PetscBool flg)
1239: {
1240: PetscMPIInt size;
1241: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1243: PetscFunctionBegin;
1244: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1245: if (flg && A->offloadmask == PETSC_OFFLOAD_GPU) {
1246: if (size > 1) {
1247: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1248: #if defined(H2OPUS_USE_MPI)
1249: if (!a->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix_gpu);
1250: else *a->dist_hmatrix = *a->dist_hmatrix_gpu;
1251: #endif
1252: } else {
1253: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1254: if (!a->hmatrix) a->hmatrix = new HMatrix(*a->hmatrix_gpu);
1255: else *a->hmatrix = *a->hmatrix_gpu;
1256: }
1257: delete a->hmatrix_gpu;
1258: delete a->dist_hmatrix_gpu;
1259: a->hmatrix_gpu = NULL;
1260: a->dist_hmatrix_gpu = NULL;
1261: A->offloadmask = PETSC_OFFLOAD_CPU;
1262: } else if (!flg && A->offloadmask == PETSC_OFFLOAD_CPU) {
1263: if (size > 1) {
1264: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1265: #if defined(H2OPUS_USE_MPI)
1266: if (!a->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
1267: else *a->dist_hmatrix_gpu = *a->dist_hmatrix;
1268: #endif
1269: } else {
1270: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1271: if (!a->hmatrix_gpu) a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
1272: else *a->hmatrix_gpu = *a->hmatrix;
1273: }
1274: delete a->hmatrix;
1275: delete a->dist_hmatrix;
1276: a->hmatrix = NULL;
1277: a->dist_hmatrix = NULL;
1278: A->offloadmask = PETSC_OFFLOAD_GPU;
1279: }
1280: PetscCall(PetscFree(A->defaultvectype));
1281: if (!flg) {
1282: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1283: } else {
1284: PetscCall(PetscStrallocpy(VECSTANDARD, &A->defaultvectype));
1285: }
1286: A->boundtocpu = flg;
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1289: #endif
1291: /*MC
1292: MATH2OPUS = "h2opus" - A matrix type for hierarchical matrices using the H2Opus package {cite}`zampinibouakaramturkiyyahkniokeyes2022`.
1294: Options Database Key:
1295: . -mat_type h2opus - matrix type to "h2opus"
1297: Level: beginner
1299: Notes:
1300: H2Opus implements hierarchical matrices in the $H^2$ flavor. It supports CPU or NVIDIA GPUs.
1302: For CPU only builds, use `./configure --download-h2opus --download-thrust` to install PETSc to use H2Opus.
1303: In order to run on NVIDIA GPUs, use `./configure --download-h2opus --download-magma --download-kblas`.
1305: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATHTOOL`, `MATDENSE`, `MatCreateH2OpusFromKernel()`, `MatCreateH2OpusFromMat()`
1306: M*/
1307: PETSC_EXTERN PetscErrorCode MatCreate_H2OPUS(Mat A)
1308: {
1309: Mat_H2OPUS *a;
1310: PetscMPIInt size;
1312: PetscFunctionBegin;
1313: #if PetscDefined(H2OPUS_USE_GPU)
1314: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
1315: #endif
1316: PetscCall(PetscNew(&a));
1317: A->data = (void *)a;
1319: a->eta = 0.9;
1320: a->leafsize = 32;
1321: a->basisord = 4;
1322: a->max_rank = 64;
1323: a->bs = 32;
1324: a->rtol = 1.e-4;
1325: a->s = 1.0;
1326: a->norm_max_samples = 10;
1327: a->resize = PETSC_TRUE; /* reallocate after compression */
1328: #if defined(H2OPUS_USE_MPI)
1329: h2opusCreateDistributedHandleComm(&a->handle, PetscObjectComm((PetscObject)A));
1330: #else
1331: h2opusCreateHandle(&a->handle);
1332: #endif
1333: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1334: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATH2OPUS));
1335: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
1337: A->ops->destroy = MatDestroy_H2OPUS;
1338: A->ops->view = MatView_H2OPUS;
1339: A->ops->assemblyend = MatAssemblyEnd_H2OPUS;
1340: A->ops->mult = MatMult_H2OPUS;
1341: A->ops->multtranspose = MatMultTranspose_H2OPUS;
1342: A->ops->multadd = MatMultAdd_H2OPUS;
1343: A->ops->multtransposeadd = MatMultTransposeAdd_H2OPUS;
1344: A->ops->scale = MatScale_H2OPUS;
1345: A->ops->duplicate = MatDuplicate_H2OPUS;
1346: A->ops->setfromoptions = MatSetFromOptions_H2OPUS;
1347: A->ops->norm = MatNorm_H2OPUS;
1348: A->ops->zeroentries = MatZeroEntries_H2OPUS;
1349: #if PetscDefined(H2OPUS_USE_GPU)
1350: A->ops->bindtocpu = MatBindToCPU_H2OPUS;
1351: #endif
1353: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", MatProductSetFromOptions_H2OPUS));
1354: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", MatProductSetFromOptions_H2OPUS));
1355: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", MatProductSetFromOptions_H2OPUS));
1356: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", MatProductSetFromOptions_H2OPUS));
1357: #if PetscDefined(H2OPUS_USE_GPU)
1358: PetscCall(PetscFree(A->defaultvectype));
1359: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1360: #endif
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@
1365: MatH2OpusOrthogonalize - Orthogonalize the basis tree of a hierarchical matrix.
1367: Input Parameter:
1368: . A - the matrix
1370: Level: intermediate
1372: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`
1373: @*/
1374: PetscErrorCode MatH2OpusOrthogonalize(Mat A)
1375: {
1376: PetscBool ish2opus;
1377: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1378: PetscMPIInt size;
1379: PetscBool boundtocpu = PETSC_TRUE;
1381: PetscFunctionBegin;
1384: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1385: if (!ish2opus) PetscFunctionReturn(PETSC_SUCCESS);
1386: if (a->orthogonal) PetscFunctionReturn(PETSC_SUCCESS);
1387: HLibProfile::clear();
1388: PetscCall(PetscLogEventBegin(MAT_H2Opus_Orthog, A, 0, 0, 0));
1389: #if PetscDefined(H2OPUS_USE_GPU)
1390: boundtocpu = A->boundtocpu;
1391: #endif
1392: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1393: if (size > 1) {
1394: if (boundtocpu) {
1395: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1396: #if defined(H2OPUS_USE_MPI)
1397: distributed_horthog(*a->dist_hmatrix, a->handle);
1398: #endif
1399: #if PetscDefined(H2OPUS_USE_GPU)
1400: A->offloadmask = PETSC_OFFLOAD_CPU;
1401: } else {
1402: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1403: PetscCall(PetscLogGpuTimeBegin());
1404: #if defined(H2OPUS_USE_MPI)
1405: distributed_horthog(*a->dist_hmatrix_gpu, a->handle);
1406: #endif
1407: PetscCall(PetscLogGpuTimeEnd());
1408: #endif
1409: }
1410: } else {
1411: #if defined(H2OPUS_USE_MPI)
1412: h2opusHandle_t handle = a->handle->handle;
1413: #else
1414: h2opusHandle_t handle = a->handle;
1415: #endif
1416: if (boundtocpu) {
1417: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1418: horthog(*a->hmatrix, handle);
1419: #if PetscDefined(H2OPUS_USE_GPU)
1420: A->offloadmask = PETSC_OFFLOAD_CPU;
1421: } else {
1422: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1423: PetscCall(PetscLogGpuTimeBegin());
1424: horthog(*a->hmatrix_gpu, handle);
1425: PetscCall(PetscLogGpuTimeEnd());
1426: #endif
1427: }
1428: }
1429: a->orthogonal = PETSC_TRUE;
1430: { /* log flops */
1431: double gops, time, perf, dev;
1432: HLibProfile::getHorthogPerf(gops, time, perf, dev);
1433: #if PetscDefined(H2OPUS_USE_GPU)
1434: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1435: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1436: #else
1437: PetscCall(PetscLogFlops(1e9 * gops));
1438: #endif
1439: }
1440: PetscCall(PetscLogEventEnd(MAT_H2Opus_Orthog, A, 0, 0, 0));
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /*@
1445: MatH2OpusCompress - Compress a hierarchical matrix.
1447: Input Parameters:
1448: + A - the matrix
1449: - tol - the absolute truncation threshold
1451: Level: intermediate
1453: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusOrthogonalize()`
1454: @*/
1455: PetscErrorCode MatH2OpusCompress(Mat A, PetscReal tol)
1456: {
1457: PetscBool ish2opus;
1458: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1459: PetscMPIInt size;
1460: PetscBool boundtocpu = PETSC_TRUE;
1462: PetscFunctionBegin;
1466: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1467: if (!ish2opus || tol <= 0.0) PetscFunctionReturn(PETSC_SUCCESS);
1468: PetscCall(MatH2OpusOrthogonalize(A));
1469: HLibProfile::clear();
1470: PetscCall(PetscLogEventBegin(MAT_H2Opus_Compress, A, 0, 0, 0));
1471: #if PetscDefined(H2OPUS_USE_GPU)
1472: boundtocpu = A->boundtocpu;
1473: #endif
1474: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1475: if (size > 1) {
1476: if (boundtocpu) {
1477: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1478: #if defined(H2OPUS_USE_MPI)
1479: distributed_hcompress(*a->dist_hmatrix, tol, a->handle);
1480: if (a->resize) {
1481: DistributedHMatrix *dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix);
1482: delete a->dist_hmatrix;
1483: a->dist_hmatrix = dist_hmatrix;
1484: }
1485: #endif
1486: #if PetscDefined(H2OPUS_USE_GPU)
1487: A->offloadmask = PETSC_OFFLOAD_CPU;
1488: } else {
1489: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1490: PetscCall(PetscLogGpuTimeBegin());
1491: #if defined(H2OPUS_USE_MPI)
1492: distributed_hcompress(*a->dist_hmatrix_gpu, tol, a->handle);
1494: if (a->resize) {
1495: DistributedHMatrix_GPU *dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix_gpu);
1496: delete a->dist_hmatrix_gpu;
1497: a->dist_hmatrix_gpu = dist_hmatrix_gpu;
1498: }
1499: #endif
1500: PetscCall(PetscLogGpuTimeEnd());
1501: #endif
1502: }
1503: } else {
1504: #if defined(H2OPUS_USE_MPI)
1505: h2opusHandle_t handle = a->handle->handle;
1506: #else
1507: h2opusHandle_t handle = a->handle;
1508: #endif
1509: if (boundtocpu) {
1510: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1511: hcompress(*a->hmatrix, tol, handle);
1513: if (a->resize) {
1514: HMatrix *hmatrix = new HMatrix(*a->hmatrix);
1515: delete a->hmatrix;
1516: a->hmatrix = hmatrix;
1517: }
1518: #if PetscDefined(H2OPUS_USE_GPU)
1519: A->offloadmask = PETSC_OFFLOAD_CPU;
1520: } else {
1521: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1522: PetscCall(PetscLogGpuTimeBegin());
1523: hcompress(*a->hmatrix_gpu, tol, handle);
1524: PetscCall(PetscLogGpuTimeEnd());
1526: if (a->resize) {
1527: HMatrix_GPU *hmatrix_gpu = new HMatrix_GPU(*a->hmatrix_gpu);
1528: delete a->hmatrix_gpu;
1529: a->hmatrix_gpu = hmatrix_gpu;
1530: }
1531: #endif
1532: }
1533: }
1534: { /* log flops */
1535: double gops, time, perf, dev;
1536: HLibProfile::getHcompressPerf(gops, time, perf, dev);
1537: #if PetscDefined(H2OPUS_USE_GPU)
1538: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1539: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1540: #else
1541: PetscCall(PetscLogFlops(1e9 * gops));
1542: #endif
1543: }
1544: PetscCall(PetscLogEventEnd(MAT_H2Opus_Compress, A, 0, 0, 0));
1545: PetscFunctionReturn(PETSC_SUCCESS);
1546: }
1548: /*@
1549: MatH2OpusSetSamplingMat - Set a matrix to be sampled from matrix-vector products on another matrix to construct a hierarchical matrix.
1551: Input Parameters:
1552: + A - the hierarchical matrix
1553: . B - the matrix to be sampled
1554: . bs - maximum number of samples to be taken concurrently
1555: - tol - relative tolerance for construction
1557: Level: intermediate
1559: Notes:
1560: You need to call `MatAssemblyBegin()` and `MatAssemblyEnd()` to update the hierarchical matrix.
1562: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`
1563: @*/
1564: PetscErrorCode MatH2OpusSetSamplingMat(Mat A, Mat B, PetscInt bs, PetscReal tol)
1565: {
1566: PetscBool ish2opus;
1568: PetscFunctionBegin;
1574: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1575: if (ish2opus) {
1576: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1578: if (!a->sampler) a->sampler = new PetscMatrixSampler();
1579: a->sampler->SetSamplingMat(B);
1580: if (bs > 0) a->bs = bs;
1581: if (tol > 0.) a->rtol = tol;
1582: delete a->kernel;
1583: }
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@
1588: MatCreateH2OpusFromKernel - Creates a `MATH2OPUS` from a user-supplied kernel.
1590: Input Parameters:
1591: + comm - MPI communicator
1592: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1593: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1594: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1595: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1596: . spacedim - dimension of the space coordinates
1597: . coords - coordinates of the points
1598: . cdist - whether or not coordinates are distributed
1599: . kernel - computational kernel (or `NULL`), see `MatH2OpusKernelFn` for the calling sequence
1600: . kernelctx - kernel context
1601: . eta - admissibility condition tolerance
1602: . leafsize - leaf size in cluster tree
1603: - basisord - approximation order for Chebychev interpolation of low-rank blocks
1605: Output Parameter:
1606: . nA - matrix
1608: Options Database Keys:
1609: + -mat_h2opus_leafsize leafsize - Leaf size of cluster tree
1610: . -mat_h2opus_eta eta - Admissibility condition tolerance
1611: . -mat_h2opus_order order - Chebychev approximation order
1612: - -mat_h2opus_normsamples maxsamples - Maximum number of samples to be used when estimating norms
1614: Level: intermediate
1616: Note:
1617: Requires `./configure` with `--download-h2opus` or `--with-h2opus-dir`
1619: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatH2OpusKernelFn`
1620: @*/
1621: PetscErrorCode MatCreateH2OpusFromKernel(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, PetscCtx kernelctx, PetscReal eta, PetscInt leafsize, PetscInt basisord, Mat *nA)
1622: {
1623: Mat A;
1624: Mat_H2OPUS *h2opus;
1625: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1627: PetscFunctionBegin;
1628: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1629: PetscCall(MatCreate(comm, &A));
1630: PetscCall(MatSetSizes(A, m, n, M, N));
1631: PetscCheck(M == N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1632: PetscCall(MatSetType(A, MATH2OPUS));
1633: PetscCall(MatBindToCPU(A, iscpu));
1634: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, kernel, kernelctx));
1636: h2opus = (Mat_H2OPUS *)A->data;
1637: if (eta > 0.) h2opus->eta = eta;
1638: if (leafsize > 0) h2opus->leafsize = leafsize;
1639: if (basisord > 0) h2opus->basisord = basisord;
1641: *nA = A;
1642: PetscFunctionReturn(PETSC_SUCCESS);
1643: }
1645: /*@
1646: MatCreateH2OpusFromMat - Creates a `MATH2OPUS` sampling from a user-supplied operator.
1648: Input Parameters:
1649: + B - the matrix to be sampled
1650: . spacedim - dimension of the space coordinates
1651: . coords - coordinates of the points
1652: . cdist - whether or not coordinates are distributed
1653: . eta - admissibility condition tolerance
1654: . leafsize - leaf size in cluster tree
1655: . maxrank - maximum rank allowed
1656: . bs - maximum number of samples to be taken concurrently
1657: - rtol - relative tolerance for construction
1659: Output Parameter:
1660: . nA - matrix
1662: Options Database Keys:
1663: + -mat_h2opus_leafsize <`PetscInt`> - Leaf size of cluster tree
1664: . -mat_h2opus_eta <`PetscReal`> - Admissibility condition tolerance
1665: . -mat_h2opus_maxrank <`PetscInt`> - Maximum rank when constructed from matvecs
1666: . -mat_h2opus_samples <`PetscInt`> - Maximum number of samples to be taken concurrently when constructing from matvecs
1667: . -mat_h2opus_rtol <`PetscReal`> - Relative tolerance for construction from sampling
1668: . -mat_h2opus_check <`PetscBool`> - Check error when constructing from sampling during MatAssemblyEnd()
1669: . -mat_h2opus_hara_verbose <`PetscBool`> - Verbose output from hara construction
1670: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be when estimating norms
1672: Level: intermediate
1674: Note:
1675: Not available in parallel
1677: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
1678: @*/
1679: PetscErrorCode MatCreateH2OpusFromMat(Mat B, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, PetscReal eta, PetscInt leafsize, PetscInt maxrank, PetscInt bs, PetscReal rtol, Mat *nA)
1680: {
1681: Mat A;
1682: Mat_H2OPUS *h2opus;
1683: MPI_Comm comm;
1684: PetscBool boundtocpu = PETSC_TRUE;
1686: PetscFunctionBegin;
1695: PetscAssertPointer(nA, 10);
1696: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1697: PetscCheck(B->rmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1698: PetscCheck(B->rmap->N == B->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1699: PetscCall(MatCreate(comm, &A));
1700: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1701: #if PetscDefined(H2OPUS_USE_GPU)
1702: {
1703: VecType vtype;
1704: PetscBool isstd, iscuda, iskok;
1706: PetscCall(MatGetVecType(B, &vtype));
1707: PetscCall(PetscStrcmpAny(vtype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
1708: PetscCall(PetscStrcmpAny(vtype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
1709: PetscCall(PetscStrcmpAny(vtype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
1710: PetscCheck(isstd || iscuda || iskok, comm, PETSC_ERR_SUP, "Not for type %s", vtype);
1711: if (iscuda && !B->boundtocpu) boundtocpu = PETSC_FALSE;
1712: if (iskok && PetscDefined(HAVE_MACRO_KOKKOS_ENABLE_CUDA)) boundtocpu = PETSC_FALSE;
1713: }
1714: #endif
1715: PetscCall(MatSetType(A, MATH2OPUS));
1716: PetscCall(MatBindToCPU(A, boundtocpu));
1717: if (spacedim) PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, NULL, NULL));
1718: PetscCall(MatPropagateSymmetryOptions(B, A));
1719: /* PetscCheck(A->symmetric,comm,PETSC_ERR_SUP,"Unsymmetric sampling does not work"); */
1721: h2opus = (Mat_H2OPUS *)A->data;
1722: h2opus->sampler = new PetscMatrixSampler(B);
1723: if (eta > 0.) h2opus->eta = eta;
1724: if (leafsize > 0) h2opus->leafsize = leafsize;
1725: if (maxrank > 0) h2opus->max_rank = maxrank;
1726: if (bs > 0) h2opus->bs = bs;
1727: if (rtol > 0.) h2opus->rtol = rtol;
1728: *nA = A;
1729: A->preallocated = PETSC_TRUE;
1730: PetscFunctionReturn(PETSC_SUCCESS);
1731: }
1733: /*@
1734: MatH2OpusGetIndexMap - Access reordering index set.
1736: Input Parameter:
1737: . A - the matrix
1739: Output Parameter:
1740: . indexmap - the index set for the reordering
1742: Level: intermediate
1744: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1745: @*/
1746: PetscErrorCode MatH2OpusGetIndexMap(Mat A, IS *indexmap)
1747: {
1748: PetscBool ish2opus;
1749: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1751: PetscFunctionBegin;
1754: PetscAssertPointer(indexmap, 2);
1755: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1756: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1757: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1758: *indexmap = a->h2opus_indexmap;
1759: PetscFunctionReturn(PETSC_SUCCESS);
1760: }
1762: /*@
1763: MatH2OpusMapVec - Maps a vector between PETSc and H2Opus ordering
1765: Input Parameters:
1766: + A - the matrix
1767: . nativetopetsc - if true, maps from H2Opus ordering to PETSc ordering. If false, applies the reverse map
1768: - in - the vector to be mapped
1770: Output Parameter:
1771: . out - the newly created mapped vector
1773: Level: intermediate
1775: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1776: @*/
1777: PetscErrorCode MatH2OpusMapVec(Mat A, PetscBool nativetopetsc, Vec in, Vec *out)
1778: {
1779: PetscBool ish2opus;
1780: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1781: PetscScalar *xin, *xout;
1782: PetscBool nm;
1784: PetscFunctionBegin;
1789: PetscAssertPointer(out, 4);
1790: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1791: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1792: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1793: nm = a->nativemult;
1794: PetscCall(MatH2OpusSetNativeMult(A, (PetscBool)!nativetopetsc));
1795: PetscCall(MatCreateVecs(A, out, NULL));
1796: PetscCall(MatH2OpusSetNativeMult(A, nm));
1797: if (!a->sf) { /* same ordering */
1798: PetscCall(VecCopy(in, *out));
1799: PetscFunctionReturn(PETSC_SUCCESS);
1800: }
1801: PetscCall(VecGetArrayRead(in, (const PetscScalar **)&xin));
1802: PetscCall(VecGetArrayWrite(*out, &xout));
1803: if (nativetopetsc) {
1804: PetscCall(PetscSFReduceBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1805: PetscCall(PetscSFReduceEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1806: } else {
1807: PetscCall(PetscSFBcastBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1808: PetscCall(PetscSFBcastEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1809: }
1810: PetscCall(VecRestoreArrayRead(in, (const PetscScalar **)&xin));
1811: PetscCall(VecRestoreArrayWrite(*out, &xout));
1812: PetscFunctionReturn(PETSC_SUCCESS);
1813: }
1815: /*@
1816: MatH2OpusLowRankUpdate - Perform a low-rank update of the form $ A = A + s * U * V^T $
1818: Input Parameters:
1819: + A - the hierarchical `MATH2OPUS` matrix
1820: . s - the scaling factor
1821: . U - the dense low-rank update matrix
1822: - V - (optional) the dense low-rank update matrix (if `NULL`, then `V` = `U` is assumed)
1824: Note:
1825: The `U` and `V` matrices must be in `MATDENSE` dense format
1827: Level: intermediate
1829: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`, `MATDENSE`
1830: @*/
1831: PetscErrorCode MatH2OpusLowRankUpdate(Mat A, Mat U, Mat V, PetscScalar s)
1832: {
1833: PetscBool flg;
1835: PetscFunctionBegin;
1838: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1840: PetscCheckSameComm(A, 1, U, 2);
1841: if (V) {
1843: PetscCheckSameComm(A, 1, V, 3);
1844: }
1847: if (!V) V = U;
1848: 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);
1849: if (!U->cmap->N) PetscFunctionReturn(PETSC_SUCCESS);
1850: PetscCall(PetscLayoutCompare(U->rmap, A->rmap, &flg));
1851: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A and U must have the same row layout");
1852: PetscCall(PetscLayoutCompare(V->rmap, A->cmap, &flg));
1853: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A column layout must match V row column layout");
1854: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &flg));
1855: if (flg) {
1856: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1857: const PetscScalar *u, *v, *uu, *vv;
1858: PetscInt ldu, ldv;
1859: PetscMPIInt size;
1860: #if defined(H2OPUS_USE_MPI)
1861: h2opusHandle_t handle = a->handle->handle;
1862: #else
1863: h2opusHandle_t handle = a->handle;
1864: #endif
1865: PetscBool usesf = (PetscBool)(a->sf && !a->nativemult);
1866: PetscSF usf, vsf;
1868: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1869: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet implemented in parallel");
1870: PetscCall(PetscLogEventBegin(MAT_H2Opus_LR, A, 0, 0, 0));
1871: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)U, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1872: PetscCheck(flg, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "Not for U of type %s", ((PetscObject)U)->type_name);
1873: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)V, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1874: PetscCheck(flg, PetscObjectComm((PetscObject)V), PETSC_ERR_SUP, "Not for V of type %s", ((PetscObject)V)->type_name);
1875: PetscCall(MatDenseGetLDA(U, &ldu));
1876: PetscCall(MatDenseGetLDA(V, &ldv));
1877: PetscCall(MatBoundToCPU(A, &flg));
1878: if (usesf) {
1879: PetscInt n;
1881: PetscCall(MatDenseGetH2OpusStridedSF(U, a->sf, &usf));
1882: PetscCall(MatDenseGetH2OpusStridedSF(V, a->sf, &vsf));
1883: PetscCall(MatH2OpusResizeBuffers_Private(A, U->cmap->N, V->cmap->N));
1884: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
1885: ldu = n;
1886: ldv = n;
1887: }
1888: if (flg) {
1889: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1890: PetscCall(MatDenseGetArrayRead(U, &u));
1891: PetscCall(MatDenseGetArrayRead(V, &v));
1892: if (usesf) {
1893: vv = MatH2OpusGetThrustPointer(*a->yy);
1894: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1895: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1896: if (U != V) {
1897: uu = MatH2OpusGetThrustPointer(*a->xx);
1898: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1899: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1900: } else uu = vv;
1901: } else {
1902: uu = u;
1903: vv = v;
1904: }
1905: hlru_global(*a->hmatrix, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1906: PetscCall(MatDenseRestoreArrayRead(U, &u));
1907: PetscCall(MatDenseRestoreArrayRead(V, &v));
1908: } else {
1909: #if PetscDefined(H2OPUS_USE_GPU)
1910: PetscBool flgU, flgV;
1912: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1913: PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &flgU, MATSEQDENSE, MATMPIDENSE, ""));
1914: if (flgU) PetscCall(MatConvert(U, MATDENSECUDA, MAT_INPLACE_MATRIX, &U));
1915: PetscCall(PetscObjectTypeCompareAny((PetscObject)V, &flgV, MATSEQDENSE, MATMPIDENSE, ""));
1916: if (flgV) PetscCall(MatConvert(V, MATDENSECUDA, MAT_INPLACE_MATRIX, &V));
1917: PetscCall(MatDenseCUDAGetArrayRead(U, &u));
1918: PetscCall(MatDenseCUDAGetArrayRead(V, &v));
1919: if (usesf) {
1920: vv = MatH2OpusGetThrustPointer(*a->yy_gpu);
1921: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1922: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1923: if (U != V) {
1924: uu = MatH2OpusGetThrustPointer(*a->xx_gpu);
1925: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1926: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1927: } else uu = vv;
1928: } else {
1929: uu = u;
1930: vv = v;
1931: }
1932: #else
1933: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
1934: #endif
1935: hlru_global(*a->hmatrix_gpu, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1936: #if PetscDefined(H2OPUS_USE_GPU)
1937: PetscCall(MatDenseCUDARestoreArrayRead(U, &u));
1938: PetscCall(MatDenseCUDARestoreArrayRead(V, &v));
1939: if (flgU) PetscCall(MatConvert(U, MATDENSE, MAT_INPLACE_MATRIX, &U));
1940: if (flgV) PetscCall(MatConvert(V, MATDENSE, MAT_INPLACE_MATRIX, &V));
1941: #endif
1942: }
1943: PetscCall(PetscLogEventEnd(MAT_H2Opus_LR, A, 0, 0, 0));
1944: a->orthogonal = PETSC_FALSE;
1945: }
1946: PetscFunctionReturn(PETSC_SUCCESS);
1947: }
1948: #endif