Actual source code: mhypre.c
1: /*
2: Creates hypre ijmatrix from PETSc matrix
3: */
5: #include <petscpkg_version.h>
6: #include <petsc/private/petschypre.h>
7: #include <petscmathypre.h>
8: #include <petsc/private/matimpl.h>
9: #include <petsc/private/deviceimpl.h>
10: #include <../src/mat/impls/hypre/mhypre.h>
11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
12: #include <../src/vec/vec/impls/hypre/vhyp.h>
13: #include <HYPRE.h>
14: #include <HYPRE_utilities.h>
15: #include <_hypre_parcsr_ls.h>
16: #include <_hypre_sstruct_ls.h>
18: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
19: #define hypre_ParCSRMatrixClone(A, B) hypre_ParCSRMatrixCompleteClone(A)
20: #endif
22: #if PETSC_PKG_HYPRE_VERSION_GE(2, 15, 0)
23: #define HYPRE_AssumedPartitionCheck() 1
24: #endif
26: static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *);
27: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix);
28: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat, HYPRE_IJMatrix);
29: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat, HYPRE_IJMatrix);
30: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool);
31: static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins);
33: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
34: {
35: PetscInt i, n_d, n_o;
36: const PetscInt *ia_d, *ia_o;
37: PetscBool done_d = PETSC_FALSE, done_o = PETSC_FALSE;
38: HYPRE_Int *nnz_d = NULL, *nnz_o = NULL;
40: PetscFunctionBegin;
41: if (A_d) { /* determine number of nonzero entries in local diagonal part */
42: PetscCall(MatGetRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, &n_d, &ia_d, NULL, &done_d));
43: if (done_d) {
44: PetscCall(PetscMalloc1(n_d, &nnz_d));
45: for (i = 0; i < n_d; i++) nnz_d[i] = (HYPRE_Int)(ia_d[i + 1] - ia_d[i]);
46: }
47: PetscCall(MatRestoreRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, NULL, &ia_d, NULL, &done_d));
48: }
49: if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
50: PetscCall(MatGetRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
51: if (done_o) {
52: PetscCall(PetscMalloc1(n_o, &nnz_o));
53: for (i = 0; i < n_o; i++) nnz_o[i] = (HYPRE_Int)(ia_o[i + 1] - ia_o[i]);
54: }
55: PetscCall(MatRestoreRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
56: }
57: if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */
58: if (!done_o) { /* only diagonal part */
59: PetscCall(PetscCalloc1(n_d, &nnz_o));
60: }
61: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
62: { /* If we don't do this, the columns of the matrix will be all zeros! */
63: hypre_AuxParCSRMatrix *aux_matrix;
64: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
65: hypre_AuxParCSRMatrixDestroy(aux_matrix);
66: hypre_IJMatrixTranslator(ij) = NULL;
67: PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(ij, nnz_d, nnz_o));
68: /* it seems they partially fixed it in 2.19.0 */
69: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
70: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
71: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
72: #endif
73: }
74: #else
75: PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(ij, nnz_d, nnz_o));
76: #endif
77: PetscCall(PetscFree(nnz_d));
78: PetscCall(PetscFree(nnz_o));
79: }
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
84: {
85: HYPRE_Int rstart, rend, cstart, cend;
87: PetscFunctionBegin;
88: PetscCall(PetscLayoutSetUp(A->rmap));
89: PetscCall(PetscLayoutSetUp(A->cmap));
90: rstart = (HYPRE_Int)A->rmap->rstart;
91: rend = (HYPRE_Int)A->rmap->rend;
92: cstart = (HYPRE_Int)A->cmap->rstart;
93: cend = (HYPRE_Int)A->cmap->rend;
94: PetscHYPREInitialize();
95: if (hA->ij) {
96: if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
97: PetscCallHYPRE(HYPRE_IJMatrixDestroy(hA->ij));
98: }
99: PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij));
100: PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
101: {
102: PetscBool same;
103: Mat A_d, A_o;
104: const PetscInt *colmap;
105: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &same));
106: if (same) {
107: PetscCall(MatMPIAIJGetSeqAIJ(A, &A_d, &A_o, &colmap));
108: PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
111: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &same));
112: if (same) {
113: PetscCall(MatMPIBAIJGetSeqBAIJ(A, &A_d, &A_o, &colmap));
114: PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
117: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &same));
118: if (same) {
119: PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
122: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQBAIJ, &same));
123: if (same) {
124: PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ(Mat A, HYPRE_IJMatrix ij)
132: {
133: PetscBool flg;
135: PetscFunctionBegin;
136: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
137: PetscCallHYPRE(HYPRE_IJMatrixInitialize(ij));
138: #else
139: PetscCallHYPRE(HYPRE_IJMatrixInitialize_v2(ij, HYPRE_MEMORY_HOST));
140: #endif
141: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg));
142: if (flg) {
143: PetscCall(MatHYPRE_IJMatrixCopyIJ_MPIAIJ(A, ij));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
146: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg));
147: if (flg) {
148: PetscCall(MatHYPRE_IJMatrixCopyIJ_SeqAIJ(A, ij));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
151: PetscCheck(PETSC_FALSE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for matrix type %s", ((PetscObject)A)->type_name);
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
156: {
157: Mat_SeqAIJ *pdiag = (Mat_SeqAIJ *)A->data;
158: HYPRE_Int type;
159: hypre_ParCSRMatrix *par_matrix;
160: hypre_AuxParCSRMatrix *aux_matrix;
161: hypre_CSRMatrix *hdiag;
162: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
164: PetscFunctionBegin;
165: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(ij, &type));
166: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
167: PetscCallHYPRE(HYPRE_IJMatrixGetObject(ij, (void **)&par_matrix));
168: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
169: /*
170: this is the Hack part where we monkey directly with the hypre datastructures
171: */
172: if (sameint) {
173: PetscCall(PetscArraycpy(hdiag->i, pdiag->i, A->rmap->n + 1));
174: PetscCall(PetscArraycpy(hdiag->j, pdiag->j, pdiag->nz));
175: } else {
176: PetscInt i;
178: for (i = 0; i < A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
179: for (i = 0; i < pdiag->nz; i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
180: }
182: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
183: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
188: {
189: Mat_MPIAIJ *pA = (Mat_MPIAIJ *)A->data;
190: Mat_SeqAIJ *pdiag, *poffd;
191: PetscInt i, *garray = pA->garray, *jj, cstart, *pjj;
192: HYPRE_Int *hjj, type;
193: hypre_ParCSRMatrix *par_matrix;
194: hypre_AuxParCSRMatrix *aux_matrix;
195: hypre_CSRMatrix *hdiag, *hoffd;
196: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
198: PetscFunctionBegin;
199: pdiag = (Mat_SeqAIJ *)pA->A->data;
200: poffd = (Mat_SeqAIJ *)pA->B->data;
201: /* cstart is only valid for square MPIAIJ laid out in the usual way */
202: PetscCall(MatGetOwnershipRange(A, &cstart, NULL));
204: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(ij, &type));
205: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
206: PetscCallHYPRE(HYPRE_IJMatrixGetObject(ij, (void **)&par_matrix));
207: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
208: hoffd = hypre_ParCSRMatrixOffd(par_matrix);
210: if (sameint) {
211: PetscCall(PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1));
212: } else {
213: for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
214: }
216: hjj = hdiag->j;
217: pjj = pdiag->j;
218: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
219: for (i = 0; i < pdiag->nz; i++) hjj[i] = (HYPRE_Int)pjj[i];
220: #else
221: for (i = 0; i < pdiag->nz; i++) hjj[i] = cstart + pjj[i];
222: #endif
223: if (sameint) {
224: PetscCall(PetscArraycpy(hoffd->i, poffd->i, pA->A->rmap->n + 1));
225: } else {
226: for (i = 0; i < pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)poffd->i[i];
227: }
229: jj = (PetscInt *)hoffd->j;
230: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
231: PetscCallHYPRE(hypre_CSRMatrixBigInitialize(hoffd));
232: jj = (PetscInt *)hoffd->big_j;
233: #endif
234: pjj = poffd->j;
235: for (i = 0; i < poffd->nz; i++) jj[i] = garray[pjj[i]];
237: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
238: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat *B)
243: {
244: Mat_HYPRE *mhA = (Mat_HYPRE *)A->data;
245: Mat lA;
246: ISLocalToGlobalMapping rl2g, cl2g;
247: IS is;
248: hypre_ParCSRMatrix *hA;
249: hypre_CSRMatrix *hdiag, *hoffd;
250: MPI_Comm comm;
251: HYPRE_Complex *hdd, *hod, *aa;
252: PetscScalar *data;
253: HYPRE_BigInt *col_map_offd;
254: HYPRE_Int *hdi, *hdj, *hoi, *hoj;
255: PetscInt *ii, *jj, *iptr, *jptr;
256: PetscInt cum, dr, dc, oc, str, stc, nnz, i, jd, jo, M, N;
257: HYPRE_Int type;
258: MatType lmattype = NULL;
259: PetscBool freeparcsr = PETSC_FALSE;
261: PetscFunctionBegin;
262: comm = PetscObjectComm((PetscObject)A);
263: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(mhA->ij, &type));
264: PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
265: PetscCallHYPRE(HYPRE_IJMatrixGetObject(mhA->ij, (void **)&hA));
266: #if defined(PETSC_HAVE_HYPRE_DEVICE)
267: if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(mhA->ij)) {
268: /* Support by copying back on the host and copy to GPU
269: Kind of inefficient, but this is the best we can do now */
270: #if defined(HYPRE_USING_HIP)
271: lmattype = MATSEQAIJHIPSPARSE;
272: #elif defined(HYPRE_USING_CUDA)
273: lmattype = MATSEQAIJCUSPARSE;
274: #endif
275: hA = hypre_ParCSRMatrixClone_v2(hA, 1, HYPRE_MEMORY_HOST);
276: freeparcsr = PETSC_TRUE;
277: }
278: #endif
279: M = hypre_ParCSRMatrixGlobalNumRows(hA);
280: N = hypre_ParCSRMatrixGlobalNumCols(hA);
281: str = hypre_ParCSRMatrixFirstRowIndex(hA);
282: stc = hypre_ParCSRMatrixFirstColDiag(hA);
283: hdiag = hypre_ParCSRMatrixDiag(hA);
284: hoffd = hypre_ParCSRMatrixOffd(hA);
285: dr = hypre_CSRMatrixNumRows(hdiag);
286: dc = hypre_CSRMatrixNumCols(hdiag);
287: nnz = hypre_CSRMatrixNumNonzeros(hdiag);
288: hdi = hypre_CSRMatrixI(hdiag);
289: hdj = hypre_CSRMatrixJ(hdiag);
290: hdd = hypre_CSRMatrixData(hdiag);
291: oc = hypre_CSRMatrixNumCols(hoffd);
292: nnz += hypre_CSRMatrixNumNonzeros(hoffd);
293: hoi = hypre_CSRMatrixI(hoffd);
294: hoj = hypre_CSRMatrixJ(hoffd);
295: hod = hypre_CSRMatrixData(hoffd);
296: if (reuse != MAT_REUSE_MATRIX) {
297: PetscInt *aux;
299: /* generate l2g maps for rows and cols */
300: PetscCall(ISCreateStride(comm, dr, str, 1, &is));
301: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
302: PetscCall(ISDestroy(&is));
303: col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
304: PetscCall(PetscMalloc1(dc + oc, &aux));
305: for (i = 0; i < dc; i++) aux[i] = i + stc;
306: for (i = 0; i < oc; i++) aux[i + dc] = col_map_offd[i];
307: PetscCall(ISCreateGeneral(comm, dc + oc, aux, PETSC_OWN_POINTER, &is));
308: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
309: PetscCall(ISDestroy(&is));
310: /* create MATIS object */
311: PetscCall(MatCreate(comm, B));
312: PetscCall(MatSetSizes(*B, dr, dc, M, N));
313: PetscCall(MatSetType(*B, MATIS));
314: PetscCall(MatSetLocalToGlobalMapping(*B, rl2g, cl2g));
315: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
316: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
318: /* allocate CSR for local matrix */
319: PetscCall(PetscMalloc1(dr + 1, &iptr));
320: PetscCall(PetscMalloc1(nnz, &jptr));
321: PetscCall(PetscMalloc1(nnz, &data));
322: } else {
323: PetscInt nr;
324: PetscBool done;
325: PetscCall(MatISGetLocalMat(*B, &lA));
326: PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&iptr, (const PetscInt **)&jptr, &done));
327: PetscCheck(nr == dr, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of rows in local mat! %" PetscInt_FMT " != %" PetscInt_FMT, nr, dr);
328: PetscCheck(iptr[nr] >= nnz, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of nonzeros in local mat! reuse %" PetscInt_FMT " requested %" PetscInt_FMT, iptr[nr], nnz);
329: PetscCall(MatSeqAIJGetArrayWrite(lA, &data));
330: }
331: /* merge local matrices */
332: ii = iptr;
333: jj = jptr;
334: aa = (HYPRE_Complex *)data; /* this cast fixes the clang error when doing the assignments below: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
335: *ii = *(hdi++) + *(hoi++);
336: for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
337: PetscScalar *aold = (PetscScalar *)aa;
338: PetscInt *jold = jj, nc = jd + jo;
339: for (; jd < *hdi; jd++) {
340: *jj++ = *hdj++;
341: *aa++ = *hdd++;
342: }
343: for (; jo < *hoi; jo++) {
344: *jj++ = *hoj++ + dc;
345: *aa++ = *hod++;
346: }
347: *(++ii) = *(hdi++) + *(hoi++);
348: PetscCall(PetscSortIntWithScalarArray(jd + jo - nc, jold, aold));
349: }
350: for (; cum < dr; cum++) *(++ii) = nnz;
351: if (reuse != MAT_REUSE_MATRIX) {
352: Mat_SeqAIJ *a;
354: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, dc + oc, iptr, jptr, data, &lA));
355: /* hack SeqAIJ */
356: a = (Mat_SeqAIJ *)lA->data;
357: a->free_a = PETSC_TRUE;
358: a->free_ij = PETSC_TRUE;
359: if (lmattype) PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));
360: PetscCall(MatISSetLocalMat(*B, lA));
361: PetscCall(MatDestroy(&lA));
362: } else {
363: PetscCall(MatSeqAIJRestoreArrayWrite(lA, &data));
364: }
365: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
366: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
367: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, B));
368: if (freeparcsr) PetscCallHYPRE(hypre_ParCSRMatrixDestroy(hA));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode MatHYPRE_DestroyCOOMat(Mat mat)
373: {
374: Mat_HYPRE *hA = (Mat_HYPRE *)mat->data;
376: PetscFunctionBegin;
377: if (hA->cooMat) { /* If cooMat is present we need to destroy the column indices */
378: PetscCall(MatDestroy(&hA->cooMat));
379: if (hA->cooMatAttached) {
380: hypre_CSRMatrix *csr;
381: hypre_ParCSRMatrix *parcsr;
382: HYPRE_MemoryLocation mem;
384: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
385: csr = hypre_ParCSRMatrixDiag(parcsr);
386: if (csr) {
387: mem = hypre_CSRMatrixMemoryLocation(csr);
388: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
389: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
390: }
391: csr = hypre_ParCSRMatrixOffd(parcsr);
392: if (csr) {
393: mem = hypre_CSRMatrixMemoryLocation(csr);
394: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
395: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
396: }
397: }
398: }
399: hA->cooMatAttached = PETSC_FALSE;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode MatHYPRE_CreateCOOMat(Mat mat)
404: {
405: MPI_Comm comm;
406: PetscMPIInt size;
407: PetscLayout rmap, cmap;
408: Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
409: MatType matType = MATAIJ; /* default type of cooMat */
411: PetscFunctionBegin;
412: /* Build an agent matrix cooMat with AIJ format
413: It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
414: */
415: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
416: PetscCallMPI(MPI_Comm_size(comm, &size));
417: PetscCall(PetscLayoutSetUp(mat->rmap));
418: PetscCall(PetscLayoutSetUp(mat->cmap));
419: PetscCall(MatGetLayouts(mat, &rmap, &cmap));
421: #if defined(PETSC_HAVE_HYPRE_DEVICE)
422: if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
423: #if defined(HYPRE_USING_HIP)
424: matType = MATAIJHIPSPARSE;
425: #elif defined(HYPRE_USING_CUDA)
426: matType = MATAIJCUSPARSE;
427: #elif defined(HYPRE_USING_SYCL) && defined(PETSC_HAVE_KOKKOS_KERNELS)
428: matType = MATAIJKOKKOS;
429: #else
430: SETERRQ(comm, PETSC_ERR_SUP, "No HYPRE device available. Suggest re-installing with Kokkos Kernels");
431: #endif
432: }
433: #endif
435: /* Do COO preallocation through cooMat */
436: PetscCall(MatHYPRE_DestroyCOOMat(mat));
437: PetscCall(MatCreate(comm, &hmat->cooMat));
438: PetscCall(MatSetType(hmat->cooMat, matType));
439: PetscCall(MatSetLayouts(hmat->cooMat, rmap, cmap));
441: /* allocate local matrices if needed */
442: PetscCall(MatMPIAIJSetPreallocation(hmat->cooMat, 0, NULL, 0, NULL));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /* Attach cooMat data array to hypre matrix.
447: When AIJCUPMSPARSE will support raw device pointers and not THRUSTARRAY
448: we should swap the arrays: i.e., attach hypre matrix array to cooMat
449: This is because hypre should be in charge of handling the memory,
450: cooMat is only a way to reuse PETSc COO code.
451: attaching the memory will then be done at MatSetValuesCOO time and it will dynamically
452: support hypre matrix migrating to host.
453: */
454: static PetscErrorCode MatHYPRE_AttachCOOMat(Mat mat)
455: {
456: Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
457: hypre_CSRMatrix *diag, *offd;
458: hypre_ParCSRMatrix *parCSR;
459: HYPRE_MemoryLocation hmem = HYPRE_MEMORY_HOST;
460: PetscMemType pmem;
461: Mat A, B;
462: PetscScalar *a;
463: PetscMPIInt size;
464: MPI_Comm comm;
466: PetscFunctionBegin;
467: PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
468: if (hmat->cooMatAttached) PetscFunctionReturn(PETSC_SUCCESS);
469: PetscCheck(hmat->cooMat->preallocated, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix is not preallocated");
470: PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
471: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
472: PetscCallMPI(MPI_Comm_size(comm, &size));
474: /* Alias cooMat's data array to IJMatrix's */
475: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hmat->ij, (void **)&parCSR));
476: diag = hypre_ParCSRMatrixDiag(parCSR);
477: offd = hypre_ParCSRMatrixOffd(parCSR);
479: A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A;
480: B = (size == 1) ? NULL : ((Mat_MPIAIJ *)hmat->cooMat->data)->B;
482: PetscCall(PetscObjectSetName((PetscObject)A, "_internal_COO_mat_for_hypre"));
483: hmem = hypre_CSRMatrixMemoryLocation(diag);
484: PetscCall(MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &a, &pmem));
485: PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
486: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hmem));
487: hypre_CSRMatrixData(diag) = (HYPRE_Complex *)a;
488: hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */
490: if (B) {
491: hmem = hypre_CSRMatrixMemoryLocation(offd);
492: PetscCall(MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &a, &pmem));
493: PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
494: PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hmem));
495: hypre_CSRMatrixData(offd) = (HYPRE_Complex *)a;
496: hypre_CSRMatrixOwnsData(offd) = 0;
497: }
498: hmat->cooMatAttached = PETSC_TRUE;
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: // Build COO's coordinate list i[], j[] based on CSR's i[], j[] arrays and the number of local rows 'n'
503: static PetscErrorCode CSRtoCOO_Private(PetscInt n, const PetscInt ii[], const PetscInt jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
504: {
505: PetscInt *cooi, *cooj;
507: PetscFunctionBegin;
508: *ncoo = ii[n];
509: PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
510: for (PetscInt i = 0; i < n; i++) {
511: for (PetscInt j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
512: }
513: PetscCall(PetscArraycpy(cooj, jj, *ncoo));
514: *coo_i = cooi;
515: *coo_j = cooj;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: // Similar to CSRtoCOO_Private, but the CSR's i[], j[] are of type HYPRE_Int
520: static PetscErrorCode CSRtoCOO_HYPRE_Int_Private(PetscInt n, const HYPRE_Int ii[], const HYPRE_Int jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
521: {
522: PetscInt *cooi, *cooj;
524: PetscFunctionBegin;
525: *ncoo = ii[n];
526: PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
527: for (PetscInt i = 0; i < n; i++) {
528: for (HYPRE_Int j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
529: }
530: for (PetscCount i = 0; i < *ncoo; i++) cooj[i] = jj[i];
531: *coo_i = cooi;
532: *coo_j = cooj;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: // Build a COO data structure for the seqaij matrix, as if the nonzeros are laid out in the same order as in the CSR
537: static PetscErrorCode MatSeqAIJGetCOO_Private(Mat A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
538: {
539: PetscInt n;
540: const PetscInt *ii, *jj;
541: PetscBool done;
543: PetscFunctionBegin;
544: PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
545: PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatGetRowIJ");
546: PetscCall(CSRtoCOO_Private(n, ii, jj, ncoo, coo_i, coo_j));
547: PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
548: PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatRestoreRowIJ");
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: // Build a COO data structure for the hypreCSRMatrix, as if the nonzeros are laid out in the same order as in the hypreCSRMatrix
553: static PetscErrorCode hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix *A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
554: {
555: PetscInt n = hypre_CSRMatrixNumRows(A);
556: HYPRE_Int *ii, *jj;
557: HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
559: PetscFunctionBegin;
560: #if defined(PETSC_HAVE_HYPRE_DEVICE)
561: mem = hypre_CSRMatrixMemoryLocation(A);
562: if (mem != HYPRE_MEMORY_HOST) {
563: PetscCount nnz = hypre_CSRMatrixNumNonzeros(A);
564: PetscCall(PetscMalloc2(n + 1, &ii, nnz, &jj));
565: hypre_TMemcpy(ii, hypre_CSRMatrixI(A), HYPRE_Int, n + 1, HYPRE_MEMORY_HOST, mem);
566: hypre_TMemcpy(jj, hypre_CSRMatrixJ(A), HYPRE_Int, nnz, HYPRE_MEMORY_HOST, mem);
567: } else {
568: #else
569: {
570: #endif
571: ii = hypre_CSRMatrixI(A);
572: jj = hypre_CSRMatrixJ(A);
573: }
574: PetscCall(CSRtoCOO_HYPRE_Int_Private(n, ii, jj, ncoo, coo_i, coo_j));
575: if (mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree2(ii, jj));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode MatSetValuesCOOFromCSRMatrix_Private(Mat A, hypre_CSRMatrix *H)
580: {
581: PetscBool iscpu = PETSC_TRUE;
582: PetscScalar *a;
583: HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
585: PetscFunctionBegin;
586: #if defined(PETSC_HAVE_HYPRE_DEVICE)
587: mem = hypre_CSRMatrixMemoryLocation(H);
588: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &iscpu));
589: #endif
590: if (iscpu && mem != HYPRE_MEMORY_HOST) {
591: PetscCount nnz = hypre_CSRMatrixNumNonzeros(H);
592: PetscCall(PetscMalloc1(nnz, &a));
593: hypre_TMemcpy(a, hypre_CSRMatrixData(H), PetscScalar, nnz, HYPRE_MEMORY_HOST, mem);
594: } else {
595: a = (PetscScalar *)hypre_CSRMatrixData(H);
596: }
597: PetscCall(MatSetValuesCOO(A, a, INSERT_VALUES));
598: if (iscpu && mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree(a));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
603: {
604: MPI_Comm comm = PetscObjectComm((PetscObject)A);
605: Mat M = NULL, dH = NULL, oH = NULL, dA = NULL, oA = NULL;
606: PetscBool ismpiaij, issbaij, isbaij, boundtocpu = PETSC_TRUE;
607: Mat_HYPRE *hA;
608: PetscMemType memtype = PETSC_MEMTYPE_HOST;
610: PetscFunctionBegin;
611: if (PetscDefined(HAVE_HYPRE_DEVICE)) {
612: PetscCall(MatGetCurrentMemType(A, &memtype));
613: PetscHYPREInitialize();
614: boundtocpu = PetscMemTypeHost(memtype) ? PETSC_TRUE : PETSC_FALSE;
615: PetscCallHYPRE(HYPRE_SetMemoryLocation(boundtocpu ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE));
616: }
618: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &issbaij, MATSEQSBAIJ, MATMPIBAIJ, ""));
619: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &isbaij, MATSEQBAIJ, MATMPIBAIJ, ""));
620: if (isbaij || issbaij) { /* handle BAIJ and SBAIJ */
621: PetscBool ismpi;
622: MatType newtype;
624: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &ismpi, MATMPISBAIJ, MATMPIBAIJ, ""));
625: newtype = ismpi ? MATMPIAIJ : MATSEQAIJ;
626: if (reuse == MAT_REUSE_MATRIX) {
627: PetscCall(MatConvert(*B, newtype, MAT_INPLACE_MATRIX, B));
628: PetscCall(MatConvert(A, newtype, MAT_REUSE_MATRIX, B));
629: PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
630: } else if (reuse == MAT_INITIAL_MATRIX) {
631: PetscCall(MatConvert(A, newtype, MAT_INITIAL_MATRIX, B));
632: PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
633: } else {
634: PetscCall(MatConvert(A, newtype, MAT_INPLACE_MATRIX, &A));
635: PetscCall(MatConvert(A, MATHYPRE, MAT_INPLACE_MATRIX, &A));
636: }
637: #if defined(PETSC_HAVE_DEVICE)
638: (*B)->boundtocpu = boundtocpu;
639: #endif
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: dA = A;
644: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
645: if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(A, &dA, &oA, NULL));
647: if (reuse != MAT_REUSE_MATRIX) {
648: PetscCount coo_n;
649: PetscInt *coo_i, *coo_j;
651: PetscCall(MatCreate(comm, &M));
652: PetscCall(MatSetType(M, MATHYPRE));
653: PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
654: PetscCall(MatSetOption(M, MAT_SORTED_FULL, PETSC_TRUE));
655: PetscCall(MatSetOption(M, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
657: hA = (Mat_HYPRE *)M->data;
658: PetscCall(MatHYPRE_CreateFromMat(A, hA));
659: PetscCall(MatHYPRE_IJMatrixCopyIJ(A, hA->ij));
661: PetscCall(MatHYPRE_CreateCOOMat(M));
663: dH = hA->cooMat;
664: PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
665: if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));
667: PetscCall(PetscObjectSetName((PetscObject)dH, "_internal_COO_mat_for_hypre"));
668: PetscCall(MatSeqAIJGetCOO_Private(dA, &coo_n, &coo_i, &coo_j));
669: PetscCall(MatSetPreallocationCOO(dH, coo_n, coo_i, coo_j));
670: PetscCall(PetscFree2(coo_i, coo_j));
671: if (oH) {
672: PetscCall(PetscLayoutDestroy(&oH->cmap));
673: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oH), oA->cmap->n, oA->cmap->n, 1, &oH->cmap));
674: PetscCall(MatSeqAIJGetCOO_Private(oA, &coo_n, &coo_i, &coo_j));
675: PetscCall(MatSetPreallocationCOO(oH, coo_n, coo_i, coo_j));
676: PetscCall(PetscFree2(coo_i, coo_j));
677: }
678: hA->cooMat->assembled = PETSC_TRUE;
680: M->preallocated = PETSC_TRUE;
681: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
682: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
684: PetscCall(MatHYPRE_AttachCOOMat(M));
685: if (reuse == MAT_INITIAL_MATRIX) *B = M;
686: } else M = *B;
688: hA = (Mat_HYPRE *)M->data;
689: PetscCheck(hA->cooMat, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
691: dH = hA->cooMat;
692: PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
693: if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));
695: PetscScalar *a;
696: PetscCall(MatSeqAIJGetCSRAndMemType(dA, NULL, NULL, &a, NULL));
697: PetscCall(MatSetValuesCOO(dH, a, INSERT_VALUES));
698: if (oH) {
699: PetscCall(MatSeqAIJGetCSRAndMemType(oA, NULL, NULL, &a, NULL));
700: PetscCall(MatSetValuesCOO(oH, a, INSERT_VALUES));
701: }
703: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
704: #if defined(PETSC_HAVE_DEVICE)
705: (*B)->boundtocpu = boundtocpu;
706: #endif
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
711: {
712: Mat M, dA = NULL, oA = NULL;
713: hypre_ParCSRMatrix *parcsr;
714: hypre_CSRMatrix *dH, *oH;
715: MPI_Comm comm;
716: PetscBool ismpiaij, isseqaij;
718: PetscFunctionBegin;
719: comm = PetscObjectComm((PetscObject)A);
720: if (reuse == MAT_REUSE_MATRIX) {
721: PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij));
722: PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij));
723: PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Only MATMPIAIJ or MATSEQAIJ base types are supported");
724: }
725: PetscCall(MatHYPREGetParCSR(A, &parcsr));
726: #if defined(PETSC_HAVE_HYPRE_DEVICE)
727: if (HYPRE_MEMORY_DEVICE == hypre_ParCSRMatrixMemoryLocation(parcsr)) {
728: PetscBool isaij;
730: PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
731: if (isaij) {
732: PetscMPIInt size;
734: PetscCallMPI(MPI_Comm_size(comm, &size));
735: #if defined(HYPRE_USING_HIP)
736: mtype = size > 1 ? MATMPIAIJHIPSPARSE : MATSEQAIJHIPSPARSE;
737: #elif defined(HYPRE_USING_CUDA)
738: mtype = size > 1 ? MATMPIAIJCUSPARSE : MATSEQAIJCUSPARSE;
739: #else
740: mtype = size > 1 ? MATMPIAIJ : MATSEQAIJ;
741: #endif
742: }
743: }
744: #endif
745: dH = hypre_ParCSRMatrixDiag(parcsr);
746: oH = hypre_ParCSRMatrixOffd(parcsr);
747: if (reuse != MAT_REUSE_MATRIX) {
748: PetscCount coo_n;
749: PetscInt *coo_i, *coo_j;
751: PetscCall(MatCreate(comm, &M));
752: PetscCall(MatSetType(M, mtype));
753: PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
754: PetscCall(MatMPIAIJSetPreallocation(M, 0, NULL, 0, NULL));
756: dA = M;
757: PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
758: if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
760: PetscCall(hypreCSRMatrixGetCOO_Private(dH, &coo_n, &coo_i, &coo_j));
761: PetscCall(MatSetPreallocationCOO(dA, coo_n, coo_i, coo_j));
762: PetscCall(PetscFree2(coo_i, coo_j));
763: if (ismpiaij) {
764: HYPRE_Int nc = hypre_CSRMatrixNumCols(oH);
766: PetscCall(PetscLayoutDestroy(&oA->cmap));
767: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oA), nc, nc, 1, &oA->cmap));
768: PetscCall(hypreCSRMatrixGetCOO_Private(oH, &coo_n, &coo_i, &coo_j));
769: PetscCall(MatSetPreallocationCOO(oA, coo_n, coo_i, coo_j));
770: PetscCall(PetscFree2(coo_i, coo_j));
772: /* garray */
773: Mat_MPIAIJ *aij = (Mat_MPIAIJ *)M->data;
774: HYPRE_BigInt *harray = hypre_ParCSRMatrixColMapOffd(parcsr);
775: PetscInt *garray;
777: PetscCall(PetscFree(aij->garray));
778: PetscCall(PetscMalloc1(nc, &garray));
779: for (HYPRE_Int i = 0; i < nc; i++) garray[i] = (PetscInt)harray[i];
780: aij->garray = garray;
781: PetscCall(MatSetUpMultiply_MPIAIJ(M));
782: }
783: if (reuse == MAT_INITIAL_MATRIX) *B = M;
784: } else M = *B;
786: dA = M;
787: PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
788: if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
789: PetscCall(MatSetValuesCOOFromCSRMatrix_Private(dA, dH));
790: if (oA) PetscCall(MatSetValuesCOOFromCSRMatrix_Private(oA, oH));
791: M->assembled = PETSC_TRUE;
792: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }
796: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
797: {
798: hypre_ParCSRMatrix *tA;
799: hypre_CSRMatrix *hdiag, *hoffd;
800: Mat_SeqAIJ *diag, *offd;
801: PetscInt *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts;
802: MPI_Comm comm = PetscObjectComm((PetscObject)A);
803: PetscBool ismpiaij, isseqaij;
804: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
805: HYPRE_Int *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL;
806: PetscInt *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL;
807: PetscBool iscuda, iship;
808: #if defined(PETSC_HAVE_DEVICE) && defined(PETSC_HAVE_HYPRE_DEVICE)
809: PetscBool boundtocpu = A->boundtocpu;
810: #else
811: PetscBool boundtocpu = PETSC_TRUE;
812: #endif
814: PetscFunctionBegin;
815: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
816: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
817: PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Unsupported type %s", ((PetscObject)A)->type_name);
818: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
819: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
820: PetscHYPREInitialize();
821: if (ismpiaij) {
822: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;
824: diag = (Mat_SeqAIJ *)a->A->data;
825: offd = (Mat_SeqAIJ *)a->B->data;
826: if (!boundtocpu && (iscuda || iship)) {
827: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
828: if (iscuda) {
829: sameint = PETSC_TRUE;
830: PetscCall(MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
831: PetscCall(MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
832: }
833: #endif
834: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
835: if (iship) {
836: sameint = PETSC_TRUE;
837: PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
838: PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
839: }
840: #endif
841: } else {
842: boundtocpu = PETSC_TRUE;
843: pdi = diag->i;
844: pdj = diag->j;
845: poi = offd->i;
846: poj = offd->j;
847: if (sameint) {
848: hdi = (HYPRE_Int *)pdi;
849: hdj = (HYPRE_Int *)pdj;
850: hoi = (HYPRE_Int *)poi;
851: hoj = (HYPRE_Int *)poj;
852: }
853: }
854: garray = a->garray;
855: noffd = a->B->cmap->N;
856: dnnz = diag->nz;
857: onnz = offd->nz;
858: } else {
859: diag = (Mat_SeqAIJ *)A->data;
860: offd = NULL;
861: if (!boundtocpu && (iscuda || iship)) {
862: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
863: if (iscuda) {
864: sameint = PETSC_TRUE;
865: PetscCall(MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
866: }
867: #endif
868: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
869: if (iship) {
870: sameint = PETSC_TRUE;
871: PetscCall(MatSeqAIJHIPSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
872: }
873: #endif
874: } else {
875: boundtocpu = PETSC_TRUE;
876: pdi = diag->i;
877: pdj = diag->j;
878: if (sameint) {
879: hdi = (HYPRE_Int *)pdi;
880: hdj = (HYPRE_Int *)pdj;
881: }
882: }
883: garray = NULL;
884: noffd = 0;
885: dnnz = diag->nz;
886: onnz = 0;
887: }
889: /* create a temporary ParCSR */
890: if (HYPRE_AssumedPartitionCheck()) {
891: PetscMPIInt myid;
893: PetscCallMPI(MPI_Comm_rank(comm, &myid));
894: row_starts = A->rmap->range + myid;
895: col_starts = A->cmap->range + myid;
896: } else {
897: row_starts = A->rmap->range;
898: col_starts = A->cmap->range;
899: }
900: tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, (HYPRE_Int)noffd, (HYPRE_Int)dnnz, (HYPRE_Int)onnz);
901: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
902: hypre_ParCSRMatrixSetRowStartsOwner(tA, 0);
903: hypre_ParCSRMatrixSetColStartsOwner(tA, 0);
904: #endif
906: /* set diagonal part */
907: hdiag = hypre_ParCSRMatrixDiag(tA);
908: if (!sameint) { /* malloc CSR pointers */
909: PetscCall(PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj));
910: for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)pdi[i];
911: for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)pdj[i];
912: }
913: hypre_CSRMatrixI(hdiag) = hdi;
914: hypre_CSRMatrixJ(hdiag) = hdj;
915: hypre_CSRMatrixData(hdiag) = (HYPRE_Complex *)diag->a;
916: hypre_CSRMatrixNumNonzeros(hdiag) = (HYPRE_Int)diag->nz;
917: hypre_CSRMatrixSetDataOwner(hdiag, 0);
919: /* set off-diagonal part */
920: hoffd = hypre_ParCSRMatrixOffd(tA);
921: if (offd) {
922: if (!sameint) { /* malloc CSR pointers */
923: PetscCall(PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj));
924: for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)poi[i];
925: for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)poj[i];
926: }
927: hypre_CSRMatrixI(hoffd) = hoi;
928: hypre_CSRMatrixJ(hoffd) = hoj;
929: hypre_CSRMatrixData(hoffd) = (HYPRE_Complex *)offd->a;
930: hypre_CSRMatrixNumNonzeros(hoffd) = (HYPRE_Int)offd->nz;
931: hypre_CSRMatrixSetDataOwner(hoffd, 0);
932: }
933: #if defined(PETSC_HAVE_HYPRE_DEVICE)
934: PetscCallHYPRE(hypre_ParCSRMatrixInitialize_v2(tA, !boundtocpu ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST));
935: #else
936: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
937: PetscCallHYPRE(hypre_ParCSRMatrixInitialize(tA));
938: #else
939: PetscCallHYPRE(hypre_ParCSRMatrixInitialize_v2(tA, HYPRE_MEMORY_HOST));
940: #endif
941: #endif
943: /* MatrixSetRownnz comes after MatrixInitialize, so the first uses the right memory location */
944: hypre_CSRMatrixSetRownnz(hdiag);
945: if (offd) hypre_CSRMatrixSetRownnz(hoffd);
947: hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST);
948: hypre_ParCSRMatrixSetNumNonzeros(tA);
949: hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray;
950: if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallHYPRE(hypre_MatvecCommPkgCreate(tA));
951: *hA = tA;
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
956: {
957: hypre_CSRMatrix *hdiag, *hoffd;
958: PetscBool ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
959: PetscBool iscuda, iship;
961: PetscFunctionBegin;
962: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
963: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
964: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
965: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
966: if (iscuda) sameint = PETSC_TRUE;
967: #elif defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
968: if (iship) sameint = PETSC_TRUE;
969: #endif
970: hdiag = hypre_ParCSRMatrixDiag(*hA);
971: hoffd = hypre_ParCSRMatrixOffd(*hA);
972: /* free temporary memory allocated by PETSc
973: set pointers to NULL before destroying tA */
974: if (!sameint) {
975: HYPRE_Int *hi, *hj;
977: hi = hypre_CSRMatrixI(hdiag);
978: hj = hypre_CSRMatrixJ(hdiag);
979: PetscCall(PetscFree2(hi, hj));
980: if (ismpiaij) {
981: hi = hypre_CSRMatrixI(hoffd);
982: hj = hypre_CSRMatrixJ(hoffd);
983: PetscCall(PetscFree2(hi, hj));
984: }
985: }
986: hypre_CSRMatrixI(hdiag) = NULL;
987: hypre_CSRMatrixJ(hdiag) = NULL;
988: hypre_CSRMatrixData(hdiag) = NULL;
989: if (ismpiaij) {
990: hypre_CSRMatrixI(hoffd) = NULL;
991: hypre_CSRMatrixJ(hoffd) = NULL;
992: hypre_CSRMatrixData(hoffd) = NULL;
993: }
994: hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
995: hypre_ParCSRMatrixDestroy(*hA);
996: *hA = NULL;
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /* calls RAP from BoomerAMG:
1001: the resulting ParCSR will not own the column and row starts
1002: It looks like we don't need to have the diagonal entries ordered first */
1003: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
1004: {
1005: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1006: HYPRE_Int P_owns_col_starts, R_owns_row_starts;
1007: #endif
1009: PetscFunctionBegin;
1010: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1011: P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
1012: R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
1013: #endif
1014: /* can be replaced by version test later */
1015: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1016: PetscStackPushExternal("hypre_ParCSRMatrixRAP");
1017: *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP);
1018: PetscStackPop;
1019: #else
1020: PetscCallHYPRE(hypre_BoomerAMGBuildCoarseOperator(hR, hA, hP, hRAP));
1021: PetscCallHYPRE(hypre_ParCSRMatrixSetNumNonzeros(*hRAP));
1022: #endif
1023: /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
1024: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1025: hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0);
1026: hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0);
1027: if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1);
1028: if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1);
1029: #endif
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C)
1034: {
1035: Mat B;
1036: hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL;
1037: Mat_Product *product = C->product;
1039: PetscFunctionBegin;
1040: PetscCall(MatAIJGetParCSR_Private(A, &hA));
1041: PetscCall(MatAIJGetParCSR_Private(P, &hP));
1042: PetscCall(MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP));
1043: PetscCall(MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B));
1045: PetscCall(MatHeaderMerge(C, &B));
1046: C->product = product;
1048: PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1049: PetscCall(MatAIJRestoreParCSR_Private(P, &hP));
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C)
1054: {
1055: PetscFunctionBegin;
1056: PetscCall(MatSetType(C, MATAIJ));
1057: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
1058: C->ops->productnumeric = MatProductNumeric_PtAP;
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C)
1063: {
1064: Mat B;
1065: Mat_HYPRE *hP;
1066: hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL;
1067: HYPRE_Int type;
1068: MPI_Comm comm = PetscObjectComm((PetscObject)A);
1069: PetscBool ishypre;
1071: PetscFunctionBegin;
1072: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1073: PetscCheck(ishypre, comm, PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1074: hP = (Mat_HYPRE *)P->data;
1075: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hP->ij, &type));
1076: PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1077: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hP->ij, (void **)&Pparcsr));
1079: PetscCall(MatAIJGetParCSR_Private(A, &hA));
1080: PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr));
1081: PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1083: /* create temporary matrix and merge to C */
1084: PetscCall(MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B));
1085: PetscCall(MatHeaderMerge(C, &B));
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C)
1090: {
1091: Mat B;
1092: hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL;
1093: Mat_HYPRE *hA, *hP;
1094: PetscBool ishypre;
1095: HYPRE_Int type;
1097: PetscFunctionBegin;
1098: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1099: PetscCheck(ishypre, PetscObjectComm((PetscObject)P), PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1100: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1101: PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1102: hA = (Mat_HYPRE *)A->data;
1103: hP = (Mat_HYPRE *)P->data;
1104: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1105: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1106: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hP->ij, &type));
1107: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1108: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&Aparcsr));
1109: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hP->ij, (void **)&Pparcsr));
1110: PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr));
1111: PetscCall(MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B));
1112: PetscCall(MatHeaderMerge(C, &B));
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: /* calls hypre_ParMatmul
1117: hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
1118: hypre_ParMatrixCreate does not duplicate the communicator
1119: It looks like we don't need to have the diagonal entries ordered first */
1120: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
1121: {
1122: PetscFunctionBegin;
1123: /* can be replaced by version test later */
1124: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1125: PetscStackPushExternal("hypre_ParCSRMatMat");
1126: *hAB = hypre_ParCSRMatMat(hA, hB);
1127: #else
1128: PetscStackPushExternal("hypre_ParMatmul");
1129: *hAB = hypre_ParMatmul(hA, hB);
1130: #endif
1131: PetscStackPop;
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C)
1136: {
1137: Mat D;
1138: hypre_ParCSRMatrix *hA, *hB, *hAB = NULL;
1139: Mat_Product *product = C->product;
1141: PetscFunctionBegin;
1142: PetscCall(MatAIJGetParCSR_Private(A, &hA));
1143: PetscCall(MatAIJGetParCSR_Private(B, &hB));
1144: PetscCall(MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB));
1145: PetscCall(MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D));
1147: PetscCall(MatHeaderMerge(C, &D));
1148: C->product = product;
1150: PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1151: PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C)
1156: {
1157: PetscFunctionBegin;
1158: PetscCall(MatSetType(C, MATAIJ));
1159: C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
1160: C->ops->productnumeric = MatProductNumeric_AB;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C)
1165: {
1166: Mat D;
1167: hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL;
1168: Mat_HYPRE *hA, *hB;
1169: PetscBool ishypre;
1170: HYPRE_Int type;
1171: Mat_Product *product;
1173: PetscFunctionBegin;
1174: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre));
1175: PetscCheck(ishypre, PetscObjectComm((PetscObject)B), PETSC_ERR_USER, "B should be of type %s", MATHYPRE);
1176: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1177: PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1178: hA = (Mat_HYPRE *)A->data;
1179: hB = (Mat_HYPRE *)B->data;
1180: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1181: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1182: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hB->ij, &type));
1183: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1184: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&Aparcsr));
1185: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hB->ij, (void **)&Bparcsr));
1186: PetscCall(MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr));
1187: PetscCall(MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D));
1189: /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1190: product = C->product; /* save it from MatHeaderReplace() */
1191: C->product = NULL;
1192: PetscCall(MatHeaderReplace(C, &D));
1193: C->product = product;
1194: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1195: C->ops->productnumeric = MatProductNumeric_AB;
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D)
1200: {
1201: Mat E;
1202: hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL;
1204: PetscFunctionBegin;
1205: PetscCall(MatAIJGetParCSR_Private(A, &hA));
1206: PetscCall(MatAIJGetParCSR_Private(B, &hB));
1207: PetscCall(MatAIJGetParCSR_Private(C, &hC));
1208: PetscCall(MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC));
1209: PetscCall(MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E));
1210: PetscCall(MatHeaderMerge(D, &E));
1211: PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1212: PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1213: PetscCall(MatAIJRestoreParCSR_Private(C, &hC));
1214: PetscFunctionReturn(PETSC_SUCCESS);
1215: }
1217: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
1218: {
1219: PetscFunctionBegin;
1220: PetscCall(MatSetType(D, MATAIJ));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1225: {
1226: PetscFunctionBegin;
1227: C->ops->productnumeric = MatProductNumeric_AB;
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1232: {
1233: Mat_Product *product = C->product;
1234: PetscBool Ahypre;
1236: PetscFunctionBegin;
1237: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre));
1238: if (Ahypre) { /* A is a Hypre matrix */
1239: PetscCall(MatSetType(C, MATHYPRE));
1240: C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1241: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1242: PetscFunctionReturn(PETSC_SUCCESS);
1243: }
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1248: {
1249: PetscFunctionBegin;
1250: C->ops->productnumeric = MatProductNumeric_PtAP;
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1255: {
1256: Mat_Product *product = C->product;
1257: PetscBool flg;
1258: PetscInt type = 0;
1259: const char *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"};
1260: PetscInt ntype = 4;
1261: Mat A = product->A;
1262: PetscBool Ahypre;
1264: PetscFunctionBegin;
1265: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre));
1266: if (Ahypre) { /* A is a Hypre matrix */
1267: PetscCall(MatSetType(C, MATHYPRE));
1268: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1269: C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1274: /* Get runtime option */
1275: if (product->api_user) {
1276: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat");
1277: PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg));
1278: PetscOptionsEnd();
1279: } else {
1280: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat");
1281: PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg));
1282: PetscOptionsEnd();
1283: }
1285: if (type == 0 || type == 1 || type == 2) {
1286: PetscCall(MatSetType(C, MATAIJ));
1287: } else if (type == 3) {
1288: PetscCall(MatSetType(C, MATHYPRE));
1289: } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported");
1290: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1291: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }
1295: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1296: {
1297: Mat_Product *product = C->product;
1299: PetscFunctionBegin;
1300: switch (product->type) {
1301: case MATPRODUCT_AB:
1302: PetscCall(MatProductSetFromOptions_HYPRE_AB(C));
1303: break;
1304: case MATPRODUCT_PtAP:
1305: PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C));
1306: break;
1307: default:
1308: break;
1309: }
1310: PetscFunctionReturn(PETSC_SUCCESS);
1311: }
1313: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1314: {
1315: PetscFunctionBegin;
1316: PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE));
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1320: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1321: {
1322: PetscFunctionBegin;
1323: PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE));
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1328: {
1329: PetscFunctionBegin;
1330: if (y != z) PetscCall(VecCopy(y, z));
1331: PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE));
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1336: {
1337: PetscFunctionBegin;
1338: if (y != z) PetscCall(VecCopy(y, z));
1339: PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE));
1340: PetscFunctionReturn(PETSC_SUCCESS);
1341: }
1343: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1344: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1345: {
1346: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1347: hypre_ParCSRMatrix *parcsr;
1348: hypre_ParVector *hx, *hy;
1350: PetscFunctionBegin;
1351: if (trans) {
1352: PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x));
1353: if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y));
1354: else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y));
1355: PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->b->ij, (void **)&hx));
1356: PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->x->ij, (void **)&hy));
1357: } else {
1358: PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x));
1359: if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y));
1360: else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y));
1361: PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->x->ij, (void **)&hx));
1362: PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->b->ij, (void **)&hy));
1363: }
1364: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1365: if (trans) {
1366: PetscCallHYPRE(hypre_ParCSRMatrixMatvecT(a, parcsr, hx, b, hy));
1367: } else {
1368: PetscCallHYPRE(hypre_ParCSRMatrixMatvec(a, parcsr, hx, b, hy));
1369: }
1370: PetscCall(VecHYPRE_IJVectorPopVec(hA->x));
1371: PetscCall(VecHYPRE_IJVectorPopVec(hA->b));
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }
1375: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1376: {
1377: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1379: PetscFunctionBegin;
1380: PetscCall(VecHYPRE_IJVectorDestroy(&hA->x));
1381: PetscCall(VecHYPRE_IJVectorDestroy(&hA->b));
1382: PetscCall(MatHYPRE_DestroyCOOMat(A)); /* must be called before destroying the individual CSR */
1383: if (hA->ij) {
1384: if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1385: PetscCallHYPRE(HYPRE_IJMatrixDestroy(hA->ij));
1386: }
1387: if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm));
1389: PetscCall(MatStashDestroy_Private(&A->stash));
1390: PetscCall(PetscFree(hA->array));
1391: if (hA->rows_d) PetscStackCallExternalVoid("hypre_Free", hypre_Free(hA->rows_d, HYPRE_MEMORY_DEVICE));
1393: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL));
1394: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL));
1395: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL));
1396: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL));
1397: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", NULL));
1398: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", NULL));
1399: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_hypre_C", NULL));
1400: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", NULL));
1401: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL));
1402: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL));
1403: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
1404: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
1405: PetscCall(PetscFree(A->data));
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1410: {
1411: PetscFunctionBegin;
1412: if (!A->preallocated) PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1417: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1418: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1419: {
1420: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1421: HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
1423: PetscFunctionBegin;
1424: A->boundtocpu = bind;
1425: if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1426: hypre_ParCSRMatrix *parcsr;
1427: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1428: PetscCallHYPRE(hypre_ParCSRMatrixMigrate(parcsr, hmem));
1429: }
1430: if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind));
1431: if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind));
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1434: #endif
1436: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1437: {
1438: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1439: PetscMPIInt n;
1440: PetscInt i, j, rstart, ncols, flg;
1441: PetscInt *row, *col;
1442: PetscScalar *val;
1444: PetscFunctionBegin;
1445: PetscCheck(mode != MAT_FLUSH_ASSEMBLY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");
1447: if (!A->nooffprocentries) {
1448: while (1) {
1449: PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
1450: if (!flg) break;
1452: for (i = 0; i < n;) {
1453: /* Now identify the consecutive vals belonging to the same row */
1454: for (j = i, rstart = row[j]; j < n; j++) {
1455: if (row[j] != rstart) break;
1456: }
1457: if (j < n) ncols = j - i;
1458: else ncols = n - i;
1459: /* Now assemble all these values with a single function call */
1460: PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode));
1462: i = j;
1463: }
1464: }
1465: PetscCall(MatStashScatterEnd_Private(&A->stash));
1466: }
1468: PetscCallHYPRE(HYPRE_IJMatrixAssemble(hA->ij));
1469: /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1470: /* If the option MAT_SORTED_FULL is set to true, the indices and values can be passed to hypre directly, so we don't need the aux_matrix */
1471: if (!A->sortedfull) {
1472: hypre_AuxParCSRMatrix *aux_matrix;
1474: /* call destroy just to make sure we do not leak anything */
1475: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1476: PetscCallHYPRE(hypre_AuxParCSRMatrixDestroy(aux_matrix));
1477: hypre_IJMatrixTranslator(hA->ij) = NULL;
1479: /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1480: PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1481: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1482: if (aux_matrix) {
1483: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1484: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1485: PetscCallHYPRE(hypre_AuxParCSRMatrixInitialize(aux_matrix));
1486: #else
1487: PetscCallHYPRE(hypre_AuxParCSRMatrixInitialize_v2(aux_matrix, HYPRE_MEMORY_HOST));
1488: #endif
1489: }
1490: }
1491: {
1492: hypre_ParCSRMatrix *parcsr;
1494: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1495: if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallHYPRE(hypre_MatvecCommPkgCreate(parcsr));
1496: }
1497: if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x));
1498: if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b));
1499: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1500: PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu));
1501: #endif
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1506: {
1507: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1509: PetscFunctionBegin;
1510: PetscCheck(hA->array_available, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Temporary space is in use");
1512: if (hA->array_size >= size) {
1513: *array = hA->array;
1514: } else {
1515: PetscCall(PetscFree(hA->array));
1516: hA->array_size = size;
1517: PetscCall(PetscMalloc(hA->array_size, &hA->array));
1518: *array = hA->array;
1519: }
1521: hA->array_available = PETSC_FALSE;
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1526: {
1527: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1529: PetscFunctionBegin;
1530: *array = NULL;
1531: hA->array_available = PETSC_TRUE;
1532: PetscFunctionReturn(PETSC_SUCCESS);
1533: }
1535: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1536: {
1537: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1538: PetscScalar *vals = (PetscScalar *)v;
1539: HYPRE_Complex *sscr;
1540: PetscInt *cscr[2];
1541: PetscInt i, nzc;
1542: PetscInt rst = A->rmap->rstart, ren = A->rmap->rend;
1543: void *array = NULL;
1545: PetscFunctionBegin;
1546: PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array));
1547: cscr[0] = (PetscInt *)array;
1548: cscr[1] = ((PetscInt *)array) + nc;
1549: sscr = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2);
1550: for (i = 0, nzc = 0; i < nc; i++) {
1551: if (cols[i] >= 0) {
1552: cscr[0][nzc] = cols[i];
1553: cscr[1][nzc++] = i;
1554: }
1555: }
1556: if (!nzc) {
1557: PetscCall(MatRestoreArray_HYPRE(A, &array));
1558: PetscFunctionReturn(PETSC_SUCCESS);
1559: }
1561: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1562: if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1563: hypre_ParCSRMatrix *parcsr;
1565: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij,(void**)&parcsr));
1566: PetscCallHYPRE(hypre_ParCSRMatrixMigrate(parcsr, HYPRE_MEMORY_HOST));
1567: }
1568: #endif
1570: if (ins == ADD_VALUES) {
1571: for (i = 0; i < nr; i++) {
1572: if (rows[i] >= 0) {
1573: PetscInt j;
1574: HYPRE_Int hnc = (HYPRE_Int)nzc;
1576: if (!nzc) continue;
1577: /* nonlocal values */
1578: if (rows[i] < rst || rows[i] >= ren) {
1579: PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1580: if (hA->donotstash) continue;
1581: }
1582: PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1583: for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1584: PetscCallHYPRE(HYPRE_IJMatrixAddToValues(hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr));
1585: }
1586: vals += nc;
1587: }
1588: } else { /* INSERT_VALUES */
1589: for (i = 0; i < nr; i++) {
1590: if (rows[i] >= 0) {
1591: PetscInt j;
1592: HYPRE_Int hnc = (HYPRE_Int)nzc;
1594: if (!nzc) continue;
1595: PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1596: for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1597: /* nonlocal values */
1598: if (rows[i] < rst || rows[i] >= ren) {
1599: PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1600: if (!hA->donotstash) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE));
1601: }
1602: /* local values */
1603: else
1604: PetscCallHYPRE(HYPRE_IJMatrixSetValues(hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr));
1605: }
1606: vals += nc;
1607: }
1608: }
1610: PetscCall(MatRestoreArray_HYPRE(A, &array));
1611: PetscFunctionReturn(PETSC_SUCCESS);
1612: }
1614: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1615: {
1616: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1617: HYPRE_Int *hdnnz, *honnz;
1618: PetscInt i, rs, re, cs, ce, bs;
1619: PetscMPIInt size;
1621: PetscFunctionBegin;
1622: PetscCall(PetscLayoutSetUp(A->rmap));
1623: PetscCall(PetscLayoutSetUp(A->cmap));
1624: rs = A->rmap->rstart;
1625: re = A->rmap->rend;
1626: cs = A->cmap->rstart;
1627: ce = A->cmap->rend;
1628: if (!hA->ij) {
1629: PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rs, re - 1, cs, ce - 1, &hA->ij));
1630: PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
1631: } else {
1632: HYPRE_BigInt hrs, hre, hcs, hce;
1633: PetscCallHYPRE(HYPRE_IJMatrixGetLocalRange(hA->ij, &hrs, &hre, &hcs, &hce));
1634: PetscCheck(hre - hrs + 1 == re - rs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local rows: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hrs, hre + 1, rs, re);
1635: PetscCheck(hce - hcs + 1 == ce - cs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local cols: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hcs, hce + 1, cs, ce);
1636: }
1637: PetscCall(MatHYPRE_DestroyCOOMat(A));
1638: PetscCall(MatGetBlockSize(A, &bs));
1639: if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
1640: if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;
1642: if (!dnnz) {
1643: PetscCall(PetscMalloc1(A->rmap->n, &hdnnz));
1644: for (i = 0; i < A->rmap->n; i++) hdnnz[i] = (HYPRE_Int)dnz;
1645: } else {
1646: hdnnz = (HYPRE_Int *)dnnz;
1647: }
1648: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1649: if (size > 1) {
1650: hypre_AuxParCSRMatrix *aux_matrix;
1651: if (!onnz) {
1652: PetscCall(PetscMalloc1(A->rmap->n, &honnz));
1653: for (i = 0; i < A->rmap->n; i++) honnz[i] = (HYPRE_Int)onz;
1654: } else honnz = (HYPRE_Int *)onnz;
1655: /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1656: they assume the user will input the entire row values, properly sorted
1657: In PETSc, we don't make such an assumption and set this flag to 1,
1658: unless the option MAT_SORTED_FULL is set to true.
1659: Also, to avoid possible memory leaks, we destroy and recreate the translator
1660: This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1661: the IJ matrix for us */
1662: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1663: hypre_AuxParCSRMatrixDestroy(aux_matrix);
1664: hypre_IJMatrixTranslator(hA->ij) = NULL;
1665: PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(hA->ij, hdnnz, honnz));
1666: aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1667: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !A->sortedfull;
1668: } else {
1669: honnz = NULL;
1670: PetscCallHYPRE(HYPRE_IJMatrixSetRowSizes(hA->ij, hdnnz));
1671: }
1673: /* reset assembled flag and call the initialize method */
1674: hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1675: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1676: PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1677: #else
1678: PetscCallHYPRE(HYPRE_IJMatrixInitialize_v2(hA->ij, HYPRE_MEMORY_HOST));
1679: #endif
1680: if (!dnnz) PetscCall(PetscFree(hdnnz));
1681: if (!onnz && honnz) PetscCall(PetscFree(honnz));
1682: /* Match AIJ logic */
1683: A->preallocated = PETSC_TRUE;
1684: A->assembled = PETSC_FALSE;
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@C
1689: MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1691: Collective
1693: Input Parameters:
1694: + A - the matrix
1695: . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix
1696: (same value is used for all local rows)
1697: . dnnz - array containing the number of nonzeros in the various rows of the
1698: DIAGONAL portion of the local submatrix (possibly different for each row)
1699: or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1700: The size of this array is equal to the number of local rows, i.e `m`.
1701: For matrices that will be factored, you must leave room for (and set)
1702: the diagonal entry even if it is zero.
1703: . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1704: submatrix (same value is used for all local rows).
1705: - onnz - array containing the number of nonzeros in the various rows of the
1706: OFF-DIAGONAL portion of the local submatrix (possibly different for
1707: each row) or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1708: structure. The size of this array is equal to the number
1709: of local rows, i.e `m`.
1711: Level: intermediate
1713: Note:
1714: If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, `onz` and `onnz` are ignored.
1716: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`, `MATAIJ`
1717: @*/
1718: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1719: {
1720: PetscFunctionBegin;
1723: PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz));
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }
1727: /*@C
1728: MatCreateFromParCSR - Creates a `Mat` from a `hypre_ParCSRMatrix`
1730: Collective
1732: Input Parameters:
1733: + parcsr - the pointer to the `hypre_ParCSRMatrix`
1734: . mtype - matrix type to be created. Currently `MATAIJ`, `MATIS` and `MATHYPRE` are supported.
1735: - copymode - PETSc copying options, see `PetscCopyMode`
1737: Output Parameter:
1738: . A - the matrix
1740: Level: intermediate
1742: .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1743: @*/
1744: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A)
1745: {
1746: Mat T;
1747: Mat_HYPRE *hA;
1748: MPI_Comm comm;
1749: PetscInt rstart, rend, cstart, cend, M, N;
1750: PetscBool isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis;
1752: PetscFunctionBegin;
1753: comm = hypre_ParCSRMatrixComm(parcsr);
1754: PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij));
1755: PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl));
1756: PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij));
1757: PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
1758: PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp));
1759: PetscCall(PetscStrcmp(mtype, MATIS, &isis));
1760: isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1761: /* TODO */
1762: PetscCheck(isaij || ishyp || isis, comm, PETSC_ERR_SUP, "Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s", mtype, MATAIJ, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, MATIS, MATHYPRE);
1763: /* access ParCSRMatrix */
1764: rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1765: rend = hypre_ParCSRMatrixLastRowIndex(parcsr);
1766: cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1767: cend = hypre_ParCSRMatrixLastColDiag(parcsr);
1768: M = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1769: N = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1771: /* create PETSc matrix with MatHYPRE */
1772: PetscCall(MatCreate(comm, &T));
1773: PetscCall(MatSetSizes(T, PetscMax(rend - rstart + 1, 0), PetscMax(cend - cstart + 1, 0), M, N));
1774: PetscCall(MatSetType(T, MATHYPRE));
1775: hA = (Mat_HYPRE *)T->data;
1777: /* create HYPRE_IJMatrix */
1778: PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rstart, rend, cstart, cend, &hA->ij));
1779: PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
1781: /* create new ParCSR object if needed */
1782: if (ishyp && copymode == PETSC_COPY_VALUES) {
1783: hypre_ParCSRMatrix *new_parcsr;
1784: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
1785: hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd;
1787: new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1788: hdiag = hypre_ParCSRMatrixDiag(parcsr);
1789: hoffd = hypre_ParCSRMatrixOffd(parcsr);
1790: ndiag = hypre_ParCSRMatrixDiag(new_parcsr);
1791: noffd = hypre_ParCSRMatrixOffd(new_parcsr);
1792: PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag)));
1793: PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd)));
1794: #else
1795: new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1);
1796: #endif
1797: parcsr = new_parcsr;
1798: copymode = PETSC_OWN_POINTER;
1799: }
1801: /* set ParCSR object */
1802: hypre_IJMatrixObject(hA->ij) = parcsr;
1803: T->preallocated = PETSC_TRUE;
1805: /* set assembled flag */
1806: hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1807: #if 0
1808: PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1809: #endif
1810: if (ishyp) {
1811: PetscMPIInt myid = 0;
1813: /* make sure we always have row_starts and col_starts available */
1814: if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid));
1815: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1816: if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1817: PetscLayout map;
1819: PetscCall(MatGetLayouts(T, NULL, &map));
1820: PetscCall(PetscLayoutSetUp(map));
1821: hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1822: }
1823: if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1824: PetscLayout map;
1826: PetscCall(MatGetLayouts(T, &map, NULL));
1827: PetscCall(PetscLayoutSetUp(map));
1828: hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1829: }
1830: #endif
1831: /* prevent from freeing the pointer */
1832: if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1833: *A = T;
1834: PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE));
1835: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
1836: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
1837: } else if (isaij) {
1838: if (copymode != PETSC_OWN_POINTER) {
1839: /* prevent from freeing the pointer */
1840: hA->inner_free = PETSC_FALSE;
1841: PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A));
1842: PetscCall(MatDestroy(&T));
1843: } else { /* AIJ return type with PETSC_OWN_POINTER */
1844: PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T));
1845: *A = T;
1846: }
1847: } else if (isis) {
1848: PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A));
1849: if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1850: PetscCall(MatDestroy(&T));
1851: }
1852: PetscFunctionReturn(PETSC_SUCCESS);
1853: }
1855: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1856: {
1857: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1858: HYPRE_Int type;
1860: PetscFunctionBegin;
1861: PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present");
1862: PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1863: PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1864: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)parcsr));
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: /*@C
1869: MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1871: Not Collective, No Fortran Support
1873: Input Parameter:
1874: . A - the `MATHYPRE` object
1876: Output Parameter:
1877: . parcsr - the pointer to the `hypre_ParCSRMatrix`
1879: Level: intermediate
1881: .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1882: @*/
1883: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1884: {
1885: PetscFunctionBegin;
1888: PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr));
1889: PetscFunctionReturn(PETSC_SUCCESS);
1890: }
1892: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1893: {
1894: hypre_ParCSRMatrix *parcsr;
1895: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1896: hypre_CSRMatrix *ha;
1897: #endif
1898: HYPRE_Complex hs;
1900: PetscFunctionBegin;
1901: PetscCall(PetscHYPREScalarCast(s, &hs));
1902: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1903: #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0)
1904: PetscCallHYPRE(hypre_ParCSRMatrixScale(parcsr, hs));
1905: #else /* diagonal part */
1906: ha = hypre_ParCSRMatrixDiag(parcsr);
1907: if (ha) {
1908: PetscInt size, i;
1909: HYPRE_Int *ii;
1910: HYPRE_Complex *a;
1912: size = hypre_CSRMatrixNumRows(ha);
1913: a = hypre_CSRMatrixData(ha);
1914: ii = hypre_CSRMatrixI(ha);
1915: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1916: }
1917: /* off-diagonal part */
1918: ha = hypre_ParCSRMatrixOffd(parcsr);
1919: if (ha) {
1920: PetscInt size, i;
1921: HYPRE_Int *ii;
1922: HYPRE_Complex *a;
1924: size = hypre_CSRMatrixNumRows(ha);
1925: a = hypre_CSRMatrixData(ha);
1926: ii = hypre_CSRMatrixI(ha);
1927: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1928: }
1929: #endif
1930: PetscFunctionReturn(PETSC_SUCCESS);
1931: }
1933: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1934: {
1935: hypre_ParCSRMatrix *parcsr;
1936: HYPRE_Int *lrows;
1937: PetscInt rst, ren, i;
1939: PetscFunctionBegin;
1940: PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented");
1941: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1942: PetscCall(PetscMalloc1(numRows, &lrows));
1943: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1944: for (i = 0; i < numRows; i++) {
1945: PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported");
1946: lrows[i] = (HYPRE_Int)(rows[i] - rst);
1947: }
1948: PetscCallHYPRE(hypre_ParCSRMatrixEliminateRowsCols(parcsr, (HYPRE_Int)numRows, lrows));
1949: PetscCall(PetscFree(lrows));
1950: PetscFunctionReturn(PETSC_SUCCESS);
1951: }
1953: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1954: {
1955: PetscFunctionBegin;
1956: if (ha) {
1957: HYPRE_Int *ii, size;
1958: HYPRE_Complex *a;
1960: size = hypre_CSRMatrixNumRows(ha);
1961: a = hypre_CSRMatrixData(ha);
1962: ii = hypre_CSRMatrixI(ha);
1964: if (a) PetscCall(PetscArrayzero(a, ii[size]));
1965: }
1966: PetscFunctionReturn(PETSC_SUCCESS);
1967: }
1969: static PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1970: {
1971: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1973: PetscFunctionBegin;
1974: if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1975: PetscCallHYPRE(HYPRE_IJMatrixSetConstantValues(hA->ij, 0.0));
1976: } else {
1977: hypre_ParCSRMatrix *parcsr;
1979: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1980: PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr)));
1981: PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr)));
1982: }
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1986: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag)
1987: {
1988: PetscInt ii;
1989: HYPRE_Int *i, *j;
1990: HYPRE_Complex *a;
1992: PetscFunctionBegin;
1993: if (!hA) PetscFunctionReturn(PETSC_SUCCESS);
1995: i = hypre_CSRMatrixI(hA);
1996: j = hypre_CSRMatrixJ(hA);
1997: a = hypre_CSRMatrixData(hA);
1998: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1999: if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hA)) {
2000: #if defined(HYPRE_USING_CUDA)
2001: PetscCall(MatZeroRows_CUDA(N, rows, i, j, a, diag));
2002: #elif defined(HYPRE_USING_HIP)
2003: PetscCall(MatZeroRows_HIP(N, rows, i, j, a, diag));
2004: #elif defined(PETSC_HAVE_KOKKOS)
2005: PetscCall(MatZeroRows_Kokkos(N, rows, i, j, a, diag));
2006: #else
2007: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for MatZeroRows on a hypre matrix in this memory location");
2008: #endif
2009: } else
2010: #endif
2011: {
2012: for (ii = 0; ii < N; ii++) {
2013: HYPRE_Int jj, ibeg, iend, irow;
2015: irow = (HYPRE_Int)rows[ii];
2016: ibeg = i[irow];
2017: iend = i[irow + 1];
2018: for (jj = ibeg; jj < iend; jj++)
2019: if (j[jj] == irow) a[jj] = diag;
2020: else a[jj] = 0.0;
2021: }
2022: }
2023: PetscFunctionReturn(PETSC_SUCCESS);
2024: }
2026: static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2027: {
2028: hypre_ParCSRMatrix *parcsr;
2029: PetscInt *lrows, len, *lrows2;
2030: HYPRE_Complex hdiag;
2032: PetscFunctionBegin;
2033: PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
2034: PetscCall(PetscHYPREScalarCast(diag, &hdiag));
2035: /* retrieve the internal matrix */
2036: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2037: /* get locally owned rows */
2038: PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));
2040: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2041: if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hypre_ParCSRMatrixDiag(parcsr))) {
2042: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2043: PetscInt m;
2044: PetscCall(MatGetLocalSize(A, &m, NULL));
2045: if (!hA->rows_d) {
2046: hA->rows_d = hypre_TAlloc(PetscInt, m, HYPRE_MEMORY_DEVICE);
2047: if (m) PetscCheck(hA->rows_d, PETSC_COMM_SELF, PETSC_ERR_MEM, "HYPRE_TAlloc failed");
2048: }
2049: PetscCheck(len <= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many rows in rows[]");
2050: PetscStackCallExternalVoid("hypre_Memcpy", hypre_Memcpy(hA->rows_d, lrows, sizeof(PetscInt) * len, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST));
2051: lrows2 = hA->rows_d;
2052: } else
2053: #endif
2054: {
2055: lrows2 = lrows;
2056: }
2058: /* zero diagonal part */
2059: PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows2, hdiag));
2060: /* zero off-diagonal part */
2061: PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows2, 0.0));
2063: PetscCall(PetscFree(lrows));
2064: PetscFunctionReturn(PETSC_SUCCESS);
2065: }
2067: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode)
2068: {
2069: PetscFunctionBegin;
2070: if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
2072: PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
2073: PetscFunctionReturn(PETSC_SUCCESS);
2074: }
2076: static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2077: {
2078: hypre_ParCSRMatrix *parcsr;
2079: HYPRE_Int hnz;
2081: PetscFunctionBegin;
2082: /* retrieve the internal matrix */
2083: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2084: /* call HYPRE API */
2085: PetscCallHYPRE(HYPRE_ParCSRMatrixGetRow(parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v));
2086: if (nz) *nz = (PetscInt)hnz;
2087: PetscFunctionReturn(PETSC_SUCCESS);
2088: }
2090: static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2091: {
2092: hypre_ParCSRMatrix *parcsr;
2093: HYPRE_Int hnz;
2095: PetscFunctionBegin;
2096: /* retrieve the internal matrix */
2097: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2098: /* call HYPRE API */
2099: hnz = nz ? (HYPRE_Int)(*nz) : 0;
2100: PetscCallHYPRE(HYPRE_ParCSRMatrixRestoreRow(parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v));
2101: PetscFunctionReturn(PETSC_SUCCESS);
2102: }
2104: static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
2105: {
2106: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2107: PetscInt i;
2109: PetscFunctionBegin;
2110: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
2111: /* Ignore negative row indices
2112: * And negative column indices should be automatically ignored in hypre
2113: * */
2114: for (i = 0; i < m; i++) {
2115: if (idxm[i] >= 0) {
2116: HYPRE_Int hn = (HYPRE_Int)n;
2117: PetscCallHYPRE(HYPRE_IJMatrixGetValues(hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n)));
2118: }
2119: }
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg)
2124: {
2125: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2127: PetscFunctionBegin;
2128: switch (op) {
2129: case MAT_NO_OFF_PROC_ENTRIES:
2130: if (flg) PetscCallHYPRE(HYPRE_IJMatrixSetMaxOffProcElmts(hA->ij, 0));
2131: break;
2132: case MAT_IGNORE_OFF_PROC_ENTRIES:
2133: hA->donotstash = flg;
2134: break;
2135: default:
2136: break;
2137: }
2138: PetscFunctionReturn(PETSC_SUCCESS);
2139: }
2141: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2142: {
2143: PetscViewerFormat format;
2145: PetscFunctionBegin;
2146: PetscCall(PetscViewerGetFormat(view, &format));
2147: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
2148: if (format != PETSC_VIEWER_NATIVE) {
2149: Mat B;
2150: hypre_ParCSRMatrix *parcsr;
2151: PetscErrorCode (*mview)(Mat, PetscViewer) = NULL;
2153: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2154: PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B));
2155: PetscCall(MatGetOperation(B, MATOP_VIEW, (PetscErrorCodeFn **)&mview));
2156: PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation");
2157: PetscCall((*mview)(B, view));
2158: PetscCall(MatDestroy(&B));
2159: } else {
2160: Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2161: PetscMPIInt size;
2162: PetscBool isascii;
2163: const char *filename;
2165: /* HYPRE uses only text files */
2166: PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
2167: PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name);
2168: PetscCall(PetscViewerFileGetName(view, &filename));
2169: PetscCallHYPRE(HYPRE_IJMatrixPrint(hA->ij, filename));
2170: PetscCallMPI(MPI_Comm_size(hA->comm, &size));
2171: if (size > 1) {
2172: PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1));
2173: } else {
2174: PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0));
2175: }
2176: }
2177: PetscFunctionReturn(PETSC_SUCCESS);
2178: }
2180: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2181: {
2182: hypre_ParCSRMatrix *acsr, *bcsr;
2184: PetscFunctionBegin;
2185: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2186: PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr));
2187: PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr));
2188: PetscCallHYPRE(hypre_ParCSRMatrixCopy(acsr, bcsr, 1));
2189: PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2190: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2191: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2192: } else {
2193: PetscCall(MatCopy_Basic(A, B, str));
2194: }
2195: PetscFunctionReturn(PETSC_SUCCESS);
2196: }
2198: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2199: {
2200: hypre_ParCSRMatrix *parcsr;
2201: hypre_CSRMatrix *dmat;
2202: HYPRE_Complex *a;
2203: PetscBool cong;
2205: PetscFunctionBegin;
2206: PetscCall(MatHasCongruentLayouts(A, &cong));
2207: PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns");
2208: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2209: dmat = hypre_ParCSRMatrixDiag(parcsr);
2210: if (dmat) {
2211: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2212: HYPRE_MemoryLocation mem = hypre_CSRMatrixMemoryLocation(dmat);
2213: #else
2214: HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
2215: #endif
2217: if (mem != HYPRE_MEMORY_HOST) PetscCall(VecGetArrayWriteAndMemType(d, (PetscScalar **)&a, NULL));
2218: else PetscCall(VecGetArrayWrite(d, (PetscScalar **)&a));
2219: hypre_CSRMatrixExtractDiagonal(dmat, a, 0);
2220: if (mem != HYPRE_MEMORY_HOST) PetscCall(VecRestoreArrayWriteAndMemType(d, (PetscScalar **)&a));
2221: else PetscCall(VecRestoreArrayWrite(d, (PetscScalar **)&a));
2222: }
2223: PetscFunctionReturn(PETSC_SUCCESS);
2224: }
2226: #include <petscblaslapack.h>
2228: static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str)
2229: {
2230: PetscFunctionBegin;
2231: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2232: {
2233: Mat B;
2234: hypre_ParCSRMatrix *x, *y, *z;
2236: PetscCall(MatHYPREGetParCSR(Y, &y));
2237: PetscCall(MatHYPREGetParCSR(X, &x));
2238: PetscCallHYPRE(hypre_ParCSRMatrixAdd(1.0, y, 1.0, x, &z));
2239: PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B));
2240: PetscCall(MatHeaderMerge(Y, &B));
2241: }
2242: #else
2243: if (str == SAME_NONZERO_PATTERN) {
2244: hypre_ParCSRMatrix *x, *y;
2245: hypre_CSRMatrix *xloc, *yloc;
2246: PetscInt xnnz, ynnz;
2247: HYPRE_Complex *xarr, *yarr;
2248: PetscBLASInt one = 1, bnz;
2250: PetscCall(MatHYPREGetParCSR(Y, &y));
2251: PetscCall(MatHYPREGetParCSR(X, &x));
2253: /* diagonal block */
2254: xloc = hypre_ParCSRMatrixDiag(x);
2255: yloc = hypre_ParCSRMatrixDiag(y);
2256: xnnz = 0;
2257: ynnz = 0;
2258: xarr = NULL;
2259: yarr = NULL;
2260: if (xloc) {
2261: xarr = hypre_CSRMatrixData(xloc);
2262: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2263: }
2264: if (yloc) {
2265: yarr = hypre_CSRMatrixData(yloc);
2266: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2267: }
2268: PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2269: PetscCall(PetscBLASIntCast(xnnz, &bnz));
2270: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2272: /* off-diagonal block */
2273: xloc = hypre_ParCSRMatrixOffd(x);
2274: yloc = hypre_ParCSRMatrixOffd(y);
2275: xnnz = 0;
2276: ynnz = 0;
2277: xarr = NULL;
2278: yarr = NULL;
2279: if (xloc) {
2280: xarr = hypre_CSRMatrixData(xloc);
2281: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2282: }
2283: if (yloc) {
2284: yarr = hypre_CSRMatrixData(yloc);
2285: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2286: }
2287: PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2288: PetscCall(PetscBLASIntCast(xnnz, &bnz));
2289: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2290: } else if (str == SUBSET_NONZERO_PATTERN) {
2291: PetscCall(MatAXPY_Basic(Y, a, X, str));
2292: } else {
2293: Mat B;
2295: PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
2296: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2297: PetscCall(MatHeaderReplace(Y, &B));
2298: }
2299: #endif
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B)
2304: {
2305: hypre_ParCSRMatrix *parcsr = NULL;
2306: PetscCopyMode cpmode;
2307: Mat_HYPRE *hA;
2309: PetscFunctionBegin;
2310: PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2311: if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2312: parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
2313: cpmode = PETSC_OWN_POINTER;
2314: } else {
2315: cpmode = PETSC_COPY_VALUES;
2316: }
2317: PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B));
2318: hA = (Mat_HYPRE *)A->data;
2319: if (hA->cooMat) {
2320: Mat_HYPRE *hB = (Mat_HYPRE *)((*B)->data);
2321: op = (op == MAT_DO_NOT_COPY_VALUES) ? op : MAT_COPY_VALUES;
2322: /* Cannot simply increase the reference count of hA->cooMat, since B needs to share cooMat's data array */
2323: PetscCall(MatDuplicate(hA->cooMat, op, &hB->cooMat));
2324: PetscCall(MatHYPRE_AttachCOOMat(*B));
2325: }
2326: PetscFunctionReturn(PETSC_SUCCESS);
2327: }
2329: static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[])
2330: {
2331: Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
2333: PetscFunctionBegin;
2334: /* Build an agent matrix cooMat with AIJ format
2335: It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2336: */
2337: PetscCall(MatHYPRE_CreateCOOMat(mat));
2338: PetscCall(MatSetOption(hmat->cooMat, MAT_IGNORE_OFF_PROC_ENTRIES, hmat->donotstash));
2339: PetscCall(MatSetOption(hmat->cooMat, MAT_NO_OFF_PROC_ENTRIES, mat->nooffprocentries));
2341: /* MatSetPreallocationCOO_SeqAIJ and MatSetPreallocationCOO_MPIAIJ uses this specific
2342: name to automatically put the diagonal entries first */
2343: PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
2344: PetscCall(MatSetPreallocationCOO(hmat->cooMat, coo_n, coo_i, coo_j));
2345: hmat->cooMat->assembled = PETSC_TRUE;
2347: /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */
2348: PetscCall(MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE));
2349: PetscCall(MatHYPRE_CreateFromMat(hmat->cooMat, hmat)); /* Create hmat->ij and preallocate it */
2350: PetscCall(MatHYPRE_IJMatrixCopyIJ(hmat->cooMat, hmat->ij)); /* Copy A's (i,j) to hmat->ij */
2352: mat->preallocated = PETSC_TRUE;
2353: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2354: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */
2356: /* Attach cooMat to mat */
2357: PetscCall(MatHYPRE_AttachCOOMat(mat));
2358: PetscFunctionReturn(PETSC_SUCCESS);
2359: }
2361: static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode)
2362: {
2363: Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
2365: PetscFunctionBegin;
2366: PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
2367: PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode));
2368: PetscCall(MatViewFromOptions(hmat->cooMat, (PetscObject)mat, "-cooMat_view"));
2369: PetscFunctionReturn(PETSC_SUCCESS);
2370: }
2372: static PetscErrorCode MatGetCurrentMemType_HYPRE(Mat A, PetscMemType *m)
2373: {
2374: PetscBool petsconcpu;
2376: PetscFunctionBegin;
2377: PetscCall(MatBoundToCPU(A, &petsconcpu));
2378: *m = petsconcpu ? PETSC_MEMTYPE_HOST : PETSC_MEMTYPE_DEVICE;
2379: PetscFunctionReturn(PETSC_SUCCESS);
2380: }
2382: /*MC
2383: MATHYPRE - "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2384: based on the hypre IJ interface.
2386: Level: intermediate
2388: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatHYPRESetPreallocation`
2389: M*/
2390: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2391: {
2392: Mat_HYPRE *hB;
2393: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2394: HYPRE_MemoryLocation memory_location;
2395: #endif
2397: PetscFunctionBegin;
2398: PetscHYPREInitialize();
2399: PetscCall(PetscNew(&hB));
2401: hB->inner_free = PETSC_TRUE;
2402: hB->array_available = PETSC_TRUE;
2404: B->data = (void *)hB;
2406: PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps)));
2407: B->ops->mult = MatMult_HYPRE;
2408: B->ops->multtranspose = MatMultTranspose_HYPRE;
2409: B->ops->multadd = MatMultAdd_HYPRE;
2410: B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
2411: B->ops->setup = MatSetUp_HYPRE;
2412: B->ops->destroy = MatDestroy_HYPRE;
2413: B->ops->assemblyend = MatAssemblyEnd_HYPRE;
2414: B->ops->assemblybegin = MatAssemblyBegin_HYPRE;
2415: B->ops->setvalues = MatSetValues_HYPRE;
2416: B->ops->scale = MatScale_HYPRE;
2417: B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
2418: B->ops->zeroentries = MatZeroEntries_HYPRE;
2419: B->ops->zerorows = MatZeroRows_HYPRE;
2420: B->ops->getrow = MatGetRow_HYPRE;
2421: B->ops->restorerow = MatRestoreRow_HYPRE;
2422: B->ops->getvalues = MatGetValues_HYPRE;
2423: B->ops->setoption = MatSetOption_HYPRE;
2424: B->ops->duplicate = MatDuplicate_HYPRE;
2425: B->ops->copy = MatCopy_HYPRE;
2426: B->ops->view = MatView_HYPRE;
2427: B->ops->getdiagonal = MatGetDiagonal_HYPRE;
2428: B->ops->axpy = MatAXPY_HYPRE;
2429: B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2430: B->ops->getcurrentmemtype = MatGetCurrentMemType_HYPRE;
2431: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2432: B->ops->bindtocpu = MatBindToCPU_HYPRE;
2433: /* Get hypre's default memory location. Users can control this using the corresponding HYPRE_SetMemoryLocation API */
2434: PetscCallHYPRE(HYPRE_GetMemoryLocation(&memory_location));
2435: B->boundtocpu = (memory_location == HYPRE_MEMORY_HOST) ? PETSC_TRUE : PETSC_FALSE;
2436: #endif
2438: /* build cache for off array entries formed */
2439: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
2441: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm));
2442: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE));
2443: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ));
2444: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS));
2445: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE));
2446: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE));
2447: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE));
2448: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE));
2449: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE));
2450: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE));
2451: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2452: #if defined(HYPRE_USING_HIP)
2453: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2454: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2455: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2456: PetscCall(MatSetVecType(B, VECHIP));
2457: #endif
2458: #if defined(HYPRE_USING_CUDA)
2459: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2460: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2461: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2462: PetscCall(MatSetVecType(B, VECCUDA));
2463: #endif
2464: #endif
2465: PetscFunctionReturn(PETSC_SUCCESS);
2466: }