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 defined(PETSC_HAVE_CUDA) && defined(H2OPUS_USE_GPU)
25: #define PETSC_H2OPUS_USE_GPU
26: #endif
27: #if defined(PETSC_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 defined(PETSC_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 defined(PETSC_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: PetscErrorCode MatH2OpusSetNativeMult(Mat A, PetscBool nm)
232: {
233: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
234: PetscBool ish2opus;
236: PetscFunctionBegin;
239: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
240: if (ish2opus) {
241: if (a->h2opus_rmap) { /* need to swap layouts for vector creation */
242: if ((!a->nativemult && nm) || (a->nativemult && !nm)) {
243: PetscLayout t;
244: t = A->rmap;
245: A->rmap = a->h2opus_rmap;
246: a->h2opus_rmap = t;
247: t = A->cmap;
248: A->cmap = a->h2opus_cmap;
249: a->h2opus_cmap = t;
250: }
251: }
252: a->nativemult = nm;
253: }
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: PetscErrorCode MatH2OpusGetNativeMult(Mat A, PetscBool *nm)
258: {
259: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
260: PetscBool ish2opus;
262: PetscFunctionBegin;
264: PetscAssertPointer(nm, 2);
265: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
266: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
267: *nm = a->nativemult;
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: PETSC_EXTERN PetscErrorCode MatNorm_H2OPUS(Mat A, NormType normtype, PetscReal *n)
272: {
273: PetscBool ish2opus;
274: PetscInt nmax = PETSC_DECIDE;
275: Mat_H2OPUS *a = NULL;
276: PetscBool mult = PETSC_FALSE;
278: PetscFunctionBegin;
279: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
280: if (ish2opus) { /* set userdefine number of samples and fastpath for mult (norms are order independent) */
281: a = (Mat_H2OPUS *)A->data;
283: nmax = a->norm_max_samples;
284: mult = a->nativemult;
285: PetscCall(MatH2OpusSetNativeMult(A, PETSC_TRUE));
286: } else {
287: PetscCall(PetscOptionsGetInt(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_approximate_norm_samples", &nmax, NULL));
288: }
289: PetscCall(MatNormApproximate(A, normtype, nmax, n));
290: if (a) PetscCall(MatH2OpusSetNativeMult(A, mult));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: static PetscErrorCode MatH2OpusResizeBuffers_Private(Mat A, PetscInt xN, PetscInt yN)
295: {
296: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
297: PetscInt n;
298: PetscBool boundtocpu = PETSC_TRUE;
300: PetscFunctionBegin;
301: #if defined(PETSC_H2OPUS_USE_GPU)
302: boundtocpu = A->boundtocpu;
303: #endif
304: PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
305: if (boundtocpu) {
306: if (h2opus->xxs < xN) {
307: h2opus->xx->resize(n * xN);
308: h2opus->xxs = xN;
309: }
310: if (h2opus->yys < yN) {
311: h2opus->yy->resize(n * yN);
312: h2opus->yys = yN;
313: }
314: }
315: #if defined(PETSC_H2OPUS_USE_GPU)
316: if (!boundtocpu) {
317: if (h2opus->xxs_gpu < xN) {
318: h2opus->xx_gpu->resize(n * xN);
319: h2opus->xxs_gpu = xN;
320: }
321: if (h2opus->yys_gpu < yN) {
322: h2opus->yy_gpu->resize(n * yN);
323: h2opus->yys_gpu = yN;
324: }
325: }
326: #endif
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode MatMultNKernel_H2OPUS(Mat A, PetscBool transA, Mat B, Mat C)
331: {
332: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
333: #if defined(H2OPUS_USE_MPI)
334: h2opusHandle_t handle = h2opus->handle->handle;
335: #else
336: h2opusHandle_t handle = h2opus->handle;
337: #endif
338: PetscBool boundtocpu = PETSC_TRUE;
339: PetscScalar *xx, *yy, *uxx, *uyy;
340: PetscInt blda, clda;
341: PetscMPIInt size;
342: PetscSF bsf, csf;
343: PetscBool usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);
345: PetscFunctionBegin;
346: HLibProfile::clear();
347: #if defined(PETSC_H2OPUS_USE_GPU)
348: boundtocpu = A->boundtocpu;
349: #endif
350: PetscCall(MatDenseGetLDA(B, &blda));
351: PetscCall(MatDenseGetLDA(C, &clda));
352: if (usesf) {
353: PetscInt n;
355: PetscCall(MatDenseGetH2OpusStridedSF(B, h2opus->sf, &bsf));
356: PetscCall(MatDenseGetH2OpusStridedSF(C, h2opus->sf, &csf));
358: PetscCall(MatH2OpusResizeBuffers_Private(A, B->cmap->N, C->cmap->N));
359: PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
360: blda = n;
361: clda = n;
362: }
363: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
364: if (boundtocpu) {
365: PetscCall(MatDenseGetArrayRead(B, (const PetscScalar **)&xx));
366: PetscCall(MatDenseGetArrayWrite(C, &yy));
367: if (usesf) {
368: uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
369: uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
370: PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
371: PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
372: } else {
373: uxx = xx;
374: uyy = yy;
375: }
376: if (size > 1) {
377: PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
378: PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
379: #if defined(H2OPUS_USE_MPI)
380: distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
381: #endif
382: } else {
383: PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
384: hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
385: }
386: PetscCall(MatDenseRestoreArrayRead(B, (const PetscScalar **)&xx));
387: if (usesf) {
388: PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
389: PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
390: }
391: PetscCall(MatDenseRestoreArrayWrite(C, &yy));
392: #if defined(PETSC_H2OPUS_USE_GPU)
393: } else {
394: PetscBool ciscuda, biscuda;
396: /* If not of type seqdensecuda, convert on the fly (i.e. allocate GPU memory) */
397: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &biscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
398: if (!biscuda) PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
399: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &ciscuda, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
400: if (!ciscuda) {
401: C->assembled = PETSC_TRUE;
402: PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
403: }
404: PetscCall(MatDenseCUDAGetArrayRead(B, (const PetscScalar **)&xx));
405: PetscCall(MatDenseCUDAGetArrayWrite(C, &yy));
406: if (usesf) {
407: uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
408: uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
409: PetscCall(PetscSFBcastBegin(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
410: PetscCall(PetscSFBcastEnd(bsf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
411: } else {
412: uxx = xx;
413: uyy = yy;
414: }
415: PetscCall(PetscLogGpuTimeBegin());
416: if (size > 1) {
417: PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
418: PetscCheck(!transA || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
419: #if defined(H2OPUS_USE_MPI)
420: distributed_hgemv(/* transA ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, h2opus->handle);
421: #endif
422: } else {
423: PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
424: hgemv(transA ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, blda, 0.0, uyy, clda, B->cmap->N, handle);
425: }
426: PetscCall(PetscLogGpuTimeEnd());
427: PetscCall(MatDenseCUDARestoreArrayRead(B, (const PetscScalar **)&xx));
428: if (usesf) {
429: PetscCall(PetscSFReduceBegin(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
430: PetscCall(PetscSFReduceEnd(csf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
431: }
432: PetscCall(MatDenseCUDARestoreArrayWrite(C, &yy));
433: if (!biscuda) PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
434: if (!ciscuda) PetscCall(MatConvert(C, MATDENSE, MAT_INPLACE_MATRIX, &C));
435: #endif
436: }
437: { /* log flops */
438: double gops, time, perf, dev;
439: HLibProfile::getHgemvPerf(gops, time, perf, dev);
440: #if defined(PETSC_H2OPUS_USE_GPU)
441: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
442: else PetscCall(PetscLogGpuFlops(1e9 * gops));
443: #else
444: PetscCall(PetscLogFlops(1e9 * gops));
445: #endif
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode MatProductNumeric_H2OPUS(Mat C)
451: {
452: Mat_Product *product = C->product;
454: PetscFunctionBegin;
455: MatCheckProduct(C, 1);
456: switch (product->type) {
457: case MATPRODUCT_AB:
458: PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_FALSE, product->B, C));
459: break;
460: case MATPRODUCT_AtB:
461: PetscCall(MatMultNKernel_H2OPUS(product->A, PETSC_TRUE, product->B, C));
462: break;
463: default:
464: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
465: }
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: static PetscErrorCode MatProductSymbolic_H2OPUS(Mat C)
470: {
471: Mat_Product *product = C->product;
472: PetscBool cisdense;
473: Mat A, B;
475: PetscFunctionBegin;
476: MatCheckProduct(C, 1);
477: A = product->A;
478: B = product->B;
479: switch (product->type) {
480: case MATPRODUCT_AB:
481: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
482: PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
483: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
484: if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
485: PetscCall(MatSetUp(C));
486: break;
487: case MATPRODUCT_AtB:
488: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
489: PetscCall(MatSetBlockSizesFromMats(C, product->A, product->B));
490: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
491: if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)product->B)->type_name));
492: PetscCall(MatSetUp(C));
493: break;
494: default:
495: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProduct type %s is not supported", MatProductTypes[product->type]);
496: }
497: C->ops->productsymbolic = NULL;
498: C->ops->productnumeric = MatProductNumeric_H2OPUS;
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static PetscErrorCode MatProductSetFromOptions_H2OPUS(Mat C)
503: {
504: PetscFunctionBegin;
505: MatCheckProduct(C, 1);
506: if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_H2OPUS;
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: static PetscErrorCode MatMultKernel_H2OPUS(Mat A, Vec x, PetscScalar sy, Vec y, PetscBool trans)
511: {
512: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
513: #if defined(H2OPUS_USE_MPI)
514: h2opusHandle_t handle = h2opus->handle->handle;
515: #else
516: h2opusHandle_t handle = h2opus->handle;
517: #endif
518: PetscBool boundtocpu = PETSC_TRUE;
519: PetscInt n;
520: PetscScalar *xx, *yy, *uxx, *uyy;
521: PetscMPIInt size;
522: PetscBool usesf = (PetscBool)(h2opus->sf && !h2opus->nativemult);
524: PetscFunctionBegin;
525: HLibProfile::clear();
526: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
527: #if defined(PETSC_H2OPUS_USE_GPU)
528: boundtocpu = A->boundtocpu;
529: #endif
530: if (usesf) PetscCall(PetscSFGetGraph(h2opus->sf, NULL, &n, NULL, NULL));
531: else n = A->rmap->n;
532: if (boundtocpu) {
533: PetscCall(VecGetArrayRead(x, (const PetscScalar **)&xx));
534: if (sy == 0.0) {
535: PetscCall(VecGetArrayWrite(y, &yy));
536: } else {
537: PetscCall(VecGetArray(y, &yy));
538: }
539: if (usesf) {
540: uxx = MatH2OpusGetThrustPointer(*h2opus->xx);
541: uyy = MatH2OpusGetThrustPointer(*h2opus->yy);
543: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
544: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
545: if (sy != 0.0) {
546: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
547: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
548: }
549: } else {
550: uxx = xx;
551: uyy = yy;
552: }
553: if (size > 1) {
554: PetscCheck(h2opus->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
555: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
556: #if defined(H2OPUS_USE_MPI)
557: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix, uxx, n, sy, uyy, n, 1, h2opus->handle);
558: #endif
559: } else {
560: PetscCheck(h2opus->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
561: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix, uxx, n, sy, uyy, n, 1, handle);
562: }
563: PetscCall(VecRestoreArrayRead(x, (const PetscScalar **)&xx));
564: if (usesf) {
565: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
566: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
567: }
568: if (sy == 0.0) {
569: PetscCall(VecRestoreArrayWrite(y, &yy));
570: } else {
571: PetscCall(VecRestoreArray(y, &yy));
572: }
573: #if defined(PETSC_H2OPUS_USE_GPU)
574: } else {
575: PetscCall(VecCUDAGetArrayRead(x, (const PetscScalar **)&xx));
576: if (sy == 0.0) {
577: PetscCall(VecCUDAGetArrayWrite(y, &yy));
578: } else {
579: PetscCall(VecCUDAGetArray(y, &yy));
580: }
581: if (usesf) {
582: uxx = MatH2OpusGetThrustPointer(*h2opus->xx_gpu);
583: uyy = MatH2OpusGetThrustPointer(*h2opus->yy_gpu);
585: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
586: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, xx, uxx, MPI_REPLACE));
587: if (sy != 0.0) {
588: PetscCall(PetscSFBcastBegin(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
589: PetscCall(PetscSFBcastEnd(h2opus->sf, MPIU_SCALAR, yy, uyy, MPI_REPLACE));
590: }
591: } else {
592: uxx = xx;
593: uyy = yy;
594: }
595: PetscCall(PetscLogGpuTimeBegin());
596: if (size > 1) {
597: PetscCheck(h2opus->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed GPU matrix");
598: PetscCheck(!trans || A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatMultTranspose not yet coded in parallel");
599: #if defined(H2OPUS_USE_MPI)
600: distributed_hgemv(/*trans ? H2Opus_Trans : H2Opus_NoTrans, */ h2opus->s, *h2opus->dist_hmatrix_gpu, uxx, n, sy, uyy, n, 1, h2opus->handle);
601: #endif
602: } else {
603: PetscCheck(h2opus->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
604: hgemv(trans ? H2Opus_Trans : H2Opus_NoTrans, h2opus->s, *h2opus->hmatrix_gpu, uxx, n, sy, uyy, n, 1, handle);
605: }
606: PetscCall(PetscLogGpuTimeEnd());
607: PetscCall(VecCUDARestoreArrayRead(x, (const PetscScalar **)&xx));
608: if (usesf) {
609: PetscCall(PetscSFReduceBegin(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
610: PetscCall(PetscSFReduceEnd(h2opus->sf, MPIU_SCALAR, uyy, yy, MPI_REPLACE));
611: }
612: if (sy == 0.0) {
613: PetscCall(VecCUDARestoreArrayWrite(y, &yy));
614: } else {
615: PetscCall(VecCUDARestoreArray(y, &yy));
616: }
617: #endif
618: }
619: { /* log flops */
620: double gops, time, perf, dev;
621: HLibProfile::getHgemvPerf(gops, time, perf, dev);
622: #if defined(PETSC_H2OPUS_USE_GPU)
623: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
624: else PetscCall(PetscLogGpuFlops(1e9 * gops));
625: #else
626: PetscCall(PetscLogFlops(1e9 * gops));
627: #endif
628: }
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode MatMultTranspose_H2OPUS(Mat A, Vec x, Vec y)
633: {
634: PetscBool xiscuda, yiscuda;
636: PetscFunctionBegin;
637: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
638: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
639: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
640: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_TRUE));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode MatMult_H2OPUS(Mat A, Vec x, Vec y)
645: {
646: PetscBool xiscuda, yiscuda;
648: PetscFunctionBegin;
649: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
650: PetscCall(PetscObjectTypeCompareAny((PetscObject)y, &yiscuda, VECSEQCUDA, VECMPICUDA, ""));
651: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !yiscuda));
652: PetscCall(MatMultKernel_H2OPUS(A, x, 0.0, y, PETSC_FALSE));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode MatMultTransposeAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
657: {
658: PetscBool xiscuda, ziscuda;
660: PetscFunctionBegin;
661: PetscCall(VecCopy(y, z));
662: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
663: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
664: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
665: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_TRUE));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: static PetscErrorCode MatMultAdd_H2OPUS(Mat A, Vec x, Vec y, Vec z)
670: {
671: PetscBool xiscuda, ziscuda;
673: PetscFunctionBegin;
674: PetscCall(VecCopy(y, z));
675: PetscCall(PetscObjectTypeCompareAny((PetscObject)x, &xiscuda, VECSEQCUDA, VECMPICUDA, ""));
676: PetscCall(PetscObjectTypeCompareAny((PetscObject)z, &ziscuda, VECSEQCUDA, VECMPICUDA, ""));
677: PetscCall(MatH2OpusUpdateIfNeeded(A, !xiscuda || !ziscuda));
678: PetscCall(MatMultKernel_H2OPUS(A, x, 1.0, z, PETSC_FALSE));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode MatScale_H2OPUS(Mat A, PetscScalar s)
683: {
684: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
686: PetscFunctionBegin;
687: a->s *= s;
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode MatSetFromOptions_H2OPUS(Mat A, PetscOptionItems PetscOptionsObject)
692: {
693: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
695: PetscFunctionBegin;
696: PetscOptionsHeadBegin(PetscOptionsObject, "H2OPUS options");
697: PetscCall(PetscOptionsInt("-mat_h2opus_leafsize", "Leaf size of cluster tree", NULL, a->leafsize, &a->leafsize, NULL));
698: PetscCall(PetscOptionsReal("-mat_h2opus_eta", "Admissibility condition tolerance", NULL, a->eta, &a->eta, NULL));
699: PetscCall(PetscOptionsInt("-mat_h2opus_order", "Basis order for off-diagonal sampling when constructed from kernel", NULL, a->basisord, &a->basisord, NULL));
700: PetscCall(PetscOptionsInt("-mat_h2opus_maxrank", "Maximum rank when constructed from matvecs", NULL, a->max_rank, &a->max_rank, NULL));
701: PetscCall(PetscOptionsInt("-mat_h2opus_samples", "Maximum number of samples to be taken concurrently when constructing from matvecs", NULL, a->bs, &a->bs, NULL));
702: PetscCall(PetscOptionsInt("-mat_h2opus_normsamples", "Maximum number of samples to be when estimating norms", NULL, a->norm_max_samples, &a->norm_max_samples, NULL));
703: PetscCall(PetscOptionsReal("-mat_h2opus_rtol", "Relative tolerance for construction from sampling", NULL, a->rtol, &a->rtol, NULL));
704: PetscCall(PetscOptionsBool("-mat_h2opus_check", "Check error when constructing from sampling during MatAssemblyEnd()", NULL, a->check_construction, &a->check_construction, NULL));
705: PetscCall(PetscOptionsBool("-mat_h2opus_hara_verbose", "Verbose output from hara construction", NULL, a->hara_verbose, &a->hara_verbose, NULL));
706: PetscCall(PetscOptionsBool("-mat_h2opus_resize", "Resize after compression", NULL, a->resize, &a->resize, NULL));
707: PetscOptionsHeadEnd();
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat, PetscInt, const PetscReal[], PetscBool, MatH2OpusKernelFn *, void *);
713: static PetscErrorCode MatH2OpusInferCoordinates_Private(Mat A)
714: {
715: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
716: Vec c;
717: PetscInt spacedim;
718: const PetscScalar *coords;
720: PetscFunctionBegin;
721: if (a->ptcloud) PetscFunctionReturn(PETSC_SUCCESS);
722: PetscCall(PetscObjectQuery((PetscObject)A, "__math2opus_coords", (PetscObject *)&c));
723: if (!c && a->sampler) {
724: Mat S = a->sampler->GetSamplingMat();
726: PetscCall(PetscObjectQuery((PetscObject)S, "__math2opus_coords", (PetscObject *)&c));
727: }
728: if (!c) {
729: PetscCall(MatH2OpusSetCoords_H2OPUS(A, -1, NULL, PETSC_FALSE, NULL, NULL));
730: } else {
731: PetscCall(VecGetArrayRead(c, &coords));
732: PetscCall(VecGetBlockSize(c, &spacedim));
733: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, PETSC_FALSE, NULL, NULL));
734: PetscCall(VecRestoreArrayRead(c, &coords));
735: }
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode MatSetUpMultiply_H2OPUS(Mat A)
740: {
741: MPI_Comm comm;
742: PetscMPIInt size;
743: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
744: PetscInt n = 0, *idx = NULL;
745: int *iidx = NULL;
746: PetscCopyMode own;
747: PetscBool rid;
749: PetscFunctionBegin;
750: if (a->multsetup) PetscFunctionReturn(PETSC_SUCCESS);
751: if (a->sf) { /* MatDuplicate_H2OPUS takes reference to the SF */
752: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
753: #if defined(PETSC_H2OPUS_USE_GPU)
754: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
755: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
756: a->xxs_gpu = 1;
757: a->yys_gpu = 1;
758: #endif
759: a->xx = new thrust::host_vector<PetscScalar>(n);
760: a->yy = new thrust::host_vector<PetscScalar>(n);
761: a->xxs = 1;
762: a->yys = 1;
763: } else {
764: IS is;
765: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
766: PetscCallMPI(MPI_Comm_size(comm, &size));
767: if (!a->h2opus_indexmap) {
768: if (size > 1) {
769: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
770: #if defined(H2OPUS_USE_MPI)
771: iidx = MatH2OpusGetThrustPointer(a->dist_hmatrix->basis_tree.basis_branch.index_map);
772: n = a->dist_hmatrix->basis_tree.basis_branch.index_map.size();
773: #endif
774: } else {
775: iidx = MatH2OpusGetThrustPointer(a->hmatrix->u_basis_tree.index_map);
776: n = a->hmatrix->u_basis_tree.index_map.size();
777: }
779: if (PetscDefined(USE_64BIT_INDICES)) {
780: PetscInt i;
782: own = PETSC_OWN_POINTER;
783: PetscCall(PetscMalloc1(n, &idx));
784: for (i = 0; i < n; i++) idx[i] = iidx[i];
785: } else {
786: own = PETSC_COPY_VALUES;
787: idx = (PetscInt *)iidx;
788: }
789: PetscCall(ISCreateGeneral(comm, n, idx, own, &is));
790: PetscCall(ISSetPermutation(is));
791: PetscCall(ISViewFromOptions(is, (PetscObject)A, "-mat_h2opus_indexmap_view"));
792: a->h2opus_indexmap = is;
793: }
794: PetscCall(ISGetLocalSize(a->h2opus_indexmap, &n));
795: PetscCall(ISGetIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
796: rid = (PetscBool)(n == A->rmap->n);
797: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &rid, 1, MPI_C_BOOL, MPI_LAND, comm));
798: if (rid) PetscCall(ISIdentity(a->h2opus_indexmap, &rid));
799: if (!rid) {
800: if (size > 1) { /* Parallel distribution may be different, save it here for fast path in MatMult (see MatH2OpusSetNativeMult) */
801: PetscCall(PetscLayoutCreate(comm, &a->h2opus_rmap));
802: PetscCall(PetscLayoutSetLocalSize(a->h2opus_rmap, n));
803: PetscCall(PetscLayoutSetUp(a->h2opus_rmap));
804: PetscCall(PetscLayoutReference(a->h2opus_rmap, &a->h2opus_cmap));
805: }
806: PetscCall(PetscSFCreate(comm, &a->sf));
807: PetscCall(PetscSFSetGraphLayout(a->sf, A->rmap, n, NULL, PETSC_OWN_POINTER, idx));
808: PetscCall(PetscSFSetUp(a->sf));
809: PetscCall(PetscSFViewFromOptions(a->sf, (PetscObject)A, "-mat_h2opus_sf_view"));
810: #if defined(PETSC_H2OPUS_USE_GPU)
811: a->xx_gpu = new thrust::device_vector<PetscScalar>(n);
812: a->yy_gpu = new thrust::device_vector<PetscScalar>(n);
813: a->xxs_gpu = 1;
814: a->yys_gpu = 1;
815: #endif
816: a->xx = new thrust::host_vector<PetscScalar>(n);
817: a->yy = new thrust::host_vector<PetscScalar>(n);
818: a->xxs = 1;
819: a->yys = 1;
820: }
821: PetscCall(ISRestoreIndices(a->h2opus_indexmap, (const PetscInt **)&idx));
822: }
823: a->multsetup = PETSC_TRUE;
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: static PetscErrorCode MatAssemblyEnd_H2OPUS(Mat A, MatAssemblyType assemblytype)
828: {
829: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
830: #if defined(H2OPUS_USE_MPI)
831: h2opusHandle_t handle = a->handle->handle;
832: #else
833: h2opusHandle_t handle = a->handle;
834: #endif
835: PetscBool kernel = PETSC_FALSE;
836: PetscBool boundtocpu = PETSC_TRUE;
837: PetscBool samplingdone = PETSC_FALSE;
838: MPI_Comm comm;
839: PetscMPIInt size;
841: PetscFunctionBegin;
842: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
843: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
844: PetscCheck(A->rmap->N == A->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
846: /* XXX */
847: a->leafsize = PetscMin(a->leafsize, PetscMin(A->rmap->N, A->cmap->N));
849: PetscCallMPI(MPI_Comm_size(comm, &size));
850: /* TODO REUSABILITY of geometric construction */
851: delete a->hmatrix;
852: delete a->dist_hmatrix;
853: #if defined(PETSC_H2OPUS_USE_GPU)
854: delete a->hmatrix_gpu;
855: delete a->dist_hmatrix_gpu;
856: #endif
857: a->orthogonal = PETSC_FALSE;
859: /* TODO: other? */
860: H2OpusBoxCenterAdmissibility adm(a->eta);
862: PetscCall(PetscLogEventBegin(MAT_H2Opus_Build, A, 0, 0, 0));
863: if (size > 1) {
864: #if defined(H2OPUS_USE_MPI)
865: a->dist_hmatrix = new DistributedHMatrix(A->rmap->n /* ,A->symmetric */);
866: #else
867: a->dist_hmatrix = NULL;
868: #endif
869: } else a->hmatrix = new HMatrix(A->rmap->n, A->symmetric == PETSC_BOOL3_TRUE);
870: PetscCall(MatH2OpusInferCoordinates_Private(A));
871: PetscCheck(a->ptcloud, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing pointcloud");
872: if (a->kernel) {
873: BoxEntryGen<PetscScalar, H2OPUS_HWTYPE_CPU, PetscFunctionGenerator<PetscScalar>> entry_gen(*a->kernel);
874: if (size > 1) {
875: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
876: #if defined(H2OPUS_USE_MPI)
877: buildDistributedHMatrix(*a->dist_hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord, a->handle);
878: #endif
879: } else {
880: buildHMatrix(*a->hmatrix, a->ptcloud, adm, entry_gen, a->leafsize, a->basisord);
881: }
882: kernel = PETSC_TRUE;
883: } else {
884: PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Construction from sampling not supported in parallel");
885: buildHMatrixStructure(*a->hmatrix, a->ptcloud, a->leafsize, adm);
886: }
887: PetscCall(MatSetUpMultiply_H2OPUS(A));
889: #if defined(PETSC_H2OPUS_USE_GPU)
890: boundtocpu = A->boundtocpu;
891: if (!boundtocpu) {
892: if (size > 1) {
893: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing distributed CPU matrix");
894: #if defined(H2OPUS_USE_MPI)
895: a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
896: #endif
897: } else {
898: a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
899: }
900: }
901: #endif
902: if (size == 1) {
903: if (!kernel && a->sampler && a->sampler->GetSamplingMat()) {
904: PetscReal Anorm;
905: bool verbose;
907: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_hara_verbose", &a->hara_verbose, NULL));
908: verbose = a->hara_verbose;
909: PetscCall(MatNormApproximate(a->sampler->GetSamplingMat(), NORM_2, a->norm_max_samples, &Anorm));
910: 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));
911: if (a->sf && !a->nativemult) a->sampler->SetIndexMap(a->hmatrix->u_basis_tree.index_map.size(), a->hmatrix->u_basis_tree.index_map.data());
912: a->sampler->SetStream(handle->getMainStream());
913: if (boundtocpu) {
914: a->sampler->SetGPUSampling(false);
915: hara(a->sampler, *a->hmatrix, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
916: #if defined(PETSC_H2OPUS_USE_GPU)
917: } else {
918: a->sampler->SetGPUSampling(true);
919: hara(a->sampler, *a->hmatrix_gpu, a->max_rank, 10 /* TODO */, a->rtol * Anorm, a->bs, handle, verbose);
920: #endif
921: }
922: samplingdone = PETSC_TRUE;
923: }
924: }
925: #if defined(PETSC_H2OPUS_USE_GPU)
926: if (!boundtocpu) {
927: delete a->hmatrix;
928: delete a->dist_hmatrix;
929: a->hmatrix = NULL;
930: a->dist_hmatrix = NULL;
931: }
932: A->offloadmask = boundtocpu ? PETSC_OFFLOAD_CPU : PETSC_OFFLOAD_GPU;
933: #endif
934: PetscCall(PetscLogEventEnd(MAT_H2Opus_Build, A, 0, 0, 0));
936: if (!a->s) a->s = 1.0;
937: A->assembled = PETSC_TRUE;
939: if (samplingdone) {
940: PetscBool check = a->check_construction;
941: PetscBool checke = PETSC_FALSE;
943: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check", &check, NULL));
944: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_check_explicit", &checke, NULL));
945: if (check) {
946: Mat E, Ae;
947: PetscReal n1, ni, n2;
948: PetscReal n1A, niA, n2A;
949: PetscErrorCodeFn *normfunc;
951: Ae = a->sampler->GetSamplingMat();
952: PetscCall(MatConvert(A, MATSHELL, MAT_INITIAL_MATRIX, &E));
953: PetscCall(MatShellSetOperation(E, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
954: PetscCall(MatAXPY(E, -1.0, Ae, DIFFERENT_NONZERO_PATTERN));
955: PetscCall(MatNorm(E, NORM_1, &n1));
956: PetscCall(MatNorm(E, NORM_INFINITY, &ni));
957: PetscCall(MatNorm(E, NORM_2, &n2));
958: if (checke) {
959: Mat eA, eE, eAe;
961: PetscCall(MatComputeOperator(A, MATAIJ, &eA));
962: PetscCall(MatComputeOperator(E, MATAIJ, &eE));
963: PetscCall(MatComputeOperator(Ae, MATAIJ, &eAe));
964: PetscCall(MatFilter(eA, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
965: PetscCall(MatFilter(eE, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
966: PetscCall(MatFilter(eAe, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
967: PetscCall(PetscObjectSetName((PetscObject)eA, "H2Mat"));
968: PetscCall(MatView(eA, NULL));
969: PetscCall(PetscObjectSetName((PetscObject)eAe, "S"));
970: PetscCall(MatView(eAe, NULL));
971: PetscCall(PetscObjectSetName((PetscObject)eE, "H2Mat - S"));
972: PetscCall(MatView(eE, NULL));
973: PetscCall(MatDestroy(&eA));
974: PetscCall(MatDestroy(&eE));
975: PetscCall(MatDestroy(&eAe));
976: }
978: PetscCall(MatGetOperation(Ae, MATOP_NORM, &normfunc));
979: PetscCall(MatSetOperation(Ae, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
980: PetscCall(MatNorm(Ae, NORM_1, &n1A));
981: PetscCall(MatNorm(Ae, NORM_INFINITY, &niA));
982: PetscCall(MatNorm(Ae, NORM_2, &n2A));
983: n1A = PetscMax(n1A, PETSC_SMALL);
984: n2A = PetscMax(n2A, PETSC_SMALL);
985: niA = PetscMax(niA, PETSC_SMALL);
986: PetscCall(MatSetOperation(Ae, MATOP_NORM, normfunc));
987: 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)));
988: PetscCall(MatDestroy(&E));
989: }
990: a->sampler->SetSamplingMat(NULL);
991: }
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: static PetscErrorCode MatZeroEntries_H2OPUS(Mat A)
996: {
997: PetscMPIInt size;
998: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1000: PetscFunctionBegin;
1001: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1002: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet supported");
1003: a->hmatrix->clearData();
1004: #if defined(PETSC_H2OPUS_USE_GPU)
1005: if (a->hmatrix_gpu) a->hmatrix_gpu->clearData();
1006: #endif
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: static PetscErrorCode MatDuplicate_H2OPUS(Mat B, MatDuplicateOption op, Mat *nA)
1011: {
1012: Mat A;
1013: Mat_H2OPUS *a, *b = (Mat_H2OPUS *)B->data;
1014: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1015: MPI_Comm comm;
1017: PetscFunctionBegin;
1018: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1019: PetscCall(MatCreate(comm, &A));
1020: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1021: PetscCall(MatSetType(A, MATH2OPUS));
1022: PetscCall(MatPropagateSymmetryOptions(B, A));
1023: a = (Mat_H2OPUS *)A->data;
1025: a->eta = b->eta;
1026: a->leafsize = b->leafsize;
1027: a->basisord = b->basisord;
1028: a->max_rank = b->max_rank;
1029: a->bs = b->bs;
1030: a->rtol = b->rtol;
1031: a->norm_max_samples = b->norm_max_samples;
1032: if (op == MAT_COPY_VALUES) a->s = b->s;
1034: a->ptcloud = new PetscPointCloud<PetscReal>(*b->ptcloud);
1035: if (op == MAT_COPY_VALUES && b->kernel) a->kernel = new PetscFunctionGenerator<PetscScalar>(*b->kernel);
1037: #if defined(H2OPUS_USE_MPI)
1038: if (b->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*b->dist_hmatrix);
1039: #if defined(PETSC_H2OPUS_USE_GPU)
1040: if (b->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*b->dist_hmatrix_gpu);
1041: #endif
1042: #endif
1043: if (b->hmatrix) {
1044: a->hmatrix = new HMatrix(*b->hmatrix);
1045: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix->clearData();
1046: }
1047: #if defined(PETSC_H2OPUS_USE_GPU)
1048: if (b->hmatrix_gpu) {
1049: a->hmatrix_gpu = new HMatrix_GPU(*b->hmatrix_gpu);
1050: if (op == MAT_DO_NOT_COPY_VALUES) a->hmatrix_gpu->clearData();
1051: }
1052: #endif
1053: if (b->sf) {
1054: PetscCall(PetscObjectReference((PetscObject)b->sf));
1055: a->sf = b->sf;
1056: }
1057: if (b->h2opus_indexmap) {
1058: PetscCall(PetscObjectReference((PetscObject)b->h2opus_indexmap));
1059: a->h2opus_indexmap = b->h2opus_indexmap;
1060: }
1062: PetscCall(MatSetUp(A));
1063: PetscCall(MatSetUpMultiply_H2OPUS(A));
1064: if (op == MAT_COPY_VALUES) {
1065: A->assembled = PETSC_TRUE;
1066: a->orthogonal = b->orthogonal;
1067: #if defined(PETSC_H2OPUS_USE_GPU)
1068: A->offloadmask = B->offloadmask;
1069: #endif
1070: }
1071: #if defined(PETSC_H2OPUS_USE_GPU)
1072: iscpu = B->boundtocpu;
1073: #endif
1074: PetscCall(MatBindToCPU(A, iscpu));
1076: *nA = A;
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: static PetscErrorCode MatView_H2OPUS(Mat A, PetscViewer view)
1081: {
1082: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1083: PetscBool isascii, vieweps;
1084: PetscMPIInt size;
1085: PetscViewerFormat format;
1087: PetscFunctionBegin;
1088: PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
1089: PetscCall(PetscViewerGetFormat(view, &format));
1090: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1091: if (isascii) {
1092: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1093: if (size == 1) {
1094: FILE *fp;
1095: PetscCall(PetscViewerASCIIGetPointer(view, &fp));
1096: dumpHMatrix(*h2opus->hmatrix, 6, fp);
1097: }
1098: } else {
1099: PetscCall(PetscViewerASCIIPrintf(view, " H-Matrix constructed from %s\n", h2opus->kernel ? "Kernel" : "Mat"));
1100: PetscCall(PetscViewerASCIIPrintf(view, " PointCloud dim %" PetscInt_FMT "\n", h2opus->ptcloud ? h2opus->ptcloud->getDimension() : 0));
1101: PetscCall(PetscViewerASCIIPrintf(view, " Admissibility parameters: leaf size %" PetscInt_FMT ", eta %g\n", h2opus->leafsize, (double)h2opus->eta));
1102: if (!h2opus->kernel) {
1103: PetscCall(PetscViewerASCIIPrintf(view, " Sampling parameters: max_rank %" PetscInt_FMT ", samples %" PetscInt_FMT ", tolerance %g\n", h2opus->max_rank, h2opus->bs, (double)h2opus->rtol));
1104: } else {
1105: PetscCall(PetscViewerASCIIPrintf(view, " Off-diagonal blocks approximation order %" PetscInt_FMT "\n", h2opus->basisord));
1106: }
1107: PetscCall(PetscViewerASCIIPrintf(view, " Number of samples for norms %" PetscInt_FMT "\n", h2opus->norm_max_samples));
1108: if (size == 1) {
1109: double dense_mem_cpu = h2opus->hmatrix ? h2opus->hmatrix->getDenseMemoryUsage() : 0;
1110: double low_rank_cpu = h2opus->hmatrix ? h2opus->hmatrix->getLowRankMemoryUsage() : 0;
1111: #if defined(PETSC_HAVE_CUDA)
1112: double dense_mem_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getDenseMemoryUsage() : 0;
1113: double low_rank_gpu = h2opus->hmatrix_gpu ? h2opus->hmatrix_gpu->getLowRankMemoryUsage() : 0;
1114: #endif
1115: 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));
1116: #if defined(PETSC_HAVE_CUDA)
1117: 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));
1118: #endif
1119: } else {
1120: #if defined(PETSC_HAVE_CUDA)
1121: double matrix_mem[4] = {0., 0., 0., 0.};
1122: PetscMPIInt rsize = 4;
1123: #else
1124: double matrix_mem[2] = {0., 0.};
1125: PetscMPIInt rsize = 2;
1126: #endif
1127: #if defined(H2OPUS_USE_MPI)
1128: matrix_mem[0] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalDenseMemoryUsage() : 0;
1129: matrix_mem[1] = h2opus->dist_hmatrix ? h2opus->dist_hmatrix->getLocalLowRankMemoryUsage() : 0;
1130: #if defined(PETSC_HAVE_CUDA)
1131: matrix_mem[2] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalDenseMemoryUsage() : 0;
1132: matrix_mem[3] = h2opus->dist_hmatrix_gpu ? h2opus->dist_hmatrix_gpu->getLocalLowRankMemoryUsage() : 0;
1133: #endif
1134: #endif
1135: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, matrix_mem, rsize, MPI_DOUBLE_PRECISION, MPI_SUM, PetscObjectComm((PetscObject)A)));
1136: 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]));
1137: #if defined(PETSC_HAVE_CUDA)
1138: 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]));
1139: #endif
1140: }
1141: }
1142: }
1143: vieweps = PETSC_FALSE;
1144: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps", &vieweps, NULL));
1145: if (vieweps) {
1146: char filename[256];
1147: const char *name;
1149: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1150: PetscCall(PetscSNPrintf(filename, sizeof(filename), "%s_structure.eps", name));
1151: PetscCall(PetscOptionsGetString(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_h2opus_vieweps_filename", filename, sizeof(filename), NULL));
1152: outputEps(*h2opus->hmatrix, filename);
1153: }
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: static PetscErrorCode MatH2OpusSetCoords_H2OPUS(Mat A, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, MatH2OpusKernelFn *kernel, void *kernelctx)
1158: {
1159: Mat_H2OPUS *h2opus = (Mat_H2OPUS *)A->data;
1160: PetscReal *gcoords;
1161: PetscInt N;
1162: MPI_Comm comm;
1163: PetscMPIInt size;
1164: PetscBool cong;
1166: PetscFunctionBegin;
1167: PetscCall(PetscLayoutSetUp(A->rmap));
1168: PetscCall(PetscLayoutSetUp(A->cmap));
1169: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1170: PetscCall(MatHasCongruentLayouts(A, &cong));
1171: PetscCheck(cong, comm, PETSC_ERR_SUP, "Only for square matrices with congruent layouts");
1172: N = A->rmap->N;
1173: PetscCallMPI(MPI_Comm_size(comm, &size));
1174: if (spacedim > 0 && size > 1 && cdist) {
1175: PetscSF sf;
1176: MPI_Datatype dtype;
1178: PetscCallMPI(MPI_Type_contiguous(spacedim, MPIU_REAL, &dtype));
1179: PetscCallMPI(MPI_Type_commit(&dtype));
1181: PetscCall(PetscSFCreate(comm, &sf));
1182: PetscCall(PetscSFSetGraphWithPattern(sf, A->rmap, PETSCSF_PATTERN_ALLGATHER));
1183: PetscCall(PetscMalloc1(spacedim * N, &gcoords));
1184: PetscCall(PetscSFBcastBegin(sf, dtype, coords, gcoords, MPI_REPLACE));
1185: PetscCall(PetscSFBcastEnd(sf, dtype, coords, gcoords, MPI_REPLACE));
1186: PetscCall(PetscSFDestroy(&sf));
1187: PetscCallMPI(MPI_Type_free(&dtype));
1188: } else gcoords = (PetscReal *)coords;
1190: delete h2opus->ptcloud;
1191: delete h2opus->kernel;
1192: h2opus->ptcloud = new PetscPointCloud<PetscReal>(spacedim, N, gcoords);
1193: if (kernel) h2opus->kernel = new PetscFunctionGenerator<PetscScalar>(kernel, spacedim, kernelctx);
1194: if (gcoords != coords) PetscCall(PetscFree(gcoords));
1195: A->preallocated = PETSC_TRUE;
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: #if defined(PETSC_H2OPUS_USE_GPU)
1200: static PetscErrorCode MatBindToCPU_H2OPUS(Mat A, PetscBool flg)
1201: {
1202: PetscMPIInt size;
1203: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1205: PetscFunctionBegin;
1206: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1207: if (flg && A->offloadmask == PETSC_OFFLOAD_GPU) {
1208: if (size > 1) {
1209: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1210: #if defined(H2OPUS_USE_MPI)
1211: if (!a->dist_hmatrix) a->dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix_gpu);
1212: else *a->dist_hmatrix = *a->dist_hmatrix_gpu;
1213: #endif
1214: } else {
1215: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1216: if (!a->hmatrix) a->hmatrix = new HMatrix(*a->hmatrix_gpu);
1217: else *a->hmatrix = *a->hmatrix_gpu;
1218: }
1219: delete a->hmatrix_gpu;
1220: delete a->dist_hmatrix_gpu;
1221: a->hmatrix_gpu = NULL;
1222: a->dist_hmatrix_gpu = NULL;
1223: A->offloadmask = PETSC_OFFLOAD_CPU;
1224: } else if (!flg && A->offloadmask == PETSC_OFFLOAD_CPU) {
1225: if (size > 1) {
1226: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1227: #if defined(H2OPUS_USE_MPI)
1228: if (!a->dist_hmatrix_gpu) a->dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix);
1229: else *a->dist_hmatrix_gpu = *a->dist_hmatrix;
1230: #endif
1231: } else {
1232: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1233: if (!a->hmatrix_gpu) a->hmatrix_gpu = new HMatrix_GPU(*a->hmatrix);
1234: else *a->hmatrix_gpu = *a->hmatrix;
1235: }
1236: delete a->hmatrix;
1237: delete a->dist_hmatrix;
1238: a->hmatrix = NULL;
1239: a->dist_hmatrix = NULL;
1240: A->offloadmask = PETSC_OFFLOAD_GPU;
1241: }
1242: PetscCall(PetscFree(A->defaultvectype));
1243: if (!flg) {
1244: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1245: } else {
1246: PetscCall(PetscStrallocpy(VECSTANDARD, &A->defaultvectype));
1247: }
1248: A->boundtocpu = flg;
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }
1251: #endif
1253: /*MC
1254: MATH2OPUS = "h2opus" - A matrix type for hierarchical matrices using the H2Opus package {cite}`zampinibouakaramturkiyyahkniokeyes2022`.
1256: Options Database Key:
1257: . -mat_type h2opus - matrix type to "h2opus"
1259: Level: beginner
1261: Notes:
1262: H2Opus implements hierarchical matrices in the $H^2$ flavor. It supports CPU or NVIDIA GPUs.
1264: For CPU only builds, use `./configure --download-h2opus --download-thrust` to install PETSc to use H2Opus.
1265: In order to run on NVIDIA GPUs, use `./configure --download-h2opus --download-magma --download-kblas`.
1267: .seealso: [](ch_matrices), `Mat`, `MATH2OPUS`, `MATHTOOL`, `MATDENSE`, `MatCreateH2OpusFromKernel()`, `MatCreateH2OpusFromMat()`
1268: M*/
1269: PETSC_EXTERN PetscErrorCode MatCreate_H2OPUS(Mat A)
1270: {
1271: Mat_H2OPUS *a;
1272: PetscMPIInt size;
1274: PetscFunctionBegin;
1275: #if defined(PETSC_H2OPUS_USE_GPU)
1276: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
1277: #endif
1278: PetscCall(PetscNew(&a));
1279: A->data = (void *)a;
1281: a->eta = 0.9;
1282: a->leafsize = 32;
1283: a->basisord = 4;
1284: a->max_rank = 64;
1285: a->bs = 32;
1286: a->rtol = 1.e-4;
1287: a->s = 1.0;
1288: a->norm_max_samples = 10;
1289: a->resize = PETSC_TRUE; /* reallocate after compression */
1290: #if defined(H2OPUS_USE_MPI)
1291: h2opusCreateDistributedHandleComm(&a->handle, PetscObjectComm((PetscObject)A));
1292: #else
1293: h2opusCreateHandle(&a->handle);
1294: #endif
1295: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1296: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATH2OPUS));
1297: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
1299: A->ops->destroy = MatDestroy_H2OPUS;
1300: A->ops->view = MatView_H2OPUS;
1301: A->ops->assemblyend = MatAssemblyEnd_H2OPUS;
1302: A->ops->mult = MatMult_H2OPUS;
1303: A->ops->multtranspose = MatMultTranspose_H2OPUS;
1304: A->ops->multadd = MatMultAdd_H2OPUS;
1305: A->ops->multtransposeadd = MatMultTransposeAdd_H2OPUS;
1306: A->ops->scale = MatScale_H2OPUS;
1307: A->ops->duplicate = MatDuplicate_H2OPUS;
1308: A->ops->setfromoptions = MatSetFromOptions_H2OPUS;
1309: A->ops->norm = MatNorm_H2OPUS;
1310: A->ops->zeroentries = MatZeroEntries_H2OPUS;
1311: #if defined(PETSC_H2OPUS_USE_GPU)
1312: A->ops->bindtocpu = MatBindToCPU_H2OPUS;
1313: #endif
1315: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdense_C", MatProductSetFromOptions_H2OPUS));
1316: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_seqdensecuda_C", MatProductSetFromOptions_H2OPUS));
1317: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidense_C", MatProductSetFromOptions_H2OPUS));
1318: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_h2opus_mpidensecuda_C", MatProductSetFromOptions_H2OPUS));
1319: #if defined(PETSC_H2OPUS_USE_GPU)
1320: PetscCall(PetscFree(A->defaultvectype));
1321: PetscCall(PetscStrallocpy(VECCUDA, &A->defaultvectype));
1322: #endif
1323: PetscFunctionReturn(PETSC_SUCCESS);
1324: }
1326: /*@
1327: MatH2OpusOrthogonalize - Orthogonalize the basis tree of a hierarchical matrix.
1329: Input Parameter:
1330: . A - the matrix
1332: Level: intermediate
1334: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`
1335: @*/
1336: PetscErrorCode MatH2OpusOrthogonalize(Mat A)
1337: {
1338: PetscBool ish2opus;
1339: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1340: PetscMPIInt size;
1341: PetscBool boundtocpu = PETSC_TRUE;
1343: PetscFunctionBegin;
1346: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1347: if (!ish2opus) PetscFunctionReturn(PETSC_SUCCESS);
1348: if (a->orthogonal) PetscFunctionReturn(PETSC_SUCCESS);
1349: HLibProfile::clear();
1350: PetscCall(PetscLogEventBegin(MAT_H2Opus_Orthog, A, 0, 0, 0));
1351: #if defined(PETSC_H2OPUS_USE_GPU)
1352: boundtocpu = A->boundtocpu;
1353: #endif
1354: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1355: if (size > 1) {
1356: if (boundtocpu) {
1357: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1358: #if defined(H2OPUS_USE_MPI)
1359: distributed_horthog(*a->dist_hmatrix, a->handle);
1360: #endif
1361: #if defined(PETSC_H2OPUS_USE_GPU)
1362: A->offloadmask = PETSC_OFFLOAD_CPU;
1363: } else {
1364: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1365: PetscCall(PetscLogGpuTimeBegin());
1366: #if defined(H2OPUS_USE_MPI)
1367: distributed_horthog(*a->dist_hmatrix_gpu, a->handle);
1368: #endif
1369: PetscCall(PetscLogGpuTimeEnd());
1370: #endif
1371: }
1372: } else {
1373: #if defined(H2OPUS_USE_MPI)
1374: h2opusHandle_t handle = a->handle->handle;
1375: #else
1376: h2opusHandle_t handle = a->handle;
1377: #endif
1378: if (boundtocpu) {
1379: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1380: horthog(*a->hmatrix, handle);
1381: #if defined(PETSC_H2OPUS_USE_GPU)
1382: A->offloadmask = PETSC_OFFLOAD_CPU;
1383: } else {
1384: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1385: PetscCall(PetscLogGpuTimeBegin());
1386: horthog(*a->hmatrix_gpu, handle);
1387: PetscCall(PetscLogGpuTimeEnd());
1388: #endif
1389: }
1390: }
1391: a->orthogonal = PETSC_TRUE;
1392: { /* log flops */
1393: double gops, time, perf, dev;
1394: HLibProfile::getHorthogPerf(gops, time, perf, dev);
1395: #if defined(PETSC_H2OPUS_USE_GPU)
1396: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1397: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1398: #else
1399: PetscCall(PetscLogFlops(1e9 * gops));
1400: #endif
1401: }
1402: PetscCall(PetscLogEventEnd(MAT_H2Opus_Orthog, A, 0, 0, 0));
1403: PetscFunctionReturn(PETSC_SUCCESS);
1404: }
1406: /*@
1407: MatH2OpusCompress - Compress a hierarchical matrix.
1409: Input Parameters:
1410: + A - the matrix
1411: - tol - the absolute truncation threshold
1413: Level: intermediate
1415: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusOrthogonalize()`
1416: @*/
1417: PetscErrorCode MatH2OpusCompress(Mat A, PetscReal tol)
1418: {
1419: PetscBool ish2opus;
1420: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1421: PetscMPIInt size;
1422: PetscBool boundtocpu = PETSC_TRUE;
1424: PetscFunctionBegin;
1428: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1429: if (!ish2opus || tol <= 0.0) PetscFunctionReturn(PETSC_SUCCESS);
1430: PetscCall(MatH2OpusOrthogonalize(A));
1431: HLibProfile::clear();
1432: PetscCall(PetscLogEventBegin(MAT_H2Opus_Compress, A, 0, 0, 0));
1433: #if defined(PETSC_H2OPUS_USE_GPU)
1434: boundtocpu = A->boundtocpu;
1435: #endif
1436: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1437: if (size > 1) {
1438: if (boundtocpu) {
1439: PetscCheck(a->dist_hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1440: #if defined(H2OPUS_USE_MPI)
1441: distributed_hcompress(*a->dist_hmatrix, tol, a->handle);
1442: if (a->resize) {
1443: DistributedHMatrix *dist_hmatrix = new DistributedHMatrix(*a->dist_hmatrix);
1444: delete a->dist_hmatrix;
1445: a->dist_hmatrix = dist_hmatrix;
1446: }
1447: #endif
1448: #if defined(PETSC_H2OPUS_USE_GPU)
1449: A->offloadmask = PETSC_OFFLOAD_CPU;
1450: } else {
1451: PetscCheck(a->dist_hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1452: PetscCall(PetscLogGpuTimeBegin());
1453: #if defined(H2OPUS_USE_MPI)
1454: distributed_hcompress(*a->dist_hmatrix_gpu, tol, a->handle);
1456: if (a->resize) {
1457: DistributedHMatrix_GPU *dist_hmatrix_gpu = new DistributedHMatrix_GPU(*a->dist_hmatrix_gpu);
1458: delete a->dist_hmatrix_gpu;
1459: a->dist_hmatrix_gpu = dist_hmatrix_gpu;
1460: }
1461: #endif
1462: PetscCall(PetscLogGpuTimeEnd());
1463: #endif
1464: }
1465: } else {
1466: #if defined(H2OPUS_USE_MPI)
1467: h2opusHandle_t handle = a->handle->handle;
1468: #else
1469: h2opusHandle_t handle = a->handle;
1470: #endif
1471: if (boundtocpu) {
1472: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1473: hcompress(*a->hmatrix, tol, handle);
1475: if (a->resize) {
1476: HMatrix *hmatrix = new HMatrix(*a->hmatrix);
1477: delete a->hmatrix;
1478: a->hmatrix = hmatrix;
1479: }
1480: #if defined(PETSC_H2OPUS_USE_GPU)
1481: A->offloadmask = PETSC_OFFLOAD_CPU;
1482: } else {
1483: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1484: PetscCall(PetscLogGpuTimeBegin());
1485: hcompress(*a->hmatrix_gpu, tol, handle);
1486: PetscCall(PetscLogGpuTimeEnd());
1488: if (a->resize) {
1489: HMatrix_GPU *hmatrix_gpu = new HMatrix_GPU(*a->hmatrix_gpu);
1490: delete a->hmatrix_gpu;
1491: a->hmatrix_gpu = hmatrix_gpu;
1492: }
1493: #endif
1494: }
1495: }
1496: { /* log flops */
1497: double gops, time, perf, dev;
1498: HLibProfile::getHcompressPerf(gops, time, perf, dev);
1499: #if defined(PETSC_H2OPUS_USE_GPU)
1500: if (boundtocpu) PetscCall(PetscLogFlops(1e9 * gops));
1501: else PetscCall(PetscLogGpuFlops(1e9 * gops));
1502: #else
1503: PetscCall(PetscLogFlops(1e9 * gops));
1504: #endif
1505: }
1506: PetscCall(PetscLogEventEnd(MAT_H2Opus_Compress, A, 0, 0, 0));
1507: PetscFunctionReturn(PETSC_SUCCESS);
1508: }
1510: /*@
1511: MatH2OpusSetSamplingMat - Set a matrix to be sampled from matrix-vector products on another matrix to construct a hierarchical matrix.
1513: Input Parameters:
1514: + A - the hierarchical matrix
1515: . B - the matrix to be sampled
1516: . bs - maximum number of samples to be taken concurrently
1517: - tol - relative tolerance for construction
1519: Level: intermediate
1521: Notes:
1522: You need to call `MatAssemblyBegin()` and `MatAssemblyEnd()` to update the hierarchical matrix.
1524: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`
1525: @*/
1526: PetscErrorCode MatH2OpusSetSamplingMat(Mat A, Mat B, PetscInt bs, PetscReal tol)
1527: {
1528: PetscBool ish2opus;
1530: PetscFunctionBegin;
1536: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1537: if (ish2opus) {
1538: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1540: if (!a->sampler) a->sampler = new PetscMatrixSampler();
1541: a->sampler->SetSamplingMat(B);
1542: if (bs > 0) a->bs = bs;
1543: if (tol > 0.) a->rtol = tol;
1544: delete a->kernel;
1545: }
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: /*@C
1550: MatCreateH2OpusFromKernel - Creates a `MATH2OPUS` from a user-supplied kernel.
1552: Input Parameters:
1553: + comm - MPI communicator
1554: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1555: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1556: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1557: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1558: . spacedim - dimension of the space coordinates
1559: . coords - coordinates of the points
1560: . cdist - whether or not coordinates are distributed
1561: . kernel - computational kernel (or `NULL`)
1562: . kernelctx - kernel context
1563: . eta - admissibility condition tolerance
1564: . leafsize - leaf size in cluster tree
1565: - basisord - approximation order for Chebychev interpolation of low-rank blocks
1567: Output Parameter:
1568: . nA - matrix
1570: Options Database Keys:
1571: + -mat_h2opus_leafsize <`PetscInt`> - Leaf size of cluster tree
1572: . -mat_h2opus_eta <`PetscReal`> - Admissibility condition tolerance
1573: . -mat_h2opus_order <`PetscInt`> - Chebychev approximation order
1574: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be used when estimating norms
1576: Level: intermediate
1578: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`
1579: @*/
1580: 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)
1581: {
1582: Mat A;
1583: Mat_H2OPUS *h2opus;
1584: PetscBool iscpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
1586: PetscFunctionBegin;
1587: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1588: PetscCall(MatCreate(comm, &A));
1589: PetscCall(MatSetSizes(A, m, n, M, N));
1590: PetscCheck(M == N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1591: PetscCall(MatSetType(A, MATH2OPUS));
1592: PetscCall(MatBindToCPU(A, iscpu));
1593: PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, kernel, kernelctx));
1595: h2opus = (Mat_H2OPUS *)A->data;
1596: if (eta > 0.) h2opus->eta = eta;
1597: if (leafsize > 0) h2opus->leafsize = leafsize;
1598: if (basisord > 0) h2opus->basisord = basisord;
1600: *nA = A;
1601: PetscFunctionReturn(PETSC_SUCCESS);
1602: }
1604: /*@
1605: MatCreateH2OpusFromMat - Creates a `MATH2OPUS` sampling from a user-supplied operator.
1607: Input Parameters:
1608: + B - the matrix to be sampled
1609: . spacedim - dimension of the space coordinates
1610: . coords - coordinates of the points
1611: . cdist - whether or not coordinates are distributed
1612: . eta - admissibility condition tolerance
1613: . leafsize - leaf size in cluster tree
1614: . maxrank - maximum rank allowed
1615: . bs - maximum number of samples to be taken concurrently
1616: - rtol - relative tolerance for construction
1618: Output Parameter:
1619: . nA - matrix
1621: Options Database Keys:
1622: + -mat_h2opus_leafsize <`PetscInt`> - Leaf size of cluster tree
1623: . -mat_h2opus_eta <`PetscReal`> - Admissibility condition tolerance
1624: . -mat_h2opus_maxrank <`PetscInt`> - Maximum rank when constructed from matvecs
1625: . -mat_h2opus_samples <`PetscInt`> - Maximum number of samples to be taken concurrently when constructing from matvecs
1626: . -mat_h2opus_rtol <`PetscReal`> - Relative tolerance for construction from sampling
1627: . -mat_h2opus_check <`PetscBool`> - Check error when constructing from sampling during MatAssemblyEnd()
1628: . -mat_h2opus_hara_verbose <`PetscBool`> - Verbose output from hara construction
1629: - -mat_h2opus_normsamples <`PetscInt`> - Maximum number of samples to be when estimating norms
1631: Level: intermediate
1633: Note:
1634: Not available in parallel
1636: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()`
1637: @*/
1638: PetscErrorCode MatCreateH2OpusFromMat(Mat B, PetscInt spacedim, const PetscReal coords[], PetscBool cdist, PetscReal eta, PetscInt leafsize, PetscInt maxrank, PetscInt bs, PetscReal rtol, Mat *nA)
1639: {
1640: Mat A;
1641: Mat_H2OPUS *h2opus;
1642: MPI_Comm comm;
1643: PetscBool boundtocpu = PETSC_TRUE;
1645: PetscFunctionBegin;
1654: PetscAssertPointer(nA, 10);
1655: PetscCall(PetscObjectGetComm((PetscObject)B, &comm));
1656: PetscCheck(B->rmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Different row and column local sizes are not supported");
1657: PetscCheck(B->rmap->N == B->cmap->N, comm, PETSC_ERR_SUP, "Rectangular matrices are not supported");
1658: PetscCall(MatCreate(comm, &A));
1659: PetscCall(MatSetSizes(A, B->rmap->n, B->cmap->n, B->rmap->N, B->cmap->N));
1660: #if defined(PETSC_H2OPUS_USE_GPU)
1661: {
1662: VecType vtype;
1663: PetscBool isstd, iscuda, iskok;
1665: PetscCall(MatGetVecType(B, &vtype));
1666: PetscCall(PetscStrcmpAny(vtype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
1667: PetscCall(PetscStrcmpAny(vtype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
1668: PetscCall(PetscStrcmpAny(vtype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
1669: PetscCheck(isstd || iscuda || iskok, comm, PETSC_ERR_SUP, "Not for type %s", vtype);
1670: if (iscuda && !B->boundtocpu) boundtocpu = PETSC_FALSE;
1671: if (iskok && PetscDefined(HAVE_MACRO_KOKKOS_ENABLE_CUDA)) boundtocpu = PETSC_FALSE;
1672: }
1673: #endif
1674: PetscCall(MatSetType(A, MATH2OPUS));
1675: PetscCall(MatBindToCPU(A, boundtocpu));
1676: if (spacedim) PetscCall(MatH2OpusSetCoords_H2OPUS(A, spacedim, coords, cdist, NULL, NULL));
1677: PetscCall(MatPropagateSymmetryOptions(B, A));
1678: /* PetscCheck(A->symmetric,comm,PETSC_ERR_SUP,"Unsymmetric sampling does not work"); */
1680: h2opus = (Mat_H2OPUS *)A->data;
1681: h2opus->sampler = new PetscMatrixSampler(B);
1682: if (eta > 0.) h2opus->eta = eta;
1683: if (leafsize > 0) h2opus->leafsize = leafsize;
1684: if (maxrank > 0) h2opus->max_rank = maxrank;
1685: if (bs > 0) h2opus->bs = bs;
1686: if (rtol > 0.) h2opus->rtol = rtol;
1687: *nA = A;
1688: A->preallocated = PETSC_TRUE;
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: MatH2OpusGetIndexMap - Access reordering index set.
1695: Input Parameter:
1696: . A - the matrix
1698: Output Parameter:
1699: . indexmap - the index set for the reordering
1701: Level: intermediate
1703: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1704: @*/
1705: PetscErrorCode MatH2OpusGetIndexMap(Mat A, IS *indexmap)
1706: {
1707: PetscBool ish2opus;
1708: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1710: PetscFunctionBegin;
1713: PetscAssertPointer(indexmap, 2);
1714: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1715: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1716: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1717: *indexmap = a->h2opus_indexmap;
1718: PetscFunctionReturn(PETSC_SUCCESS);
1719: }
1721: /*@
1722: MatH2OpusMapVec - Maps a vector between PETSc and H2Opus ordering
1724: Input Parameters:
1725: + A - the matrix
1726: . nativetopetsc - if true, maps from H2Opus ordering to PETSc ordering. If false, applies the reverse map
1727: - in - the vector to be mapped
1729: Output Parameter:
1730: . out - the newly created mapped vector
1732: Level: intermediate
1734: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`
1735: @*/
1736: PetscErrorCode MatH2OpusMapVec(Mat A, PetscBool nativetopetsc, Vec in, Vec *out)
1737: {
1738: PetscBool ish2opus;
1739: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1740: PetscScalar *xin, *xout;
1741: PetscBool nm;
1743: PetscFunctionBegin;
1748: PetscAssertPointer(out, 4);
1749: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1750: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &ish2opus));
1751: PetscCheck(ish2opus, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
1752: nm = a->nativemult;
1753: PetscCall(MatH2OpusSetNativeMult(A, (PetscBool)!nativetopetsc));
1754: PetscCall(MatCreateVecs(A, out, NULL));
1755: PetscCall(MatH2OpusSetNativeMult(A, nm));
1756: if (!a->sf) { /* same ordering */
1757: PetscCall(VecCopy(in, *out));
1758: PetscFunctionReturn(PETSC_SUCCESS);
1759: }
1760: PetscCall(VecGetArrayRead(in, (const PetscScalar **)&xin));
1761: PetscCall(VecGetArrayWrite(*out, &xout));
1762: if (nativetopetsc) {
1763: PetscCall(PetscSFReduceBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1764: PetscCall(PetscSFReduceEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1765: } else {
1766: PetscCall(PetscSFBcastBegin(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1767: PetscCall(PetscSFBcastEnd(a->sf, MPIU_SCALAR, xin, xout, MPI_REPLACE));
1768: }
1769: PetscCall(VecRestoreArrayRead(in, (const PetscScalar **)&xin));
1770: PetscCall(VecRestoreArrayWrite(*out, &xout));
1771: PetscFunctionReturn(PETSC_SUCCESS);
1772: }
1774: /*@
1775: MatH2OpusLowRankUpdate - Perform a low-rank update of the form $ A = A + s * U * V^T $
1777: Input Parameters:
1778: + A - the hierarchical `MATH2OPUS` matrix
1779: . s - the scaling factor
1780: . U - the dense low-rank update matrix
1781: - V - (optional) the dense low-rank update matrix (if `NULL`, then `V` = `U` is assumed)
1783: Note:
1784: The `U` and `V` matrices must be in `MATDENSE` dense format
1786: Level: intermediate
1788: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MATH2OPUS`, `MatCreateH2OpusFromMat()`, `MatCreateH2OpusFromKernel()`, `MatH2OpusCompress()`, `MatH2OpusOrthogonalize()`, `MATDENSE`
1789: @*/
1790: PetscErrorCode MatH2OpusLowRankUpdate(Mat A, Mat U, Mat V, PetscScalar s)
1791: {
1792: PetscBool flg;
1794: PetscFunctionBegin;
1797: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1799: PetscCheckSameComm(A, 1, U, 2);
1800: if (V) {
1802: PetscCheckSameComm(A, 1, V, 3);
1803: }
1806: if (!V) V = U;
1807: 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);
1808: if (!U->cmap->N) PetscFunctionReturn(PETSC_SUCCESS);
1809: PetscCall(PetscLayoutCompare(U->rmap, A->rmap, &flg));
1810: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A and U must have the same row layout");
1811: PetscCall(PetscLayoutCompare(V->rmap, A->cmap, &flg));
1812: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "A column layout must match V row column layout");
1813: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATH2OPUS, &flg));
1814: if (flg) {
1815: Mat_H2OPUS *a = (Mat_H2OPUS *)A->data;
1816: const PetscScalar *u, *v, *uu, *vv;
1817: PetscInt ldu, ldv;
1818: PetscMPIInt size;
1819: #if defined(H2OPUS_USE_MPI)
1820: h2opusHandle_t handle = a->handle->handle;
1821: #else
1822: h2opusHandle_t handle = a->handle;
1823: #endif
1824: PetscBool usesf = (PetscBool)(a->sf && !a->nativemult);
1825: PetscSF usf, vsf;
1827: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1828: PetscCheck(size <= 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not yet implemented in parallel");
1829: PetscCall(PetscLogEventBegin(MAT_H2Opus_LR, A, 0, 0, 0));
1830: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)U, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1831: PetscCheck(flg, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "Not for U of type %s", ((PetscObject)U)->type_name);
1832: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)V, &flg, MATSEQDENSE, MATMPIDENSE, ""));
1833: PetscCheck(flg, PetscObjectComm((PetscObject)V), PETSC_ERR_SUP, "Not for V of type %s", ((PetscObject)V)->type_name);
1834: PetscCall(MatDenseGetLDA(U, &ldu));
1835: PetscCall(MatDenseGetLDA(V, &ldv));
1836: PetscCall(MatBoundToCPU(A, &flg));
1837: if (usesf) {
1838: PetscInt n;
1840: PetscCall(MatDenseGetH2OpusStridedSF(U, a->sf, &usf));
1841: PetscCall(MatDenseGetH2OpusStridedSF(V, a->sf, &vsf));
1842: PetscCall(MatH2OpusResizeBuffers_Private(A, U->cmap->N, V->cmap->N));
1843: PetscCall(PetscSFGetGraph(a->sf, NULL, &n, NULL, NULL));
1844: ldu = n;
1845: ldv = n;
1846: }
1847: if (flg) {
1848: PetscCheck(a->hmatrix, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing CPU matrix");
1849: PetscCall(MatDenseGetArrayRead(U, &u));
1850: PetscCall(MatDenseGetArrayRead(V, &v));
1851: if (usesf) {
1852: vv = MatH2OpusGetThrustPointer(*a->yy);
1853: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1854: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1855: if (U != V) {
1856: uu = MatH2OpusGetThrustPointer(*a->xx);
1857: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1858: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1859: } else uu = vv;
1860: } else {
1861: uu = u;
1862: vv = v;
1863: }
1864: hlru_global(*a->hmatrix, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1865: PetscCall(MatDenseRestoreArrayRead(U, &u));
1866: PetscCall(MatDenseRestoreArrayRead(V, &v));
1867: } else {
1868: #if defined(PETSC_H2OPUS_USE_GPU)
1869: PetscBool flgU, flgV;
1871: PetscCheck(a->hmatrix_gpu, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing GPU matrix");
1872: PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &flgU, MATSEQDENSE, MATMPIDENSE, ""));
1873: if (flgU) PetscCall(MatConvert(U, MATDENSECUDA, MAT_INPLACE_MATRIX, &U));
1874: PetscCall(PetscObjectTypeCompareAny((PetscObject)V, &flgV, MATSEQDENSE, MATMPIDENSE, ""));
1875: if (flgV) PetscCall(MatConvert(V, MATDENSECUDA, MAT_INPLACE_MATRIX, &V));
1876: PetscCall(MatDenseCUDAGetArrayRead(U, &u));
1877: PetscCall(MatDenseCUDAGetArrayRead(V, &v));
1878: if (usesf) {
1879: vv = MatH2OpusGetThrustPointer(*a->yy_gpu);
1880: PetscCall(PetscSFBcastBegin(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1881: PetscCall(PetscSFBcastEnd(vsf, MPIU_SCALAR, v, (PetscScalar *)vv, MPI_REPLACE));
1882: if (U != V) {
1883: uu = MatH2OpusGetThrustPointer(*a->xx_gpu);
1884: PetscCall(PetscSFBcastBegin(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1885: PetscCall(PetscSFBcastEnd(usf, MPIU_SCALAR, u, (PetscScalar *)uu, MPI_REPLACE));
1886: } else uu = vv;
1887: } else {
1888: uu = u;
1889: vv = v;
1890: }
1891: #else
1892: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
1893: #endif
1894: hlru_global(*a->hmatrix_gpu, uu, ldu, vv, ldv, U->cmap->N, s, handle);
1895: #if defined(PETSC_H2OPUS_USE_GPU)
1896: PetscCall(MatDenseCUDARestoreArrayRead(U, &u));
1897: PetscCall(MatDenseCUDARestoreArrayRead(V, &v));
1898: if (flgU) PetscCall(MatConvert(U, MATDENSE, MAT_INPLACE_MATRIX, &U));
1899: if (flgV) PetscCall(MatConvert(V, MATDENSE, MAT_INPLACE_MATRIX, &V));
1900: #endif
1901: }
1902: PetscCall(PetscLogEventEnd(MAT_H2Opus_LR, A, 0, 0, 0));
1903: a->orthogonal = PETSC_FALSE;
1904: }
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1907: #endif