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