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) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
517: PetscCall(MatSetUp(C));
518: break;
519: case MATPRODUCT_AtB:
520: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
521: PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
522: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
523: if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
524: PetscCall(MatSetUp(C));
525: break;
526: default:
527: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
528: }
529: C->ops->productsymbolic = NULL;
530: C->ops->productnumeric = MatProductNumeric_H2OPUS;
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: static PetscErrorCode MatProductSetFromOptions_H2OPUS(Mat C)
535: {
536: PetscFunctionBegin;
537: MatCheckProduct(C, 1);
538: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_H2OPUS;
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: static PetscErrorCode MatMultKernel_H2OPUS(Mat A, Vec x, PetscScalar sy, Vec y, PetscBool trans)
543: {
544: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
545: #if defined(H2OPUS_USE_MPI)
546: h2opusHandle_t handle = h2opus->handle->handle;
547: #else
548: h2opusHandle_t handle = h2opus->handle;
549: #endif
550: PetscBool boundtocpu = PETSC_TRUE;
551: PetscInt n;
552: PetscScalar *xx, *yy, *uxx, *uyy;
553: PetscMPIInt size;
554: PetscBool usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);
556: PetscFunctionBegin;
557: HLibProfile::clear();
558: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
559: #if PetscDefined(H2OPUS_USE_GPU)
560: boundtocpu = A->boundtocpu;
561: #endif
562: if (usesf) PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
563: else n = A->rmap->n;
564: if (boundtocpu) {
565: PetscCall(VecGetArrayRead(x, (const PetscScalar **)&xx));
566: if (sy == 0.0) {
567: PetscCall(VecGetArrayWrite(y, &yy));
568: } else {
569: PetscCall(VecGetArray(y, &yy));
570: }
571: if (usesf) {
572: uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
573: uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
575: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
576: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
577: if (sy != 0.0) {
578: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
579: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
580: }
581: } else {
582: uxx = xx;
583: uyy = yy;
584: }
585: if (size > 1) {
586: PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
587: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
588: #if defined(H2OPUS_USE_MPI)
589: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, n, sy, uyy, n, 1, h2opus->handle);
590: #endif
591: } else {
592: PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
593: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, n, sy, uyy, n, 1, handle);
594: }
595: PetscCall(VecRestoreArrayRead(x, (const PetscScalar **)&xx));
596: if (usesf) {
597: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
598: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
599: }
600: if (sy == 0.0) {
601: PetscCall(VecRestoreArrayWrite(y, &yy));
602: } else {
603: PetscCall(VecRestoreArray(y, &yy));
604: }
605: #if PetscDefined(H2OPUS_USE_GPU)
606: } else {
607: PetscCall(VecCUDAGetArrayRead(x, (const PetscScalar **)&xx));
608: if (sy == 0.0) {
609: PetscCall(VecCUDAGetArrayWrite(y, &yy));
610: } else {
611: PetscCall(VecCUDAGetArray(y, &yy));
612: }
613: if (usesf) {
614: uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
615: uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
617: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
618: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
619: if (sy != 0.0) {
620: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
621: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
622: }
623: } else {
624: uxx = xx;
625: uyy = yy;
626: }
627: PetscCall(PetscLogGpuTimeBegin());
628: if (size > 1) {
629: PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
630: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
631: #if defined(H2OPUS_USE_MPI)
632: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, n, sy, uyy, n, 1, h2opus->handle);
633: #endif
634: } else {
635: PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
636: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, n, sy, uyy, n, 1, handle);
637: }
638: PetscCall(PetscLogGpuTimeEnd());
639: PetscCall(VecCUDARestoreArrayRead(x, (const PetscScalar **)&xx));
640: if (usesf) {
641: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
642: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
643: }
644: if (sy == 0.0) {
645: PetscCall(VecCUDARestoreArrayWrite(y, &yy));
646: } else {
647: PetscCall(VecCUDARestoreArray(y, &yy));
648: }
649: #endif
650: }
651: { /* log flops */
652: double gops, time, perf, dev;
653: HLibProfile::getHgemvPerf(gops, time, perf, dev);
654: #if PetscDefined(H2OPUS_USE_GPU)
655: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
656: else PetscCall(PetscLogGpuFlops(1e9 * gops));
657: #else
658: PetscCall(PetscLogFlops(1e9 * gops));
659: #endif
660: }
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: static PetscErrorCode MatMultTranspose_H2OPUS(Mat A, Vec x, Vec y)
665: {
666: PetscBool xiscuda, yiscuda;
668: PetscFunctionBegin;
669: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
670: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
671: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
672: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_TRUE));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: static PetscErrorCode MatMult_H2OPUS(Mat A, Vec x, Vec y)
677: {
678: PetscBool xiscuda, yiscuda;
680: PetscFunctionBegin;
681: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
682: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
683: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
684: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_FALSE));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: static PetscErrorCode MatMultTransposeAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
689: {
690: PetscBool xiscuda, ziscuda;
692: PetscFunctionBegin;
693: PetscCall(VecCopy(y, z));
694: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
695: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
696: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
697: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_TRUE));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode MatMultAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
702: {
703: PetscBool xiscuda, ziscuda;
705: PetscFunctionBegin;
706: PetscCall(VecCopy(y, z));
707: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
708: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
709: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
710: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_FALSE));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static PetscErrorCode MatScale_H2OPUS(Mat A, PetscScalar s)
715: {
716: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
718: PetscFunctionBegin;
719: a->s *= s;
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: static PetscErrorCode MatSetFromOptions_H2OPUS(Mat A, PetscOptionItems PetscOptionsObject)
724: {
725: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
727: PetscFunctionBegin;
728: PetscOptionsHeadBegin(PetscOptionsObject, "H2OPUS options");
729: PetscCall(PetscOptionsInt("-mat_h2opus_leafsize", "Leaf size of cluster tree", NULL, a->leafsize, &a->leafsize, NULL));
730: PetscCall(PetscOptionsReal("-mat_h2opus_eta", "Admissibility condition tolerance", NULL, a->eta, &a->eta, NULL));
731: PetscCall(PetscOptionsInt("-mat_h2opus_order", "Basis order for off-diagonal sampling when constructed from kernel", NULL, a->basisord, &a->basisord, NULL));
732: PetscCall(PetscOptionsInt("-mat_h2opus_maxrank", "Maximum rank when constructed from matvecs", NULL, a->max_rank, &a->max_rank, NULL));
733: PetscCall(PetscOptionsInt("-mat_h2opus_samples", "Maximum number of samples to be taken concurrently when constructing from matvecs", NULL, a->bs, &a->bs, NULL));
734: PetscCall(PetscOptionsInt("-mat_h2opus_normsamples", "Maximum number of samples to be when estimating norms", NULL, a->norm_max_samples, &a->norm_max_samples, NULL));
735: PetscCall(PetscOptionsReal("-mat_h2opus_rtol", "Relative tolerance for construction from sampling", NULL, a->rtol, &a->rtol, NULL));
736: PetscCall(PetscOptionsBool("-mat_h2opus_check", "Check error when constructing from sampling during MatAssemblyEnd()", NULL, a->check_construction, &a->check_construction, NULL));
737: PetscCall(PetscOptionsBool("-mat_h2opus_hara_verbose", "Verbose output from hara construction", NULL, a->hara_verbose, &a->hara_verbose, NULL));
738: PetscCall(PetscOptionsBool("-mat_h2opus_resize", "Resize after compression", NULL, a->resize, &a->resize, NULL));
739: PetscOptionsHeadEnd();
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat, PetscInt, const PetscReal[], PetscBool, MatH2OpusKernelFn *, void *);
745: static PetscErrorCode MatH2OpusInferCoordinates_Private(Mat A)
746: {
747: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
748: Vec c;
749: PetscInt spacedim;
750: const PetscScalar *coords;
752: PetscFunctionBegin;
753: if (a->ptcloud) PetscFunctionReturn(PETSC_SUCCESS);
754: PetscCall(PetscObjectQuery((PetscObject)A, "__math2opus_coords", (PetscObject *)&c));
755: if (!c && a->sampler) {
756: Mat S = a->sampler->GetSamplingMat();
758: PetscCall(PetscObjectQuery((PetscObject)S, "__math2opus_coords", (PetscObject *)&c));
759: }
760: if (!c) {
761: PetscCall(MatH2OpusSetCoords_H2OPUS(A, -1, NULL, PETSC_FALSE, NULL, NULL));
762: } else {
763: PetscCall(VecGetArrayRead(c, &coords));
764: PetscCall(VecGetBlockSize(c, &spacedim));
765: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, PETSC_FALSE, NULL, NULL));
766: PetscCall(VecRestoreArrayRead(c, &coords));
767: }
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: static PetscErrorCode MatSetUpMultiply_H2OPUS(Mat A)
772: {
773: MPI_Comm comm;
774: PetscMPIInt size;
775: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
776: PetscInt n = 0, *idx = NULL;
777: int *iidx = NULL;
778: PetscCopyMode own;
779: PetscBool rid;
781: PetscFunctionBegin;
782: if (a->multsetup) PetscFunctionReturn(PETSC_SUCCESS);
783: if (a->sf) { /* MatDuplicate_H2OPUS takes reference to the SF */
784: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
785: #if PetscDefined(H2OPUS_USE_GPU)
786: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
787: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
788: a->xxs_gpu = 1;
789: a->yys_gpu = 1;
790: #endif
791: a->xx = new thrust::host_vector<PetscScalar>(n);
792: a->yy = new thrust::host_vector<PetscScalar>(n);
793: a->xxs = 1;
794: a->yys = 1;
795: } else {
796: IS is;
797: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
798: PetscCallMPI(MPI_Comm_size(comm, &size));
799: if (!a->h2opus_indexmap) {
800: if (size > 1) {
801: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
802: #if defined(H2OPUS_USE_MPI)
803: iidx = MatH2OpusGetThrustPointer(a->dist_hmatrix->basis_tree.basis_branch.index_map);
804: n = a->dist_hmatrix->basis_tree.basis_branch.index_map.size();
805: #endif
806: } else {
807: iidx = MatH2OpusGetThrustPointer(a->hmatrix->u_basis_tree.index_map);
808: n = a->hmatrix->u_basis_tree.index_map.size();
809: }
811: if (PetscDefined(USE_64BIT_INDICES)) {
812: PetscInt i;
814: own = PETSC_OWN_POINTER;
815: PetscCall(PetscMalloc1(n, &idx));
816: for (i = 0; i < n; i++) idx[i] = iidx[i];
817: } else {
818: own = PETSC_COPY_VALUES;
819: idx = (PetscInt *)iidx;
820: }
821: PetscCall(ISCreateGeneral(comm, n, idx, own, &is));
822: PetscCall(ISSetPermutation(is));
823: PetscCall(ISViewFromOptions(is, (PetscObject)A, "-mat_h2opus_indexmap_view"));
824: a->h2opus_indexmap = is;
825: }
826: PetscCall(ISGetLocalSize(a->h2opus_indexmap, &n));
827: PetscCall(ISGetIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
828: rid = (PetscBool)(n == A->rmap->n);
829: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &rid, 1, MPI_C_BOOL, MPI_LAND, comm));
830: if (rid) PetscCall(ISIdentity(a->h2opus_indexmap, &rid));
831: if (!rid) {
832: if (size > 1) { /* Parallel distribution may be different, save it here for fast path in MatMult (see MatH2OpusSetNativeMult) */
833: PetscCall(PetscLayoutCreate(comm, &a->h2opus_rmap));
834: PetscCall(PetscLayoutSetLocalSize(a->h2opus_rmap, n));
835: PetscCall(PetscLayoutSetUp(a->h2opus_rmap));
836: PetscCall(PetscLayoutReference(a->h2opus_rmap, &a->h2opus_cmap));
837: }
838: PetscCall(PetscSFCreate(comm, &a->sf));
839: PetscCall(PetscSFSetGraphLayout(a->sf, A->rmap, n, NULL, PETSC_OWN_POINTER, idx));
840: PetscCall(PetscSFSetUp(a->sf));
841: PetscCall(PetscSFViewFromOptions(a->sf, (PetscObject)A, "-mat_h2opus_sf_view"));
842: #if PetscDefined(H2OPUS_USE_GPU)
843: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
844: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
845: a->xxs_gpu = 1;
846: a->yys_gpu = 1;
847: #endif
848: a->xx = new thrust::host_vector<PetscScalar>(n);
849: a->yy = new thrust::host_vector<PetscScalar>(n);
850: a->xxs = 1;
851: a->yys = 1;
852: }
853: PetscCall(ISRestoreIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
854: }
855: a->multsetup = PETSC_TRUE;
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode MatAssemblyEnd_H2OPUS(Mat A, MatAssemblyType assemblytype)
860: {
861: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
862: #if defined(H2OPUS_USE_MPI)
863: h2opusHandle_t handle = a->handle->handle;
864: #else
865: h2opusHandle_t handle = a->handle;
866: #endif
867: PetscBool kernel = PETSC_FALSE;
868: PetscBool boundtocpu = PETSC_TRUE;
869: PetscBool samplingdone = PETSC_FALSE;
870: MPI_Comm comm;
871: PetscMPIInt size;
873: PetscFunctionBegin;
874: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
875: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
876: PetscCheck(A->rmap->N == A->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
878: /* XXX */
879: a->leafsize = PetscMin(a->leafsize, PetscMin(A->rmap->N, A->cmap->N));
881: PetscCallMPI(MPI_Comm_size(comm, &size));
882: /* TODO REUSABILITY of geometric construction */
883: delete a->hmatrix;
884: delete a->dist_hmatrix;
885: #if PetscDefined(H2OPUS_USE_GPU)
886: delete a->hmatrix_gpu;
887: delete a->dist_hmatrix_gpu;
888: #endif
889: a->orthogonal = PETSC_FALSE;
891: /* TODO: other? */
892: H2OpusBoxCenterAdmissibility adm(a->eta);
894: PetscCall(PetscLogEventBegin(MAT_H2Opus_Build, A, 0, 0, 0));
895: if (size > 1) {
896: #if defined(H2OPUS_USE_MPI)
897: a->dist_hmatrix = new DistributedHMatrix(A->rmap->n /* ,A->symmetric */);
898: #else
899: a->dist_hmatrix = NULL;
900: #endif
901: } else a->hmatrix = new HMatrix(A->rmap->n, A->symmetric == PETSC_BOOL3_TRUE);
902: PetscCall(MatH2OpusInferCoordinates_Private(A));
903: PetscCheck(a->ptcloud, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing pointcloud");
904: if (a->kernel) {
905: BoxEntryGen<PetscScalar, H2OPUS_HWTYPE_CPU, PetscFunctionGenerator<PetscScalar>> entry_gen(*a->kernel);
906: if (size > 1) {
907: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
908: #if defined(H2OPUS_USE_MPI)
909: buildDistributedHMatrix(*a->dist_hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord, a->handle);
910: #endif
911: } else {
912: buildHMatrix(*a->hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord);
913: }
914: kernel = PETSC_TRUE;
915: } else {
916: PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Construction from sampling not supported in parallel");
917: buildHMatrixStructure(*a->hmatrix, a->ptcloud, a->leafsize, adm);
918: }
919: PetscCall(MatSetUpMultiply_H2OPUS(A));
921: #if PetscDefined(H2OPUS_USE_GPU)
922: boundtocpu = A->boundtocpu;
923: if (!boundtocpu) {
924: if (size > 1) {
925: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
926: #if defined(H2OPUS_USE_MPI)
927: a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
928: #endif
929: } else {
930: a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
931: }
932: }
933: #endif
934: if (size == 1) {
935: if (!kernel && a->sampler && a->sampler->GetSamplingMat()) {
936: PetscReal Anorm;
937: bool verbose;
939: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_hara_verbose", &a->hara_verbose, NULL));
940: verbose = a->hara_verbose;
941: PetscCall(MatNormApproximate(a->sampler->GetSamplingMat(), NORM_2, a->norm_max_samples, &Anorm));
942: 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));
943: if (a->sf && !a->nativemult) a->sampler->SetIndexMap(a->hmatrix->u_basis_tree.index_map.size(), a->hmatrix->u_basis_tree.index_map.data());
944: a->sampler->SetStream(handle->getMainStream());
945: if (boundtocpu) {
946: a->sampler->SetGPUSampling(false);
947: hara(a->sampler, *a->hmatrix, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
948: #if PetscDefined(H2OPUS_USE_GPU)
949: } else {
950: a->sampler->SetGPUSampling(true);
951: hara(a->sampler, *a->hmatrix_gpu, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
952: #endif
953: }
954: samplingdone = PETSC_TRUE;
955: }
956: }
957: #if PetscDefined(H2OPUS_USE_GPU)
958: if (!boundtocpu) {
959: delete a->hmatrix;
960: delete a->dist_hmatrix;
961: a->hmatrix = NULL;
962: a->dist_hmatrix = NULL;
963: }
964: A->offloadmask = boundtocpu ? PETSC_OFFLOAD_CPU : PETSC_OFFLOAD_GPU;
965: #endif
966: PetscCall(PetscLogEventEnd(MAT_H2Opus_Build, A, 0, 0, 0));
968: if (!a->s) a->s = 1.0;
969: A->assembled = PETSC_TRUE;
971: if (samplingdone) {
972: PetscBool check = a->check_construction;
973: PetscBool checke = PETSC_FALSE;
975: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check", &check, NULL));
976: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check_explicit", &checke, NULL));
977: if (check) {
978: Mat E, Ae;
979: PetscReal n1, ni, n2;
980: PetscReal n1A, niA, n2A;
981: PetscErrorCodeFn *normfunc;
983: Ae = a->sampler->GetSamplingMat();
984: PetscCall(MatConvert(A, MATSHELL, MAT_INITIAL_MATRIX, &E));
985: PetscCall(MatShellSetOperation(E, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
986: PetscCall(MatAXPY(E, -1.0, Ae, DIFFERENT_NONZERO_PATTERN));
987: PetscCall(MatNorm(E, NORM_1, &n1));
988: PetscCall(MatNorm(E, NORM_INFINITY, &ni));
989: PetscCall(MatNorm(E, NORM_2, &n2));
990: if (checke) {
991: Mat eA, eE, eAe;
993: PetscCall(MatComputeOperator(A, MATAIJ, &eA));
994: PetscCall(MatComputeOperator(E, MATAIJ, &eE));
995: PetscCall(MatComputeOperator(Ae, MATAIJ, &eAe));
996: PetscCall(MatFilter(eA, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
997: PetscCall(MatFilter(eE, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
998: PetscCall(MatFilter(eAe, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
999: PetscCall(PetscObjectSetName((PetscObject)eA, "H2Mat"));
1000: PetscCall(MatView(eA, NULL));
1001: PetscCall(PetscObjectSetName((PetscObject)eAe, "S"));
1002: PetscCall(MatView(eAe, NULL));
1003: PetscCall(PetscObjectSetName((PetscObject)eE, "H2Mat - S"));
1004: PetscCall(MatView(eE, NULL));
1005: PetscCall(MatDestroy(&eA));
1006: PetscCall(MatDestroy(&eE));
1007: PetscCall(MatDestroy(&eAe));
1008: }
1010: PetscCall(MatGetOperation(Ae, MATOP_NORM, &normfunc));
1011: PetscCall(MatSetOperation(Ae, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
1012: PetscCall(MatNorm(Ae, NORM_1, &n1A));
1013: PetscCall(MatNorm(Ae, NORM_INFINITY, &niA));
1014: PetscCall(MatNorm(Ae, NORM_2, &n2A));
1015: n1A = PetscMax(n1A, PETSC_SMALL);
1016: n2A = PetscMax(n2A, PETSC_SMALL);
1017: niA = PetscMax(niA, PETSC_SMALL);
1018: PetscCall(MatSetOperation(Ae, MATOP_NORM, normfunc));
1019: 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)));
1020: PetscCall(MatDestroy(&E));
1021: }
1022: a->sampler->SetSamplingMat(NULL);
1023: }
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: static PetscErrorCode MatZeroEntries_H2OPUS(Mat A)
1028: {
1029: PetscMPIInt size;
1030: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1032: PetscFunctionBegin;
1033: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1034: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet supported");
1035: a->hmatrix->clearData();
1036: #if PetscDefined(H2OPUS_USE_GPU)
1037: if (a->hmatrix_gpu) a->hmatrix_gpu->clearData();
1038: #endif
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: static PetscErrorCode MatDuplicate_H2OPUS(Mat B, MatDuplicateOption op, Mat *nA)
1043: {
1044: Mat A;
1045: Mat_H2OPUS *a, *b = (Mat_H2OPUS *)B->data;
1046: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1047: MPI_Comm comm;
1049: PetscFunctionBegin;
1050: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1051: PetscCall(MatCreate(comm, &A));
1052: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1053: PetscCall(MatSetType(A, MATH2OPUS));
1054: PetscCall(MatPropagateSymmetryOptions(B, A));
1055: a = (Mat_H2OPUS *)A->data;
1057: a->eta = b->eta;
1058: a->leafsize = b->leafsize;
1059: a->basisord = b->basisord;
1060: a->max_rank = b->max_rank;
1061: a->bs = b->bs;
1062: a->rtol = b->rtol;
1063: a->norm_max_samples = b->norm_max_samples;
1064: if (op == MAT_COPY_VALUES) a->s = b->s;
1066: a->ptcloud = new PetscPointCloud<PetscReal>(*b->ptcloud);
1067: if (op == MAT_COPY_VALUES && b->kernel) a->kernel = new PetscFunctionGenerator<PetscScalar>(*b->kernel);
1069: #if defined(H2OPUS_USE_MPI)
1070: if (b->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*b->dist_hmatrix);
1071: #if PetscDefined(H2OPUS_USE_GPU)
1072: if (b->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*b->dist_hmatrix_gpu);
1073: #endif
1074: #endif
1075: if (b->hmatrix) {
1076: a->hmatrix = new HMatrix(*b->hmatrix);
1077: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix->clearData();
1078: }
1079: #if PetscDefined(H2OPUS_USE_GPU)
1080: if (b->hmatrix_gpu) {
1081: a->hmatrix_gpu = new HMatrix_GPU(*b->hmatrix_gpu);
1082: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix_gpu->clearData();
1083: }
1084: #endif
1085: if (b->sf) {
1086: PetscCall(PetscObjectReference((PetscObject)b->sf));
1087: a->sf = b->sf;
1088: }
1089: if (b->h2opus_indexmap) {
1090: PetscCall(PetscObjectReference((PetscObject)b->h2opus_indexmap));
1091: a->h2opus_indexmap = b->h2opus_indexmap;
1092: }
1094: PetscCall(MatSetUp(A));
1095: PetscCall(MatSetUpMultiply_H2OPUS(A));
1096: if (op == MAT_COPY_VALUES) {
1097: A->assembled = PETSC_TRUE;
1098: a->orthogonal = b->orthogonal;
1099: #if PetscDefined(H2OPUS_USE_GPU)
1100: A->offloadmask = B->offloadmask;
1101: #endif
1102: }
1103: #if PetscDefined(H2OPUS_USE_GPU)
1104: iscpu = B->boundtocpu;
1105: #endif
1106: PetscCall(MatBindToCPU(A, iscpu));
1108: *nA = A;
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: static PetscErrorCode MatView_H2OPUS(Mat A, PetscViewer view)
1113: {
1114: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1115: PetscBool isascii, vieweps;
1116: PetscMPIInt size;
1117: PetscViewerFormat format;
1119: PetscFunctionBegin;
1120: PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
1121: PetscCall(PetscViewerGetFormat(view, &format));
1122: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1123: if (isascii) {
1124: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1125: if (size == 1) {
1126: FILE *fp;
1127: PetscCall(PetscViewerASCIIGetPointer(view, &fp));
1128: dumpHMatrix(*h2opus->hmatrix, 6, fp);
1129: }
1130: } else {
1131: PetscCall(PetscViewerASCIIPrintf(view, " H-Matrix constructed from %s\n", h2opus->kernel ? "Kernel" : "Mat"));
1132: PetscCall(PetscViewerASCIIPrintf(view, " PointCloud dim %" PetscInt_FMT "\n", h2opus->ptcloud ? h2opus->ptcloud->getDimension() : 0));
1133: PetscCall(PetscViewerASCIIPrintf(view, " Admissibility parameters: leaf size %" PetscInt_FMT ", eta %g\n", h2opus->leafsize, (double)h2opus->eta));
1134: if (!h2opus->kernel) {
1135: PetscCall(PetscViewerASCIIPrintf(view, " Sampling parameters: max_rank %" PetscInt_FMT ", samples %" PetscInt_FMT ", tolerance %g\n", h2opus->max_rank, h2opus->bs, (double)h2opus->rtol));
1136: } else {
1137: PetscCall(PetscViewerASCIIPrintf(view, " Off-diagonal blocks approximation order %" PetscInt_FMT "\n", h2opus->basisord));
1138: }
1139: PetscCall(PetscViewerASCIIPrintf(view, " Number of samples for norms %" PetscInt_FMT "\n", h2opus->norm_max_samples));
1140: if (size == 1) {
1141: double dense_mem_cpu = h2opus->hmatrix ? h2opus->hmatrix->getDenseMemoryUsage() : 0;
1142: double low_rank_cpu = h2opus->hmatrix ? h2opus->hmatrix->getLowRankMemoryUsage() : 0;
1143: #if PetscDefined(HAVE_CUDA)
1144: double dense_mem_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getDenseMemoryUsage() : 0;
1145: double low_rank_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getLowRankMemoryUsage() : 0;
1146: #endif
1147: 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));
1148: #if PetscDefined(HAVE_CUDA)
1149: 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));
1150: #endif
1151: } else {
1152: #if PetscDefined(HAVE_CUDA)
1153: double matrix_mem[4] = {0., 0., 0., 0.};
1154: PetscMPIInt rsize = 4;
1155: #else
1156: double matrix_mem[2] = {0., 0.};
1157: PetscMPIInt rsize = 2;
1158: #endif
1159: #if defined(H2OPUS_USE_MPI)
1160: matrix_mem[0] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalDenseMemoryUsage() : 0;
1161: matrix_mem[1] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalLowRankMemoryUsage() : 0;
1162: #if PetscDefined(HAVE_CUDA)
1163: matrix_mem[2] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalDenseMemoryUsage() : 0;
1164: matrix_mem[3] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalLowRankMemoryUsage() : 0;
1165: #endif
1166: #endif
1167: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, matrix_mem, rsize, MPI_DOUBLE_PRECISION, MPI_SUM, PetscObjectComm((PetscObject)A)));
1168: 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]));
1169: #if PetscDefined(HAVE_CUDA)
1170: 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]));
1171: #endif
1172: }
1173: }
1174: }
1175: vieweps = PETSC_FALSE;
1176: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps", &vieweps, NULL));
1177: if (vieweps) {
1178: char filename[256];
1179: const char *name;
1181: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1182: PetscCall(PetscSNPrintf(filename, sizeof(filename), "%s_structure.eps", name));
1183: PetscCall(PetscOptionsGetString(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps_filename", filename, sizeof(filename), NULL));
1184: outputEps(*h2opus->hmatrix, filename);
1185: }
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat A, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, void *kernelctx)
1190: {
1191: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1192: PetscReal *gcoords;
1193: PetscInt N;
1194: MPI_Comm comm;
1195: PetscMPIInt size;
1196: PetscBool cong;
1198: PetscFunctionBegin;
1199: PetscCall(PetscLayoutSetUp(A->rmap));
1200: PetscCall(PetscLayoutSetUp(A->cmap));
1201: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1202: PetscCall(MatHasCongruentLayouts(A, &cong));
1203: PetscCheck(cong, comm, PETSC_ERR_SUP, "Only for square matrices with congruent layouts");
1204: N = A->rmap->N;
1205: PetscCallMPI(MPI_Comm_size(comm, &size));
1206: if (spacedim > 0 && size > 1 && cdist) {
1207: PetscSF sf;
1208: MPI_Datatype dtype;
1210: PetscCallMPI(MPI_Type_contiguous(spacedim, MPIU_REAL, &dtype));
1211: PetscCallMPI(MPI_Type_commit(&dtype));
1213: PetscCall(PetscSFCreate(comm, &sf));
1214: PetscCall(PetscSFSetGraphWithPattern(sf, A->rmap, PETSCSF_PATTERN_ALLGATHER));
1215: PetscCall(PetscMalloc1(spacedim * N, &gcoords));
1216: PetscCall(PetscSFBcastBegin(sf, dtype, coords, gcoords, MPI_REPLACE));
1217: PetscCall(PetscSFBcastEnd(sf, dtype, coords, gcoords, MPI_REPLACE));
1218: PetscCall(PetscSFDestroy(&sf));
1219: PetscCallMPI(MPI_Type_free(&dtype));
1220: } else gcoords = (PetscReal *)coords;
1222: delete h2opus->ptcloud;
1223: delete h2opus->kernel;
1224: h2opus->ptcloud = new PetscPointCloud<PetscReal>(spacedim, N, gcoords);
1225: if (kernel) h2opus->kernel = new PetscFunctionGenerator<PetscScalar>(kernel, spacedim, kernelctx);
1226: if (gcoords != coords) PetscCall(PetscFree(gcoords));
1227: A->preallocated = PETSC_TRUE;
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: #if PetscDefined(H2OPUS_USE_GPU)
1232: static PetscErrorCode MatBindToCPU_H2OPUS(Mat A, PetscBool flg)
1233: {
1234: PetscMPIInt size;
1235: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1237: PetscFunctionBegin;
1238: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1239: if (flg && A->offloadmask == PETSC_OFFLOAD_GPU) {
1240: if (size > 1) {
1241: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1242: #if defined(H2OPUS_USE_MPI)
1243: if (!a->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix_gpu);
1244: else *a->dist_hmatrix = *a->dist_hmatrix_gpu;
1245: #endif
1246: } else {
1247: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1248: if (!a->hmatrix) a->hmatrix = new HMatrix(*a->hmatrix_gpu);
1249: else *a->hmatrix = *a->hmatrix_gpu;
1250: }
1251: delete a->hmatrix_gpu;
1252: delete a->dist_hmatrix_gpu;
1253: a->hmatrix_gpu = NULL;
1254: a->dist_hmatrix_gpu = NULL;
1255: A->offloadmask = PETSC_OFFLOAD_CPU;
1256: } else if (!flg && A->offloadmask == PETSC_OFFLOAD_CPU) {
1257: if (size > 1) {
1258: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1259: #if defined(H2OPUS_USE_MPI)
1260: if (!a->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
1261: else *a->dist_hmatrix_gpu = *a->dist_hmatrix;
1262: #endif
1263: } else {
1264: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1265: if (!a->hmatrix_gpu) a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
1266: else *a->hmatrix_gpu = *a->hmatrix;
1267: }
1268: delete a->hmatrix;
1269: delete a->dist_hmatrix;
1270: a->hmatrix = NULL;
1271: a->dist_hmatrix = NULL;
1272: A->offloadmask = PETSC_OFFLOAD_GPU;
1273: }
1274: PetscCall(PetscFree(A->defaultvectype));
1275: if (!flg) {
1276: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1277: } else {
1278: PetscCall(PetscStrallocpy(VECSTANDARD, &A->defaultvectype));
1279: }
1280: A->boundtocpu = flg;
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1283: #endif
1285: /*MC
1286: MATH2OPUS = "h2opus" - A matrix type for hierarchical matrices using the H2Opus package {cite}`zampinibouakaramturkiyyahkniokeyes2022`.
1288: Options Database Key:
1289: . -mat_type h2opus - matrix type to "h2opus"
1291: Level: beginner
1293: Notes:
1294: H2Opus implements hierarchical matrices in the $H^2$ flavor. It supports CPU or NVIDIA GPUs.
1296: For CPU only builds, use `./configure --download-h2opus --download-thrust` to install PETSc to use H2Opus.
1297: In order to run on NVIDIA GPUs, use `./configure --download-h2opus --download-magma --download-kblas`.
1299: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATHTOOL`, `MATDENSE`, `MatCreateH2OpusFromKernel()`, `MatCreateH2OpusFromMat()`
1300: M*/
1301: PETSC_EXTERN PetscErrorCode MatCreate_H2OPUS(Mat A)
1302: {
1303: Mat_H2OPUS *a;
1304: PetscMPIInt size;
1306: PetscFunctionBegin;
1307: #if PetscDefined(H2OPUS_USE_GPU)
1308: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
1309: #endif
1310: PetscCall(PetscNew(&a));
1311: A->data = (void *)a;
1313: a->eta = 0.9;
1314: a->leafsize = 32;
1315: a->basisord = 4;
1316: a->max_rank = 64;
1317: a->bs = 32;
1318: a->rtol = 1.e-4;
1319: a->s = 1.0;
1320: a->norm_max_samples = 10;
1321: a->resize = PETSC_TRUE; /* reallocate after compression */
1322: #if defined(H2OPUS_USE_MPI)
1323: h2opusCreateDistributedHandleComm(&a->handle, PetscObjectComm((PetscObject)A));
1324: #else
1325: h2opusCreateHandle(&a->handle);
1326: #endif
1327: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1328: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATH2OPUS));
1329: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
1331: A->ops->destroy = MatDestroy_H2OPUS;
1332: A->ops->view = MatView_H2OPUS;
1333: A->ops->assemblyend = MatAssemblyEnd_H2OPUS;
1334: A->ops->mult = MatMult_H2OPUS;
1335: A->ops->multtranspose = MatMultTranspose_H2OPUS;
1336: A->ops->multadd = MatMultAdd_H2OPUS;
1337: A->ops->multtransposeadd = MatMultTransposeAdd_H2OPUS;
1338: A->ops->scale = MatScale_H2OPUS;
1339: A->ops->duplicate = MatDuplicate_H2OPUS;
1340: A->ops->setfromoptions = MatSetFromOptions_H2OPUS;
1341: A->ops->norm = MatNorm_H2OPUS;
1342: A->ops->zeroentries = MatZeroEntries_H2OPUS;
1343: #if PetscDefined(H2OPUS_USE_GPU)
1344: A->ops->bindtocpu = MatBindToCPU_H2OPUS;
1345: #endif
1347: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", MatProductSetFromOptions_H2OPUS));
1348: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", MatProductSetFromOptions_H2OPUS));
1349: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", MatProductSetFromOptions_H2OPUS));
1350: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", MatProductSetFromOptions_H2OPUS));
1351: #if PetscDefined(H2OPUS_USE_GPU)
1352: PetscCall(PetscFree(A->defaultvectype));
1353: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1354: #endif
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: /*@
1359: MatH2OpusOrthogonalize - Orthogonalize the basis tree of a hierarchical matrix.
1361: Input Parameter:
1362: . A - the matrix
1364: Level: intermediate
1366: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`
1367: @*/
1368: PetscErrorCode MatH2OpusOrthogonalize(Mat A)
1369: {
1370: PetscBool ish2opus;
1371: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1372: PetscMPIInt size;
1373: PetscBool boundtocpu = PETSC_TRUE;
1375: PetscFunctionBegin;
1378: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1379: if (!ish2opus) PetscFunctionReturn(PETSC_SUCCESS);
1380: if (a->orthogonal) PetscFunctionReturn(PETSC_SUCCESS);
1381: HLibProfile::clear();
1382: PetscCall(PetscLogEventBegin(MAT_H2Opus_Orthog, A, 0, 0, 0));
1383: #if PetscDefined(H2OPUS_USE_GPU)
1384: boundtocpu = A->boundtocpu;
1385: #endif
1386: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1387: if (size > 1) {
1388: if (boundtocpu) {
1389: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1390: #if defined(H2OPUS_USE_MPI)
1391: distributed_horthog(*a->dist_hmatrix, a->handle);
1392: #endif
1393: #if PetscDefined(H2OPUS_USE_GPU)
1394: A->offloadmask = PETSC_OFFLOAD_CPU;
1395: } else {
1396: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1397: PetscCall(PetscLogGpuTimeBegin());
1398: #if defined(H2OPUS_USE_MPI)
1399: distributed_horthog(*a->dist_hmatrix_gpu, a->handle);
1400: #endif
1401: PetscCall(PetscLogGpuTimeEnd());
1402: #endif
1403: }
1404: } else {
1405: #if defined(H2OPUS_USE_MPI)
1406: h2opusHandle_t handle = a->handle->handle;
1407: #else
1408: h2opusHandle_t handle = a->handle;
1409: #endif
1410: if (boundtocpu) {
1411: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1412: horthog(*a->hmatrix, handle);
1413: #if PetscDefined(H2OPUS_USE_GPU)
1414: A->offloadmask = PETSC_OFFLOAD_CPU;
1415: } else {
1416: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1417: PetscCall(PetscLogGpuTimeBegin());
1418: horthog(*a->hmatrix_gpu, handle);
1419: PetscCall(PetscLogGpuTimeEnd());
1420: #endif
1421: }
1422: }
1423: a->orthogonal = PETSC_TRUE;
1424: { /* log flops */
1425: double gops, time, perf, dev;
1426: HLibProfile::getHorthogPerf(gops, time, perf, dev);
1427: #if PetscDefined(H2OPUS_USE_GPU)
1428: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1429: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1430: #else
1431: PetscCall(PetscLogFlops(1e9 * gops));
1432: #endif
1433: }
1434: PetscCall(PetscLogEventEnd(MAT_H2Opus_Orthog, A, 0, 0, 0));
1435: PetscFunctionReturn(PETSC_SUCCESS);
1436: }
1438: /*@
1439: MatH2OpusCompress - Compress a hierarchical matrix.
1441: Input Parameters:
1442: + A - the matrix
1443: - tol - the absolute truncation threshold
1445: Level: intermediate
1447: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusOrthogonalize()`
1448: @*/
1449: PetscErrorCode MatH2OpusCompress(Mat A, PetscReal tol)
1450: {
1451: PetscBool ish2opus;
1452: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1453: PetscMPIInt size;
1454: PetscBool boundtocpu = PETSC_TRUE;
1456: PetscFunctionBegin;
1460: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1461: if (!ish2opus || tol <= 0.0) PetscFunctionReturn(PETSC_SUCCESS);
1462: PetscCall(MatH2OpusOrthogonalize(A));
1463: HLibProfile::clear();
1464: PetscCall(PetscLogEventBegin(MAT_H2Opus_Compress, A, 0, 0, 0));
1465: #if PetscDefined(H2OPUS_USE_GPU)
1466: boundtocpu = A->boundtocpu;
1467: #endif
1468: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1469: if (size > 1) {
1470: if (boundtocpu) {
1471: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1472: #if defined(H2OPUS_USE_MPI)
1473: distributed_hcompress(*a->dist_hmatrix, tol, a->handle);
1474: if (a->resize) {
1475: DistributedHMatrix *dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix);
1476: delete a->dist_hmatrix;
1477: a->dist_hmatrix = dist_hmatrix;
1478: }
1479: #endif
1480: #if PetscDefined(H2OPUS_USE_GPU)
1481: A->offloadmask = PETSC_OFFLOAD_CPU;
1482: } else {
1483: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1484: PetscCall(PetscLogGpuTimeBegin());
1485: #if defined(H2OPUS_USE_MPI)
1486: distributed_hcompress(*a->dist_hmatrix_gpu, tol, a->handle);
1488: if (a->resize) {
1489: DistributedHMatrix_GPU *dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix_gpu);
1490: delete a->dist_hmatrix_gpu;
1491: a->dist_hmatrix_gpu = dist_hmatrix_gpu;
1492: }
1493: #endif
1494: PetscCall(PetscLogGpuTimeEnd());
1495: #endif
1496: }
1497: } else {
1498: #if defined(H2OPUS_USE_MPI)
1499: h2opusHandle_t handle = a->handle->handle;
1500: #else
1501: h2opusHandle_t handle = a->handle;
1502: #endif
1503: if (boundtocpu) {
1504: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1505: hcompress(*a->hmatrix, tol, handle);
1507: if (a->resize) {
1508: HMatrix *hmatrix = new HMatrix(*a->hmatrix);
1509: delete a->hmatrix;
1510: a->hmatrix = hmatrix;
1511: }
1512: #if PetscDefined(H2OPUS_USE_GPU)
1513: A->offloadmask = PETSC_OFFLOAD_CPU;
1514: } else {
1515: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1516: PetscCall(PetscLogGpuTimeBegin());
1517: hcompress(*a->hmatrix_gpu, tol, handle);
1518: PetscCall(PetscLogGpuTimeEnd());
1520: if (a->resize) {
1521: HMatrix_GPU *hmatrix_gpu = new HMatrix_GPU(*a->hmatrix_gpu);
1522: delete a->hmatrix_gpu;
1523: a->hmatrix_gpu = hmatrix_gpu;
1524: }
1525: #endif
1526: }
1527: }
1528: { /* log flops */
1529: double gops, time, perf, dev;
1530: HLibProfile::getHcompressPerf(gops, time, perf, dev);
1531: #if PetscDefined(H2OPUS_USE_GPU)
1532: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1533: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1534: #else
1535: PetscCall(PetscLogFlops(1e9 * gops));
1536: #endif
1537: }
1538: PetscCall(PetscLogEventEnd(MAT_H2Opus_Compress, A, 0, 0, 0));
1539: PetscFunctionReturn(PETSC_SUCCESS);
1540: }
1542: /*@
1543: MatH2OpusSetSamplingMat - Set a matrix to be sampled from matrix-vector products on another matrix to construct a hierarchical matrix.
1545: Input Parameters:
1546: + A - the hierarchical matrix
1547: . B - the matrix to be sampled
1548: . bs - maximum number of samples to be taken concurrently
1549: - tol - relative tolerance for construction
1551: Level: intermediate
1553: Notes:
1554: You need to call `MatAssemblyBegin()` and `MatAssemblyEnd()` to update the hierarchical matrix.
1556: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`
1557: @*/
1558: PetscErrorCode MatH2OpusSetSamplingMat(Mat A, Mat B, PetscInt bs, PetscReal tol)
1559: {
1560: PetscBool ish2opus;
1562: PetscFunctionBegin;
1568: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1569: if (ish2opus) {
1570: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1572: if (!a->sampler) a->sampler = new PetscMatrixSampler();
1573: a->sampler->SetSamplingMat(B);
1574: if (bs > 0) a->bs = bs;
1575: if (tol > 0.) a->rtol = tol;
1576: delete a->kernel;
1577: }
1578: PetscFunctionReturn(PETSC_SUCCESS);
1579: }
1581: /*@C
1582: MatCreateH2OpusFromKernel - Creates a `MATH2OPUS` from a user-supplied kernel.
1584: Input Parameters:
1585: + comm - MPI communicator
1586: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1587: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1588: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1589: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1590: . spacedim - dimension of the space coordinates
1591: . coords - coordinates of the points
1592: . cdist - whether or not coordinates are distributed
1593: . kernel - computational kernel (or `NULL`)
1594: . kernelctx - kernel context
1595: . eta - admissibility condition tolerance
1596: . leafsize - leaf size in cluster tree
1597: - basisord - approximation order for Chebychev interpolation of low-rank blocks
1599: Output Parameter:
1600: . nA - matrix
1602: Options Database Keys:
1603: + -mat_h2opus_leafsize <`PetscInt`> - Leaf size of cluster tree
1604: . -mat_h2opus_eta <`PetscReal`> - Admissibility condition tolerance
1605: . -mat_h2opus_order <`PetscInt`> - Chebychev approximation order
1606: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be used when estimating norms
1608: Level: intermediate
1610: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`
1611: @*/
1612: 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)
1613: {
1614: Mat A;
1615: Mat_H2OPUS *h2opus;
1616: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1618: PetscFunctionBegin;
1619: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1620: PetscCall(MatCreate(comm, &A));
1621: PetscCall(MatSetSizes(A, m, n, M, N));
1622: PetscCheck(M == N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1623: PetscCall(MatSetType(A, MATH2OPUS));
1624: PetscCall(MatBindToCPU(A, iscpu));
1625: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, kernel, kernelctx));
1627: h2opus = (Mat_H2OPUS *)A->data;
1628: if (eta > 0.) h2opus->eta = eta;
1629: if (leafsize > 0) h2opus->leafsize = leafsize;
1630: if (basisord > 0) h2opus->basisord = basisord;
1632: *nA = A;
1633: PetscFunctionReturn(PETSC_SUCCESS);
1634: }
1636: /*@
1637: MatCreateH2OpusFromMat - Creates a `MATH2OPUS` sampling from a user-supplied operator.
1639: Input Parameters:
1640: + B - the matrix to be sampled
1641: . spacedim - dimension of the space coordinates
1642: . coords - coordinates of the points
1643: . cdist - whether or not coordinates are distributed
1644: . eta - admissibility condition tolerance
1645: . leafsize - leaf size in cluster tree
1646: . maxrank - maximum rank allowed
1647: . bs - maximum number of samples to be taken concurrently
1648: - rtol - relative tolerance for construction
1650: Output Parameter:
1651: . nA - matrix
1653: Options Database Keys:
1654: + -mat_h2opus_leafsize <`PetscInt`> - Leaf size of cluster tree
1655: . -mat_h2opus_eta <`PetscReal`> - Admissibility condition tolerance
1656: . -mat_h2opus_maxrank <`PetscInt`> - Maximum rank when constructed from matvecs
1657: . -mat_h2opus_samples <`PetscInt`> - Maximum number of samples to be taken concurrently when constructing from matvecs
1658: . -mat_h2opus_rtol <`PetscReal`> - Relative tolerance for construction from sampling
1659: . -mat_h2opus_check <`PetscBool`> - Check error when constructing from sampling during MatAssemblyEnd()
1660: . -mat_h2opus_hara_verbose <`PetscBool`> - Verbose output from hara construction
1661: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be when estimating norms
1663: Level: intermediate
1665: Note:
1666: Not available in parallel
1668: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
1669: @*/
1670: PetscErrorCode MatCreateH2OpusFromMat(Mat B, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, PetscReal eta, PetscInt leafsize, PetscInt maxrank, PetscInt bs, PetscReal rtol, Mat *nA)
1671: {
1672: Mat A;
1673: Mat_H2OPUS *h2opus;
1674: MPI_Comm comm;
1675: PetscBool boundtocpu = PETSC_TRUE;
1677: PetscFunctionBegin;
1686: PetscAssertPointer(nA, 10);
1687: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1688: PetscCheck(B->rmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1689: PetscCheck(B->rmap->N == B->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1690: PetscCall(MatCreate(comm, &A));
1691: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1692: #if PetscDefined(H2OPUS_USE_GPU)
1693: {
1694: VecType vtype;
1695: PetscBool isstd, iscuda, iskok;
1697: PetscCall(MatGetVecType(B, &vtype));
1698: PetscCall(PetscStrcmpAny(vtype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
1699: PetscCall(PetscStrcmpAny(vtype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
1700: PetscCall(PetscStrcmpAny(vtype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
1701: PetscCheck(isstd || iscuda || iskok, comm, PETSC_ERR_SUP, "Not for type %s", vtype);
1702: if (iscuda && !B->boundtocpu) boundtocpu = PETSC_FALSE;
1703: if (iskok && PetscDefined(HAVE_MACRO_KOKKOS_ENABLE_CUDA)) boundtocpu = PETSC_FALSE;
1704: }
1705: #endif
1706: PetscCall(MatSetType(A, MATH2OPUS));
1707: PetscCall(MatBindToCPU(A, boundtocpu));
1708: if (spacedim) PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, NULL, NULL));
1709: PetscCall(MatPropagateSymmetryOptions(B, A));
1710: /* PetscCheck(A->symmetric,comm,PETSC_ERR_SUP,"Unsymmetric sampling does not work"); */
1712: h2opus = (Mat_H2OPUS *)A->data;
1713: h2opus->sampler = new PetscMatrixSampler(B);
1714: if (eta > 0.) h2opus->eta = eta;
1715: if (leafsize > 0) h2opus->leafsize = leafsize;
1716: if (maxrank > 0) h2opus->max_rank = maxrank;
1717: if (bs > 0) h2opus->bs = bs;
1718: if (rtol > 0.) h2opus->rtol = rtol;
1719: *nA = A;
1720: A->preallocated = PETSC_TRUE;
1721: PetscFunctionReturn(PETSC_SUCCESS);
1722: }
1724: /*@
1725: MatH2OpusGetIndexMap - Access reordering index set.
1727: Input Parameter:
1728: . A - the matrix
1730: Output Parameter:
1731: . indexmap - the index set for the reordering
1733: Level: intermediate
1735: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1736: @*/
1737: PetscErrorCode MatH2OpusGetIndexMap(Mat A, IS *indexmap)
1738: {
1739: PetscBool ish2opus;
1740: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1742: PetscFunctionBegin;
1745: PetscAssertPointer(indexmap, 2);
1746: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1747: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1748: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1749: *indexmap = a->h2opus_indexmap;
1750: PetscFunctionReturn(PETSC_SUCCESS);
1751: }
1753: /*@
1754: MatH2OpusMapVec - Maps a vector between PETSc and H2Opus ordering
1756: Input Parameters:
1757: + A - the matrix
1758: . nativetopetsc - if true, maps from H2Opus ordering to PETSc ordering. If false, applies the reverse map
1759: - in - the vector to be mapped
1761: Output Parameter:
1762: . out - the newly created mapped vector
1764: Level: intermediate
1766: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1767: @*/
1768: PetscErrorCode MatH2OpusMapVec(Mat A, PetscBool nativetopetsc, Vec in, Vec *out)
1769: {
1770: PetscBool ish2opus;
1771: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1772: PetscScalar *xin, *xout;
1773: PetscBool nm;
1775: PetscFunctionBegin;
1780: PetscAssertPointer(out, 4);
1781: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1782: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1783: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1784: nm = a->nativemult;
1785: PetscCall(MatH2OpusSetNativeMult(A, (PetscBool)!nativetopetsc));
1786: PetscCall(MatCreateVecs(A, out, NULL));
1787: PetscCall(MatH2OpusSetNativeMult(A, nm));
1788: if (!a->sf) { /* same ordering */
1789: PetscCall(VecCopy(in, *out));
1790: PetscFunctionReturn(PETSC_SUCCESS);
1791: }
1792: PetscCall(VecGetArrayRead(in, (const PetscScalar **)&xin));
1793: PetscCall(VecGetArrayWrite(*out, &xout));
1794: if (nativetopetsc) {
1795: PetscCall(PetscSFReduceBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1796: PetscCall(PetscSFReduceEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1797: } else {
1798: PetscCall(PetscSFBcastBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1799: PetscCall(PetscSFBcastEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1800: }
1801: PetscCall(VecRestoreArrayRead(in, (const PetscScalar **)&xin));
1802: PetscCall(VecRestoreArrayWrite(*out, &xout));
1803: PetscFunctionReturn(PETSC_SUCCESS);
1804: }
1806: /*@
1807: MatH2OpusLowRankUpdate - Perform a low-rank update of the form $ A = A + s * U * V^T $
1809: Input Parameters:
1810: + A - the hierarchical `MATH2OPUS` matrix
1811: . s - the scaling factor
1812: . U - the dense low-rank update matrix
1813: - V - (optional) the dense low-rank update matrix (if `NULL`, then `V` = `U` is assumed)
1815: Note:
1816: The `U` and `V` matrices must be in `MATDENSE` dense format
1818: Level: intermediate
1820: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`, `MATDENSE`
1821: @*/
1822: PetscErrorCode MatH2OpusLowRankUpdate(Mat A, Mat U, Mat V, PetscScalar s)
1823: {
1824: PetscBool flg;
1826: PetscFunctionBegin;
1829: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1831: PetscCheckSameComm(A, 1, U, 2);
1832: if (V) {
1834: PetscCheckSameComm(A, 1, V, 3);
1835: }
1838: if (!V) V = U;
1839: 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);
1840: if (!U->cmap->N) PetscFunctionReturn(PETSC_SUCCESS);
1841: PetscCall(PetscLayoutCompare(U->rmap, A->rmap, &flg));
1842: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A and U must have the same row layout");
1843: PetscCall(PetscLayoutCompare(V->rmap, A->cmap, &flg));
1844: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A column layout must match V row column layout");
1845: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &flg));
1846: if (flg) {
1847: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1848: const PetscScalar *u, *v, *uu, *vv;
1849: PetscInt ldu, ldv;
1850: PetscMPIInt size;
1851: #if defined(H2OPUS_USE_MPI)
1852: h2opusHandle_t handle = a->handle->handle;
1853: #else
1854: h2opusHandle_t handle = a->handle;
1855: #endif
1856: PetscBool usesf = (PetscBool)(a->sf && !a->nativemult);
1857: PetscSF usf, vsf;
1859: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1860: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet implemented in parallel");
1861: PetscCall(PetscLogEventBegin(MAT_H2Opus_LR, A, 0, 0, 0));
1862: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)U, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1863: PetscCheck(flg, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "Not for U of type %s", ((PetscObject)U)->type_name);
1864: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)V, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1865: PetscCheck(flg, PetscObjectComm((PetscObject)V), PETSC_ERR_SUP, "Not for V of type %s", ((PetscObject)V)->type_name);
1866: PetscCall(MatDenseGetLDA(U, &ldu));
1867: PetscCall(MatDenseGetLDA(V, &ldv));
1868: PetscCall(MatBoundToCPU(A, &flg));
1869: if (usesf) {
1870: PetscInt n;
1872: PetscCall(MatDenseGetH2OpusStridedSF(U, a->sf, &usf));
1873: PetscCall(MatDenseGetH2OpusStridedSF(V, a->sf, &vsf));
1874: PetscCall(MatH2OpusResizeBuffers_Private(A, U->cmap->N, V->cmap->N));
1875: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
1876: ldu = n;
1877: ldv = n;
1878: }
1879: if (flg) {
1880: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1881: PetscCall(MatDenseGetArrayRead(U, &u));
1882: PetscCall(MatDenseGetArrayRead(V, &v));
1883: if (usesf) {
1884: vv = MatH2OpusGetThrustPointer(*a->yy);
1885: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1886: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1887: if (U != V) {
1888: uu = MatH2OpusGetThrustPointer(*a->xx);
1889: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1890: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1891: } else uu = vv;
1892: } else {
1893: uu = u;
1894: vv = v;
1895: }
1896: hlru_global(*a->hmatrix, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1897: PetscCall(MatDenseRestoreArrayRead(U, &u));
1898: PetscCall(MatDenseRestoreArrayRead(V, &v));
1899: } else {
1900: #if PetscDefined(H2OPUS_USE_GPU)
1901: PetscBool flgU, flgV;
1903: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1904: PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &flgU, MATSEQDENSE, MATMPIDENSE, ""));
1905: if (flgU) PetscCall(MatConvert(U, MATDENSECUDA, MAT_INPLACE_MATRIX, &U));
1906: PetscCall(PetscObjectTypeCompareAny((PetscObject)V, &flgV, MATSEQDENSE, MATMPIDENSE, ""));
1907: if (flgV) PetscCall(MatConvert(V, MATDENSECUDA, MAT_INPLACE_MATRIX, &V));
1908: PetscCall(MatDenseCUDAGetArrayRead(U, &u));
1909: PetscCall(MatDenseCUDAGetArrayRead(V, &v));
1910: if (usesf) {
1911: vv = MatH2OpusGetThrustPointer(*a->yy_gpu);
1912: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1913: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1914: if (U != V) {
1915: uu = MatH2OpusGetThrustPointer(*a->xx_gpu);
1916: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1917: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1918: } else uu = vv;
1919: } else {
1920: uu = u;
1921: vv = v;
1922: }
1923: #else
1924: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
1925: #endif
1926: hlru_global(*a->hmatrix_gpu, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1927: #if PetscDefined(H2OPUS_USE_GPU)
1928: PetscCall(MatDenseCUDARestoreArrayRead(U, &u));
1929: PetscCall(MatDenseCUDARestoreArrayRead(V, &v));
1930: if (flgU) PetscCall(MatConvert(U, MATDENSE, MAT_INPLACE_MATRIX, &U));
1931: if (flgV) PetscCall(MatConvert(V, MATDENSE, MAT_INPLACE_MATRIX, &V));
1932: #endif
1933: }
1934: PetscCall(PetscLogEventEnd(MAT_H2Opus_LR, A, 0, 0, 0));
1935: a->orthogonal = PETSC_FALSE;
1936: }
1937: PetscFunctionReturn(PETSC_SUCCESS);
1938: }
1939: #endif