Actual source code: mhypre.c
2: /*
3: Creates hypre ijmatrix from PETSc matrix
4: */
6: #include <petscpkg_version.h>
7: #include <petsc/private/petschypre.h>
8: #include <petscmathypre.h>
9: #include <petsc/private/matimpl.h>
10: #include <petsc/private/deviceimpl.h>
11: #include <../src/mat/impls/hypre/mhypre.h>
12: #include <../src/mat/impls/aij/mpi/mpiaij.h>
13: #include <../src/vec/vec/impls/hypre/vhyp.h>
14: #include <HYPRE.h>
15: #include <HYPRE_utilities.h>
16: #include <_hypre_parcsr_ls.h>
17: #include <_hypre_sstruct_ls.h>
19: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
20: #define hypre_ParCSRMatrixClone(A,B) hypre_ParCSRMatrixCompleteClone(A)
21: #endif
23: static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
24: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
25: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
26: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
27: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,HYPRE_Complex,Vec,HYPRE_Complex,Vec,PetscBool);
28: static PetscErrorCode hypre_array_destroy(void*);
29: static PetscErrorCode MatSetValues_HYPRE(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode ins);
31: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
32: {
33: PetscInt i,n_d,n_o;
34: const PetscInt *ia_d,*ia_o;
35: PetscBool done_d=PETSC_FALSE,done_o=PETSC_FALSE;
36: HYPRE_Int *nnz_d=NULL,*nnz_o=NULL;
38: if (A_d) { /* determine number of nonzero entries in local diagonal part */
39: MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);
40: if (done_d) {
41: PetscMalloc1(n_d,&nnz_d);
42: for (i=0; i<n_d; i++) {
43: nnz_d[i] = ia_d[i+1] - ia_d[i];
44: }
45: }
46: MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);
47: }
48: if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
49: MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
50: if (done_o) {
51: PetscMalloc1(n_o,&nnz_o);
52: for (i=0; i<n_o; i++) {
53: nnz_o[i] = ia_o[i+1] - ia_o[i];
54: }
55: }
56: MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
57: }
58: if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */
59: if (!done_o) { /* only diagonal part */
60: PetscCalloc1(n_d,&nnz_o);
61: }
62: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
63: { /* If we don't do this, the columns of the matrix will be all zeros! */
64: hypre_AuxParCSRMatrix *aux_matrix;
65: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
66: hypre_AuxParCSRMatrixDestroy(aux_matrix);
67: hypre_IJMatrixTranslator(ij) = NULL;
68: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,ij,nnz_d,nnz_o);
69: /* it seems they partially fixed it in 2.19.0 */
70: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
71: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
72: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
73: #endif
74: }
75: #else
76: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,ij,nnz_d,nnz_o);
77: #endif
78: PetscFree(nnz_d);
79: PetscFree(nnz_o);
80: }
81: return 0;
82: }
84: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
85: {
86: PetscInt rstart,rend,cstart,cend;
88: PetscLayoutSetUp(A->rmap);
89: PetscLayoutSetUp(A->cmap);
90: rstart = A->rmap->rstart;
91: rend = A->rmap->rend;
92: cstart = A->cmap->rstart;
93: cend = A->cmap->rend;
94: PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij);
95: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,hA->ij,HYPRE_PARCSR);
96: {
97: PetscBool same;
98: Mat A_d,A_o;
99: const PetscInt *colmap;
100: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&same);
101: if (same) {
102: MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);
103: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
104: return 0;
105: }
106: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&same);
107: if (same) {
108: MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);
109: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
110: return 0;
111: }
112: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&same);
113: if (same) {
114: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
115: return 0;
116: }
117: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQBAIJ,&same);
118: if (same) {
119: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
120: return 0;
121: }
122: }
123: return 0;
124: }
126: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
127: {
128: PetscInt i,rstart,rend,ncols,nr,nc;
129: const PetscScalar *values;
130: const PetscInt *cols;
131: PetscBool flg;
133: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
134: PetscStackCallStandard(HYPRE_IJMatrixInitialize,ij);
135: #else
136: PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,ij,HYPRE_MEMORY_HOST);
137: #endif
138: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
139: MatGetSize(A,&nr,&nc);
140: if (flg && nr == nc) {
141: MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
142: return 0;
143: }
144: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
145: if (flg) {
146: MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
147: return 0;
148: }
150: /* Do not need Aux since we have done precise i[],j[] allocation in MatHYPRE_CreateFromMat() */
151: hypre_AuxParCSRMatrixNeedAux((hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij)) = 0;
153: MatGetOwnershipRange(A,&rstart,&rend);
154: for (i=rstart; i<rend; i++) {
155: MatGetRow(A,i,&ncols,&cols,&values);
156: if (ncols) {
157: HYPRE_Int nc = (HYPRE_Int)ncols;
160: PetscStackCallStandard(HYPRE_IJMatrixSetValues,ij,1,&nc,(HYPRE_BigInt *)&i,(HYPRE_BigInt *)cols,(HYPRE_Complex *)values);
161: }
162: MatRestoreRow(A,i,&ncols,&cols,&values);
163: }
164: return 0;
165: }
167: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
168: {
169: Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data;
170: HYPRE_Int type;
171: hypre_ParCSRMatrix *par_matrix;
172: hypre_AuxParCSRMatrix *aux_matrix;
173: hypre_CSRMatrix *hdiag;
174: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
175: const PetscScalar *pa;
177: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,ij,&type);
179: PetscStackCallStandard(HYPRE_IJMatrixGetObject,ij,(void**)&par_matrix);
180: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
181: /*
182: this is the Hack part where we monkey directly with the hypre datastructures
183: */
184: if (sameint) {
185: PetscArraycpy(hdiag->i,pdiag->i,A->rmap->n + 1);
186: PetscArraycpy(hdiag->j,pdiag->j,pdiag->nz);
187: } else {
188: PetscInt i;
190: for (i=0;i<A->rmap->n + 1;i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
191: for (i=0;i<pdiag->nz;i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
192: }
194: MatSeqAIJGetArrayRead(A,&pa);
195: PetscArraycpy(hdiag->data,pa,pdiag->nz);
196: MatSeqAIJRestoreArrayRead(A,&pa);
198: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
199: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
200: return 0;
201: }
203: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
204: {
205: Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data;
206: Mat_SeqAIJ *pdiag,*poffd;
207: PetscInt i,*garray = pA->garray,*jj,cstart,*pjj;
208: HYPRE_Int *hjj,type;
209: hypre_ParCSRMatrix *par_matrix;
210: hypre_AuxParCSRMatrix *aux_matrix;
211: hypre_CSRMatrix *hdiag,*hoffd;
212: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
213: const PetscScalar *pa;
215: pdiag = (Mat_SeqAIJ*) pA->A->data;
216: poffd = (Mat_SeqAIJ*) pA->B->data;
217: /* cstart is only valid for square MPIAIJ layed out in the usual way */
218: MatGetOwnershipRange(A,&cstart,NULL);
220: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,ij,&type);
222: PetscStackCallStandard(HYPRE_IJMatrixGetObject,ij,(void**)&par_matrix);
223: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
224: hoffd = hypre_ParCSRMatrixOffd(par_matrix);
226: /*
227: this is the Hack part where we monkey directly with the hypre datastructures
228: */
229: if (sameint) {
230: PetscArraycpy(hdiag->i,pdiag->i,pA->A->rmap->n + 1);
231: } else {
232: for (i=0; i<pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)(pdiag->i[i]);
233: }
234: /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
235: hjj = hdiag->j;
236: pjj = pdiag->j;
237: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
238: for (i=0; i<pdiag->nz; i++) hjj[i] = pjj[i];
239: #else
240: for (i=0; i<pdiag->nz; i++) hjj[i] = cstart + pjj[i];
241: #endif
242: MatSeqAIJGetArrayRead(pA->A,&pa);
243: PetscArraycpy(hdiag->data,pa,pdiag->nz);
244: MatSeqAIJRestoreArrayRead(pA->A,&pa);
245: if (sameint) {
246: PetscArraycpy(hoffd->i,poffd->i,pA->A->rmap->n + 1);
247: } else {
248: for (i=0; i<pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)(poffd->i[i]);
249: }
251: /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
252: If we hacked a hypre a bit more we might be able to avoid this step */
253: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
254: PetscStackCallStandard(hypre_CSRMatrixBigInitialize,hoffd);
255: jj = (PetscInt*) hoffd->big_j;
256: #else
257: jj = (PetscInt*) hoffd->j;
258: #endif
259: pjj = poffd->j;
260: for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
262: MatSeqAIJGetArrayRead(pA->B,&pa);
263: PetscArraycpy(hoffd->data,pa,poffd->nz);
264: MatSeqAIJRestoreArrayRead(pA->B,&pa);
266: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
267: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
268: return 0;
269: }
271: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
272: {
273: Mat_HYPRE* mhA = (Mat_HYPRE*)(A->data);
274: Mat lA;
275: ISLocalToGlobalMapping rl2g,cl2g;
276: IS is;
277: hypre_ParCSRMatrix *hA;
278: hypre_CSRMatrix *hdiag,*hoffd;
279: MPI_Comm comm;
280: HYPRE_Complex *hdd,*hod,*aa;
281: PetscScalar *data;
282: HYPRE_BigInt *col_map_offd;
283: HYPRE_Int *hdi,*hdj,*hoi,*hoj;
284: PetscInt *ii,*jj,*iptr,*jptr;
285: PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
286: HYPRE_Int type;
288: comm = PetscObjectComm((PetscObject)A);
289: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,mhA->ij,&type);
291: PetscStackCallStandard(HYPRE_IJMatrixGetObject,mhA->ij,(void**)&hA);
292: M = hypre_ParCSRMatrixGlobalNumRows(hA);
293: N = hypre_ParCSRMatrixGlobalNumCols(hA);
294: str = hypre_ParCSRMatrixFirstRowIndex(hA);
295: stc = hypre_ParCSRMatrixFirstColDiag(hA);
296: hdiag = hypre_ParCSRMatrixDiag(hA);
297: hoffd = hypre_ParCSRMatrixOffd(hA);
298: dr = hypre_CSRMatrixNumRows(hdiag);
299: dc = hypre_CSRMatrixNumCols(hdiag);
300: nnz = hypre_CSRMatrixNumNonzeros(hdiag);
301: hdi = hypre_CSRMatrixI(hdiag);
302: hdj = hypre_CSRMatrixJ(hdiag);
303: hdd = hypre_CSRMatrixData(hdiag);
304: oc = hypre_CSRMatrixNumCols(hoffd);
305: nnz += hypre_CSRMatrixNumNonzeros(hoffd);
306: hoi = hypre_CSRMatrixI(hoffd);
307: hoj = hypre_CSRMatrixJ(hoffd);
308: hod = hypre_CSRMatrixData(hoffd);
309: if (reuse != MAT_REUSE_MATRIX) {
310: PetscInt *aux;
312: /* generate l2g maps for rows and cols */
313: ISCreateStride(comm,dr,str,1,&is);
314: ISLocalToGlobalMappingCreateIS(is,&rl2g);
315: ISDestroy(&is);
316: col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
317: PetscMalloc1(dc+oc,&aux);
318: for (i=0; i<dc; i++) aux[i] = i+stc;
319: for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
320: ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
321: ISLocalToGlobalMappingCreateIS(is,&cl2g);
322: ISDestroy(&is);
323: /* create MATIS object */
324: MatCreate(comm,B);
325: MatSetSizes(*B,dr,dc,M,N);
326: MatSetType(*B,MATIS);
327: MatSetLocalToGlobalMapping(*B,rl2g,cl2g);
328: ISLocalToGlobalMappingDestroy(&rl2g);
329: ISLocalToGlobalMappingDestroy(&cl2g);
331: /* allocate CSR for local matrix */
332: PetscMalloc1(dr+1,&iptr);
333: PetscMalloc1(nnz,&jptr);
334: PetscMalloc1(nnz,&data);
335: } else {
336: PetscInt nr;
337: PetscBool done;
338: MatISGetLocalMat(*B,&lA);
339: MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);
342: MatSeqAIJGetArray(lA,&data);
343: }
344: /* merge local matrices */
345: ii = iptr;
346: jj = jptr;
347: 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++ */
348: *ii = *(hdi++) + *(hoi++);
349: for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
350: PetscScalar *aold = (PetscScalar*)aa;
351: PetscInt *jold = jj,nc = jd+jo;
352: for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; }
353: for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
354: *(++ii) = *(hdi++) + *(hoi++);
355: PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);
356: }
357: for (; cum<dr; cum++) *(++ii) = nnz;
358: if (reuse != MAT_REUSE_MATRIX) {
359: Mat_SeqAIJ* a;
361: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);
362: MatISSetLocalMat(*B,lA);
363: /* hack SeqAIJ */
364: a = (Mat_SeqAIJ*)(lA->data);
365: a->free_a = PETSC_TRUE;
366: a->free_ij = PETSC_TRUE;
367: MatDestroy(&lA);
368: }
369: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
370: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
371: if (reuse == MAT_INPLACE_MATRIX) {
372: MatHeaderReplace(A,B);
373: }
374: return 0;
375: }
377: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
378: {
379: Mat M = NULL;
380: Mat_HYPRE *hB;
381: MPI_Comm comm = PetscObjectComm((PetscObject)A);
383: if (reuse == MAT_REUSE_MATRIX) {
384: /* always destroy the old matrix and create a new memory;
385: hope this does not churn the memory too much. The problem
386: is I do not know if it is possible to put the matrix back to
387: its initial state so that we can directly copy the values
388: the second time through. */
389: hB = (Mat_HYPRE*)((*B)->data);
390: PetscStackCallStandard(HYPRE_IJMatrixDestroy,hB->ij);
391: } else {
392: MatCreate(comm,&M);
393: MatSetType(M,MATHYPRE);
394: MatSetSizes(M,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
395: hB = (Mat_HYPRE*)(M->data);
396: if (reuse == MAT_INITIAL_MATRIX) *B = M;
397: }
398: MatSetOption(*B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
399: MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
400: MatHYPRE_CreateFromMat(A,hB);
401: MatHYPRE_IJMatrixCopy(A,hB->ij);
402: if (reuse == MAT_INPLACE_MATRIX) {
403: MatHeaderReplace(A,&M);
404: }
405: (*B)->preallocated = PETSC_TRUE;
406: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
407: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
408: return 0;
409: }
411: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
412: {
413: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
414: hypre_ParCSRMatrix *parcsr;
415: hypre_CSRMatrix *hdiag,*hoffd;
416: MPI_Comm comm;
417: PetscScalar *da,*oa,*aptr;
418: PetscInt *dii,*djj,*oii,*ojj,*iptr;
419: PetscInt i,dnnz,onnz,m,n;
420: HYPRE_Int type;
421: PetscMPIInt size;
422: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
424: comm = PetscObjectComm((PetscObject)A);
425: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type);
427: if (reuse == MAT_REUSE_MATRIX) {
428: PetscBool ismpiaij,isseqaij;
429: PetscObjectBaseTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
430: PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
432: }
433: #if defined(PETSC_HAVE_HYPRE_DEVICE)
435: #endif
436: MPI_Comm_size(comm,&size);
438: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
439: hdiag = hypre_ParCSRMatrixDiag(parcsr);
440: hoffd = hypre_ParCSRMatrixOffd(parcsr);
441: m = hypre_CSRMatrixNumRows(hdiag);
442: n = hypre_CSRMatrixNumCols(hdiag);
443: dnnz = hypre_CSRMatrixNumNonzeros(hdiag);
444: onnz = hypre_CSRMatrixNumNonzeros(hoffd);
445: if (reuse == MAT_INITIAL_MATRIX) {
446: PetscMalloc1(m+1,&dii);
447: PetscMalloc1(dnnz,&djj);
448: PetscMalloc1(dnnz,&da);
449: } else if (reuse == MAT_REUSE_MATRIX) {
450: PetscInt nr;
451: PetscBool done;
452: if (size > 1) {
453: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
455: MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
458: MatSeqAIJGetArray(b->A,&da);
459: } else {
460: MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
463: MatSeqAIJGetArray(*B,&da);
464: }
465: } else { /* MAT_INPLACE_MATRIX */
466: if (!sameint) {
467: PetscMalloc1(m+1,&dii);
468: PetscMalloc1(dnnz,&djj);
469: } else {
470: dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
471: djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
472: }
473: da = (PetscScalar*)hypre_CSRMatrixData(hdiag);
474: }
476: if (!sameint) {
477: if (reuse != MAT_REUSE_MATRIX) { for (i=0;i<m+1;i++) dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]); }
478: for (i=0;i<dnnz;i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]);
479: } else {
480: if (reuse != MAT_REUSE_MATRIX) PetscArraycpy(dii,hypre_CSRMatrixI(hdiag),m+1);
481: PetscArraycpy(djj,hypre_CSRMatrixJ(hdiag),dnnz);
482: }
483: PetscArraycpy(da,hypre_CSRMatrixData(hdiag),dnnz);
484: iptr = djj;
485: aptr = da;
486: for (i=0; i<m; i++) {
487: PetscInt nc = dii[i+1]-dii[i];
488: PetscSortIntWithScalarArray(nc,iptr,aptr);
489: iptr += nc;
490: aptr += nc;
491: }
492: if (size > 1) {
493: HYPRE_BigInt *coffd;
494: HYPRE_Int *offdj;
496: if (reuse == MAT_INITIAL_MATRIX) {
497: PetscMalloc1(m+1,&oii);
498: PetscMalloc1(onnz,&ojj);
499: PetscMalloc1(onnz,&oa);
500: } else if (reuse == MAT_REUSE_MATRIX) {
501: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
502: PetscInt nr,hr = hypre_CSRMatrixNumRows(hoffd);
503: PetscBool done;
505: MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
508: MatSeqAIJGetArray(b->B,&oa);
509: } else { /* MAT_INPLACE_MATRIX */
510: if (!sameint) {
511: PetscMalloc1(m+1,&oii);
512: PetscMalloc1(onnz,&ojj);
513: } else {
514: oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
515: ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
516: }
517: oa = (PetscScalar*)hypre_CSRMatrixData(hoffd);
518: }
519: if (reuse != MAT_REUSE_MATRIX) {
520: if (!sameint) {
521: for (i=0;i<m+1;i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]);
522: } else {
523: PetscArraycpy(oii,hypre_CSRMatrixI(hoffd),m+1);
524: }
525: }
526: PetscArraycpy(oa,hypre_CSRMatrixData(hoffd),onnz);
528: offdj = hypre_CSRMatrixJ(hoffd);
529: coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
530: /* we only need the permutation to be computed properly, I don't know if HYPRE
531: messes up with the ordering. Just in case, allocate some memory and free it
532: later */
533: if (reuse == MAT_REUSE_MATRIX) {
534: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
535: PetscInt mnz;
537: MatSeqAIJGetMaxRowNonzeros(b->B,&mnz);
538: PetscMalloc1(mnz,&ojj);
539: } else for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
540: iptr = ojj;
541: aptr = oa;
542: for (i=0; i<m; i++) {
543: PetscInt nc = oii[i+1]-oii[i];
544: if (reuse == MAT_REUSE_MATRIX) {
545: PetscInt j;
547: iptr = ojj;
548: for (j=0; j<nc; j++) iptr[j] = coffd[offdj[oii[i] + j]];
549: }
550: PetscSortIntWithScalarArray(nc,iptr,aptr);
551: iptr += nc;
552: aptr += nc;
553: }
554: if (reuse == MAT_REUSE_MATRIX) PetscFree(ojj);
555: if (reuse == MAT_INITIAL_MATRIX) {
556: Mat_MPIAIJ *b;
557: Mat_SeqAIJ *d,*o;
559: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
560: /* hack MPIAIJ */
561: b = (Mat_MPIAIJ*)((*B)->data);
562: d = (Mat_SeqAIJ*)b->A->data;
563: o = (Mat_SeqAIJ*)b->B->data;
564: d->free_a = PETSC_TRUE;
565: d->free_ij = PETSC_TRUE;
566: o->free_a = PETSC_TRUE;
567: o->free_ij = PETSC_TRUE;
568: } else if (reuse == MAT_INPLACE_MATRIX) {
569: Mat T;
571: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
572: if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
573: hypre_CSRMatrixI(hdiag) = NULL;
574: hypre_CSRMatrixJ(hdiag) = NULL;
575: hypre_CSRMatrixI(hoffd) = NULL;
576: hypre_CSRMatrixJ(hoffd) = NULL;
577: } else { /* Hack MPIAIJ -> free ij but not a */
578: Mat_MPIAIJ *b = (Mat_MPIAIJ*)(T->data);
579: Mat_SeqAIJ *d = (Mat_SeqAIJ*)(b->A->data);
580: Mat_SeqAIJ *o = (Mat_SeqAIJ*)(b->B->data);
582: d->free_ij = PETSC_TRUE;
583: o->free_ij = PETSC_TRUE;
584: }
585: hypre_CSRMatrixData(hdiag) = NULL;
586: hypre_CSRMatrixData(hoffd) = NULL;
587: MatHeaderReplace(A,&T);
588: }
589: } else {
590: oii = NULL;
591: ojj = NULL;
592: oa = NULL;
593: if (reuse == MAT_INITIAL_MATRIX) {
594: Mat_SeqAIJ* b;
596: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
597: /* hack SeqAIJ */
598: b = (Mat_SeqAIJ*)((*B)->data);
599: b->free_a = PETSC_TRUE;
600: b->free_ij = PETSC_TRUE;
601: } else if (reuse == MAT_INPLACE_MATRIX) {
602: Mat T;
604: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
605: if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
606: hypre_CSRMatrixI(hdiag) = NULL;
607: hypre_CSRMatrixJ(hdiag) = NULL;
608: } else { /* free ij but not a */
609: Mat_SeqAIJ* b = (Mat_SeqAIJ*)(T->data);
611: b->free_ij = PETSC_TRUE;
612: }
613: hypre_CSRMatrixData(hdiag) = NULL;
614: MatHeaderReplace(A,&T);
615: }
616: }
618: /* we have to use hypre_Tfree to free the HYPRE arrays
619: that PETSc now onws */
620: if (reuse == MAT_INPLACE_MATRIX) {
621: PetscInt nh;
622: void *ptrs[6] = {da,oa,dii,djj,oii,ojj};
623: const char *names[6] = {"_hypre_csr_da",
624: "_hypre_csr_oa",
625: "_hypre_csr_dii",
626: "_hypre_csr_djj",
627: "_hypre_csr_oii",
628: "_hypre_csr_ojj"};
629: nh = sameint ? 6 : 2;
630: for (i=0; i<nh; i++) {
631: PetscContainer c;
633: PetscContainerCreate(comm,&c);
634: PetscContainerSetPointer(c,ptrs[i]);
635: PetscContainerSetUserDestroy(c,hypre_array_destroy);
636: PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
637: PetscContainerDestroy(&c);
638: }
639: }
640: return 0;
641: }
643: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
644: {
645: hypre_ParCSRMatrix *tA;
646: hypre_CSRMatrix *hdiag,*hoffd;
647: Mat_SeqAIJ *diag,*offd;
648: PetscInt *garray,i,noffd,dnnz,onnz,*row_starts,*col_starts;
649: MPI_Comm comm = PetscObjectComm((PetscObject)A);
650: PetscBool ismpiaij,isseqaij;
651: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
652: HYPRE_Int *hdi = NULL,*hdj = NULL,*hoi = NULL,*hoj = NULL;
653: PetscInt *pdi = NULL,*pdj = NULL,*poi = NULL,*poj = NULL;
654: #if defined(PETSC_HAVE_HYPRE_DEVICE)
655: PetscBool iscuda = PETSC_FALSE;
656: #endif
658: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
659: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
661: if (ismpiaij) {
662: Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);
664: diag = (Mat_SeqAIJ*)a->A->data;
665: offd = (Mat_SeqAIJ*)a->B->data;
666: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE) && defined(HYPRE_USING_CUDA)
667: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJCUSPARSE,&iscuda);
668: if (iscuda && !A->boundtocpu) {
669: sameint = PETSC_TRUE;
670: MatSeqAIJCUSPARSEGetIJ(a->A,PETSC_FALSE,(const HYPRE_Int**)&hdi,(const HYPRE_Int**)&hdj);
671: MatSeqAIJCUSPARSEGetIJ(a->B,PETSC_FALSE,(const HYPRE_Int**)&hoi,(const HYPRE_Int**)&hoj);
672: } else {
673: #else
674: {
675: #endif
676: pdi = diag->i;
677: pdj = diag->j;
678: poi = offd->i;
679: poj = offd->j;
680: if (sameint) {
681: hdi = (HYPRE_Int*)pdi;
682: hdj = (HYPRE_Int*)pdj;
683: hoi = (HYPRE_Int*)poi;
684: hoj = (HYPRE_Int*)poj;
685: }
686: }
687: garray = a->garray;
688: noffd = a->B->cmap->N;
689: dnnz = diag->nz;
690: onnz = offd->nz;
691: } else {
692: diag = (Mat_SeqAIJ*)A->data;
693: offd = NULL;
694: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE)
695: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJCUSPARSE,&iscuda);
696: if (iscuda && !A->boundtocpu) {
697: sameint = PETSC_TRUE;
698: MatSeqAIJCUSPARSEGetIJ(A,PETSC_FALSE,(const HYPRE_Int**)&hdi,(const HYPRE_Int**)&hdj);
699: } else {
700: #else
701: {
702: #endif
703: pdi = diag->i;
704: pdj = diag->j;
705: if (sameint) {
706: hdi = (HYPRE_Int*)pdi;
707: hdj = (HYPRE_Int*)pdj;
708: }
709: }
710: garray = NULL;
711: noffd = 0;
712: dnnz = diag->nz;
713: onnz = 0;
714: }
716: /* create a temporary ParCSR */
717: if (HYPRE_AssumedPartitionCheck()) {
718: PetscMPIInt myid;
720: MPI_Comm_rank(comm,&myid);
721: row_starts = A->rmap->range + myid;
722: col_starts = A->cmap->range + myid;
723: } else {
724: row_starts = A->rmap->range;
725: col_starts = A->cmap->range;
726: }
727: tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_BigInt*)row_starts,(HYPRE_BigInt*)col_starts,noffd,dnnz,onnz);
728: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
729: hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
730: hypre_ParCSRMatrixSetColStartsOwner(tA,0);
731: #endif
733: /* set diagonal part */
734: hdiag = hypre_ParCSRMatrixDiag(tA);
735: if (!sameint) { /* malloc CSR pointers */
736: PetscMalloc2(A->rmap->n+1,&hdi,dnnz,&hdj);
737: for (i = 0; i < A->rmap->n+1; i++) hdi[i] = (HYPRE_Int)(pdi[i]);
738: for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)(pdj[i]);
739: }
740: hypre_CSRMatrixI(hdiag) = hdi;
741: hypre_CSRMatrixJ(hdiag) = hdj;
742: hypre_CSRMatrixData(hdiag) = (HYPRE_Complex*)diag->a;
743: hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
744: hypre_CSRMatrixSetRownnz(hdiag);
745: hypre_CSRMatrixSetDataOwner(hdiag,0);
747: /* set offdiagonal part */
748: hoffd = hypre_ParCSRMatrixOffd(tA);
749: if (offd) {
750: if (!sameint) { /* malloc CSR pointers */
751: PetscMalloc2(A->rmap->n+1,&hoi,onnz,&hoj);
752: for (i = 0; i < A->rmap->n+1; i++) hoi[i] = (HYPRE_Int)(poi[i]);
753: for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)(poj[i]);
754: }
755: hypre_CSRMatrixI(hoffd) = hoi;
756: hypre_CSRMatrixJ(hoffd) = hoj;
757: hypre_CSRMatrixData(hoffd) = (HYPRE_Complex*)offd->a;
758: hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
759: hypre_CSRMatrixSetRownnz(hoffd);
760: hypre_CSRMatrixSetDataOwner(hoffd,0);
761: }
762: #if defined(PETSC_HAVE_HYPRE_DEVICE)
763: PetscStackCallStandard(hypre_ParCSRMatrixInitialize_v2,tA,iscuda ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST);
764: #else
765: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
766: PetscStackCallStandard(hypre_ParCSRMatrixInitialize,tA);
767: #else
768: PetscStackCallStandard(hypre_ParCSRMatrixInitialize_v2,tA,HYPRE_MEMORY_HOST);
769: #endif
770: #endif
771: hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA),HYPRE_MEMORY_HOST);
772: hypre_ParCSRMatrixSetNumNonzeros(tA);
773: hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt*)garray;
774: if (!hypre_ParCSRMatrixCommPkg(tA)) PetscStackCallStandard(hypre_MatvecCommPkgCreate,tA);
775: *hA = tA;
776: return 0;
777: }
779: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
780: {
781: hypre_CSRMatrix *hdiag,*hoffd;
782: PetscBool ismpiaij,sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
783: #if defined(PETSC_HAVE_HYPRE_DEVICE)
784: PetscBool iscuda = PETSC_FALSE;
785: #endif
787: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
788: #if defined(PETSC_HAVE_HYPRE_DEVICE)
789: PetscObjectTypeCompareAny((PetscObject)A,&iscuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
790: if (iscuda) sameint = PETSC_TRUE;
791: #endif
792: hdiag = hypre_ParCSRMatrixDiag(*hA);
793: hoffd = hypre_ParCSRMatrixOffd(*hA);
794: /* free temporary memory allocated by PETSc
795: set pointers to NULL before destroying tA */
796: if (!sameint) {
797: HYPRE_Int *hi,*hj;
799: hi = hypre_CSRMatrixI(hdiag);
800: hj = hypre_CSRMatrixJ(hdiag);
801: PetscFree2(hi,hj);
802: if (ismpiaij) {
803: hi = hypre_CSRMatrixI(hoffd);
804: hj = hypre_CSRMatrixJ(hoffd);
805: PetscFree2(hi,hj);
806: }
807: }
808: hypre_CSRMatrixI(hdiag) = NULL;
809: hypre_CSRMatrixJ(hdiag) = NULL;
810: hypre_CSRMatrixData(hdiag) = NULL;
811: if (ismpiaij) {
812: hypre_CSRMatrixI(hoffd) = NULL;
813: hypre_CSRMatrixJ(hoffd) = NULL;
814: hypre_CSRMatrixData(hoffd) = NULL;
815: }
816: hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
817: hypre_ParCSRMatrixDestroy(*hA);
818: *hA = NULL;
819: return 0;
820: }
822: /* calls RAP from BoomerAMG:
823: the resulting ParCSR will not own the column and row starts
824: It looks like we don't need to have the diagonal entries ordered first */
825: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
826: {
827: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
828: HYPRE_Int P_owns_col_starts,R_owns_row_starts;
829: #endif
831: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
832: P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
833: R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
834: #endif
835: /* can be replaced by version test later */
836: #if defined(PETSC_HAVE_HYPRE_DEVICE)
837: PetscStackPush("hypre_ParCSRMatrixRAP");
838: *hRAP = hypre_ParCSRMatrixRAP(hR,hA,hP);
839: PetscStackPop;
840: #else
841: PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,hR,hA,hP,hRAP);
842: PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,*hRAP);
843: #endif
844: /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
845: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
846: hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
847: hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
848: if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
849: if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
850: #endif
851: return 0;
852: }
854: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
855: {
856: Mat B;
857: hypre_ParCSRMatrix *hA,*hP,*hPtAP = NULL;
858: Mat_Product *product=C->product;
860: MatAIJGetParCSR_Private(A,&hA);
861: MatAIJGetParCSR_Private(P,&hP);
862: MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);
863: MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);
865: MatHeaderMerge(C,&B);
866: C->product = product;
868: MatAIJRestoreParCSR_Private(A,&hA);
869: MatAIJRestoreParCSR_Private(P,&hP);
870: return 0;
871: }
873: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat C)
874: {
875: MatSetType(C,MATAIJ);
876: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
877: C->ops->productnumeric = MatProductNumeric_PtAP;
878: return 0;
879: }
881: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
882: {
883: Mat B;
884: Mat_HYPRE *hP;
885: hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr = NULL;
886: HYPRE_Int type;
887: MPI_Comm comm = PetscObjectComm((PetscObject)A);
888: PetscBool ishypre;
890: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
892: hP = (Mat_HYPRE*)P->data;
893: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hP->ij,&type);
895: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hP->ij,(void**)&Pparcsr);
897: MatAIJGetParCSR_Private(A,&hA);
898: MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
899: MatAIJRestoreParCSR_Private(A,&hA);
901: /* create temporary matrix and merge to C */
902: MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
903: MatHeaderMerge(C,&B);
904: return 0;
905: }
907: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
908: {
909: Mat B;
910: hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr = NULL;
911: Mat_HYPRE *hA,*hP;
912: PetscBool ishypre;
913: HYPRE_Int type;
915: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
917: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
919: hA = (Mat_HYPRE*)A->data;
920: hP = (Mat_HYPRE*)P->data;
921: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type);
923: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hP->ij,&type);
925: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&Aparcsr);
926: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hP->ij,(void**)&Pparcsr);
927: MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
928: MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
929: MatHeaderMerge(C,&B);
930: return 0;
931: }
933: /* calls hypre_ParMatmul
934: hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
935: hypre_ParMatrixCreate does not duplicate the communicator
936: It looks like we don't need to have the diagonal entries ordered first */
937: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
938: {
939: /* can be replaced by version test later */
940: #if defined(PETSC_HAVE_HYPRE_DEVICE)
941: PetscStackPush("hypre_ParCSRMatMat");
942: *hAB = hypre_ParCSRMatMat(hA,hB);
943: #else
944: PetscStackPush("hypre_ParMatmul");
945: *hAB = hypre_ParMatmul(hA,hB);
946: #endif
947: PetscStackPop;
948: return 0;
949: }
951: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
952: {
953: Mat D;
954: hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
955: Mat_Product *product=C->product;
957: MatAIJGetParCSR_Private(A,&hA);
958: MatAIJGetParCSR_Private(B,&hB);
959: MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
960: MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);
962: MatHeaderMerge(C,&D);
963: C->product = product;
965: MatAIJRestoreParCSR_Private(A,&hA);
966: MatAIJRestoreParCSR_Private(B,&hB);
967: return 0;
968: }
970: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat C)
971: {
972: MatSetType(C,MATAIJ);
973: C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
974: C->ops->productnumeric = MatProductNumeric_AB;
975: return 0;
976: }
978: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
979: {
980: Mat D;
981: hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
982: Mat_HYPRE *hA,*hB;
983: PetscBool ishypre;
984: HYPRE_Int type;
985: Mat_Product *product;
987: PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
989: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
991: hA = (Mat_HYPRE*)A->data;
992: hB = (Mat_HYPRE*)B->data;
993: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type);
995: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hB->ij,&type);
997: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&Aparcsr);
998: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hB->ij,(void**)&Bparcsr);
999: MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
1000: MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);
1002: /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1003: product = C->product; /* save it from MatHeaderReplace() */
1004: C->product = NULL;
1005: MatHeaderReplace(C,&D);
1006: C->product = product;
1007: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1008: C->ops->productnumeric = MatProductNumeric_AB;
1009: return 0;
1010: }
1012: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
1013: {
1014: Mat E;
1015: hypre_ParCSRMatrix *hA,*hB,*hC,*hABC = NULL;
1017: MatAIJGetParCSR_Private(A,&hA);
1018: MatAIJGetParCSR_Private(B,&hB);
1019: MatAIJGetParCSR_Private(C,&hC);
1020: MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
1021: MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
1022: MatHeaderMerge(D,&E);
1023: MatAIJRestoreParCSR_Private(A,&hA);
1024: MatAIJRestoreParCSR_Private(B,&hB);
1025: MatAIJRestoreParCSR_Private(C,&hC);
1026: return 0;
1027: }
1029: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
1030: {
1031: MatSetType(D,MATAIJ);
1032: return 0;
1033: }
1035: /* ---------------------------------------------------- */
1036: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1037: {
1038: C->ops->productnumeric = MatProductNumeric_AB;
1039: return 0;
1040: }
1042: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1043: {
1044: Mat_Product *product = C->product;
1045: PetscBool Ahypre;
1047: PetscObjectTypeCompare((PetscObject)product->A,MATHYPRE,&Ahypre);
1048: if (Ahypre) { /* A is a Hypre matrix */
1049: MatSetType(C,MATHYPRE);
1050: C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1051: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1052: return 0;
1053: }
1054: return 0;
1055: }
1057: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1058: {
1059: C->ops->productnumeric = MatProductNumeric_PtAP;
1060: return 0;
1061: }
1063: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1064: {
1065: Mat_Product *product = C->product;
1066: PetscBool flg;
1067: PetscInt type = 0;
1068: const char *outTypes[4] = {"aij","seqaij","mpiaij","hypre"};
1069: PetscInt ntype = 4;
1070: Mat A = product->A;
1071: PetscBool Ahypre;
1074: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&Ahypre);
1075: if (Ahypre) { /* A is a Hypre matrix */
1076: MatSetType(C,MATHYPRE);
1077: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1078: C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
1079: return 0;
1080: }
1082: /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1083: /* Get runtime option */
1084: if (product->api_user) {
1085: PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatPtAP_HYPRE","Mat");
1086: PetscOptionsEList("-matptap_hypre_outtype","MatPtAP outtype","MatPtAP outtype",outTypes,ntype,outTypes[type],&type,&flg);
1087: PetscOptionsEnd();
1088: } else {
1089: PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_PtAP_HYPRE","Mat");
1090: PetscOptionsEList("-mat_product_algorithm_hypre_outtype","MatProduct_PtAP outtype","MatProduct_PtAP",outTypes,ntype,outTypes[type],&type,&flg);
1091: PetscOptionsEnd();
1092: }
1094: if (type == 0 || type == 1 || type == 2) {
1095: MatSetType(C,MATAIJ);
1096: } else if (type == 3) {
1097: MatSetType(C,MATHYPRE);
1098: } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatPtAP outtype is not supported");
1099: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1100: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
1101: return 0;
1102: }
1104: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1105: {
1106: Mat_Product *product = C->product;
1108: switch (product->type) {
1109: case MATPRODUCT_AB:
1110: MatProductSetFromOptions_HYPRE_AB(C);
1111: break;
1112: case MATPRODUCT_PtAP:
1113: MatProductSetFromOptions_HYPRE_PtAP(C);
1114: break;
1115: default:
1116: break;
1117: }
1118: return 0;
1119: }
1121: /* -------------------------------------------------------- */
1123: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1124: {
1125: MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);
1126: return 0;
1127: }
1129: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1130: {
1131: MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);
1132: return 0;
1133: }
1135: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1136: {
1137: if (y != z) {
1138: VecCopy(y,z);
1139: }
1140: MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);
1141: return 0;
1142: }
1144: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1145: {
1146: if (y != z) {
1147: VecCopy(y,z);
1148: }
1149: MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);
1150: return 0;
1151: }
1153: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1154: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1155: {
1156: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1157: hypre_ParCSRMatrix *parcsr;
1158: hypre_ParVector *hx,*hy;
1160: if (trans) {
1161: VecHYPRE_IJVectorPushVecRead(hA->b,x);
1162: if (b != 0.0) VecHYPRE_IJVectorPushVec(hA->x,y);
1163: else VecHYPRE_IJVectorPushVecWrite(hA->x,y);
1164: PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->b->ij,(void**)&hx);
1165: PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->x->ij,(void**)&hy);
1166: } else {
1167: VecHYPRE_IJVectorPushVecRead(hA->x,x);
1168: if (b != 0.0) VecHYPRE_IJVectorPushVec(hA->b,y);
1169: else VecHYPRE_IJVectorPushVecWrite(hA->b,y);
1170: PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->x->ij,(void**)&hx);
1171: PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->b->ij,(void**)&hy);
1172: }
1173: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1174: if (trans) {
1175: PetscStackCallStandard(hypre_ParCSRMatrixMatvecT,a,parcsr,hx,b,hy);
1176: } else {
1177: PetscStackCallStandard(hypre_ParCSRMatrixMatvec,a,parcsr,hx,b,hy);
1178: }
1179: VecHYPRE_IJVectorPopVec(hA->x);
1180: VecHYPRE_IJVectorPopVec(hA->b);
1181: return 0;
1182: }
1184: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1185: {
1186: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1188: VecHYPRE_IJVectorDestroy(&hA->x);
1189: VecHYPRE_IJVectorDestroy(&hA->b);
1190: if (hA->ij) {
1191: if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1192: PetscStackCallStandard(HYPRE_IJMatrixDestroy,hA->ij);
1193: }
1194: if (hA->comm) PetscCommRestoreComm(PetscObjectComm((PetscObject)A),&hA->comm);
1196: MatStashDestroy_Private(&A->stash);
1197: PetscFree(hA->array);
1199: if (hA->cooMat) {
1200: MatDestroy(&hA->cooMat);
1201: PetscStackCall("hypre_TFree",hypre_TFree(hA->diagJ,hA->memType));
1202: PetscStackCall("hypre_TFree",hypre_TFree(hA->offdJ,hA->memType));
1203: PetscStackCall("hypre_TFree",hypre_TFree(hA->diag,hA->memType));
1204: }
1206: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
1207: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
1208: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_hypre_C",NULL);
1209: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_mpiaij_hypre_C",NULL);
1210: PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
1211: PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
1212: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL);
1213: PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL);
1214: PetscFree(A->data);
1215: return 0;
1216: }
1218: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1219: {
1220: MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
1221: return 0;
1222: }
1224: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1225: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1226: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1227: {
1228: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1229: HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
1231: A->boundtocpu = bind;
1232: if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1233: hypre_ParCSRMatrix *parcsr;
1234: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1235: PetscStackCallStandard(hypre_ParCSRMatrixMigrate,parcsr, hmem);
1236: }
1237: if (hA->x) VecHYPRE_IJBindToCPU(hA->x,bind);
1238: if (hA->b) VecHYPRE_IJBindToCPU(hA->b,bind);
1239: return 0;
1240: }
1241: #endif
1243: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1244: {
1245: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1246: PetscMPIInt n;
1247: PetscInt i,j,rstart,ncols,flg;
1248: PetscInt *row,*col;
1249: PetscScalar *val;
1253: if (!A->nooffprocentries) {
1254: while (1) {
1255: MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
1256: if (!flg) break;
1258: for (i=0; i<n;) {
1259: /* Now identify the consecutive vals belonging to the same row */
1260: for (j=i,rstart=row[j]; j<n; j++) {
1261: if (row[j] != rstart) break;
1262: }
1263: if (j < n) ncols = j-i;
1264: else ncols = n-i;
1265: /* Now assemble all these values with a single function call */
1266: MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);
1268: i = j;
1269: }
1270: }
1271: MatStashScatterEnd_Private(&A->stash);
1272: }
1274: PetscStackCallStandard(HYPRE_IJMatrixAssemble,hA->ij);
1275: /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1276: /* 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 */
1277: if (!hA->sorted_full) {
1278: hypre_AuxParCSRMatrix *aux_matrix;
1280: /* call destroy just to make sure we do not leak anything */
1281: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1282: PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,aux_matrix);
1283: hypre_IJMatrixTranslator(hA->ij) = NULL;
1285: /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1286: PetscStackCallStandard(HYPRE_IJMatrixInitialize,hA->ij);
1287: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1288: if (aux_matrix) {
1289: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1290: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1291: PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,aux_matrix);
1292: #else
1293: PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize_v2,aux_matrix,HYPRE_MEMORY_HOST);
1294: #endif
1295: }
1296: }
1297: {
1298: hypre_ParCSRMatrix *parcsr;
1300: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1301: if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscStackCallStandard(hypre_MatvecCommPkgCreate,parcsr);
1302: }
1303: if (!hA->x) VecHYPRE_IJVectorCreate(A->cmap,&hA->x);
1304: if (!hA->b) VecHYPRE_IJVectorCreate(A->rmap,&hA->b);
1305: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1306: MatBindToCPU_HYPRE(A,A->boundtocpu);
1307: #endif
1308: return 0;
1309: }
1311: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1312: {
1313: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1317: if (hA->size >= size) {
1318: *array = hA->array;
1319: } else {
1320: PetscFree(hA->array);
1321: hA->size = size;
1322: PetscMalloc(hA->size,&hA->array);
1323: *array = hA->array;
1324: }
1326: hA->available = PETSC_FALSE;
1327: return 0;
1328: }
1330: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1331: {
1332: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1334: *array = NULL;
1335: hA->available = PETSC_TRUE;
1336: return 0;
1337: }
1339: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1340: {
1341: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1342: PetscScalar *vals = (PetscScalar *)v;
1343: HYPRE_Complex *sscr;
1344: PetscInt *cscr[2];
1345: PetscInt i,nzc;
1346: void *array = NULL;
1348: MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(HYPRE_Complex)*nc*nr,&array);
1349: cscr[0] = (PetscInt*)array;
1350: cscr[1] = ((PetscInt*)array)+nc;
1351: sscr = (HYPRE_Complex*)(((PetscInt*)array)+nc*2);
1352: for (i=0,nzc=0;i<nc;i++) {
1353: if (cols[i] >= 0) {
1354: cscr[0][nzc ] = cols[i];
1355: cscr[1][nzc++] = i;
1356: }
1357: }
1358: if (!nzc) {
1359: MatRestoreArray_HYPRE(A,&array);
1360: return 0;
1361: }
1363: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1364: if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1365: hypre_ParCSRMatrix *parcsr;
1367: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1368: PetscStackCallStandard(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST);
1369: }
1370: #endif
1372: if (ins == ADD_VALUES) {
1373: for (i=0;i<nr;i++) {
1374: if (rows[i] >= 0) {
1375: PetscInt j;
1376: HYPRE_Int hnc = (HYPRE_Int)nzc;
1379: for (j=0;j<nzc;j++) PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]);
1380: PetscStackCallStandard(HYPRE_IJMatrixAddToValues,hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr);
1381: }
1382: vals += nc;
1383: }
1384: } else { /* INSERT_VALUES */
1385: PetscInt rst,ren;
1387: MatGetOwnershipRange(A,&rst,&ren);
1388: for (i=0;i<nr;i++) {
1389: if (rows[i] >= 0) {
1390: PetscInt j;
1391: HYPRE_Int hnc = (HYPRE_Int)nzc;
1394: for (j=0;j<nzc;j++) PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]);
1395: /* nonlocal values */
1396: if (rows[i] < rst || rows[i] >= ren) MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],(PetscScalar*)sscr,PETSC_FALSE);
1397: /* local values */
1398: else PetscStackCallStandard(HYPRE_IJMatrixSetValues,hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr);
1399: }
1400: vals += nc;
1401: }
1402: }
1404: MatRestoreArray_HYPRE(A,&array);
1405: return 0;
1406: }
1408: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1409: {
1410: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1411: HYPRE_Int *hdnnz,*honnz;
1412: PetscInt i,rs,re,cs,ce,bs;
1413: PetscMPIInt size;
1415: PetscLayoutSetUp(A->rmap);
1416: PetscLayoutSetUp(A->cmap);
1417: rs = A->rmap->rstart;
1418: re = A->rmap->rend;
1419: cs = A->cmap->rstart;
1420: ce = A->cmap->rend;
1421: if (!hA->ij) {
1422: PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rs,re-1,cs,ce-1,&hA->ij);
1423: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,hA->ij,HYPRE_PARCSR);
1424: } else {
1425: HYPRE_BigInt hrs,hre,hcs,hce;
1426: PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,hA->ij,&hrs,&hre,&hcs,&hce);
1429: }
1430: MatGetBlockSize(A,&bs);
1431: if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1432: if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;
1434: if (!dnnz) {
1435: PetscMalloc1(A->rmap->n,&hdnnz);
1436: for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1437: } else {
1438: hdnnz = (HYPRE_Int*)dnnz;
1439: }
1440: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1441: if (size > 1) {
1442: hypre_AuxParCSRMatrix *aux_matrix;
1443: if (!onnz) {
1444: PetscMalloc1(A->rmap->n,&honnz);
1445: for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1446: } else honnz = (HYPRE_Int*)onnz;
1447: /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1448: they assume the user will input the entire row values, properly sorted
1449: In PETSc, we don't make such an assumption and set this flag to 1,
1450: unless the option MAT_SORTED_FULL is set to true.
1451: Also, to avoid possible memory leaks, we destroy and recreate the translator
1452: This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1453: the IJ matrix for us */
1454: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1455: hypre_AuxParCSRMatrixDestroy(aux_matrix);
1456: hypre_IJMatrixTranslator(hA->ij) = NULL;
1457: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,hA->ij,hdnnz,honnz);
1458: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1459: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full;
1460: } else {
1461: honnz = NULL;
1462: PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,hA->ij,hdnnz);
1463: }
1465: /* reset assembled flag and call the initialize method */
1466: hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1467: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1468: PetscStackCallStandard(HYPRE_IJMatrixInitialize,hA->ij);
1469: #else
1470: PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,hA->ij,HYPRE_MEMORY_HOST);
1471: #endif
1472: if (!dnnz) {
1473: PetscFree(hdnnz);
1474: }
1475: if (!onnz && honnz) {
1476: PetscFree(honnz);
1477: }
1478: /* Match AIJ logic */
1479: A->preallocated = PETSC_TRUE;
1480: A->assembled = PETSC_FALSE;
1481: return 0;
1482: }
1484: /*@C
1485: MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1487: Collective on Mat
1489: Input Parameters:
1490: + A - the matrix
1491: . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix
1492: (same value is used for all local rows)
1493: . dnnz - array containing the number of nonzeros in the various rows of the
1494: DIAGONAL portion of the local submatrix (possibly different for each row)
1495: or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1496: The size of this array is equal to the number of local rows, i.e 'm'.
1497: For matrices that will be factored, you must leave room for (and set)
1498: the diagonal entry even if it is zero.
1499: . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1500: submatrix (same value is used for all local rows).
1501: - onnz - array containing the number of nonzeros in the various rows of the
1502: OFF-DIAGONAL portion of the local submatrix (possibly different for
1503: each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1504: structure. The size of this array is equal to the number
1505: of local rows, i.e 'm'.
1507: Notes:
1508: If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored.
1510: Level: intermediate
1512: .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE
1513: @*/
1514: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1515: {
1518: PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1519: return 0;
1520: }
1522: /*
1523: MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix
1525: Collective
1527: Input Parameters:
1528: + parcsr - the pointer to the hypre_ParCSRMatrix
1529: . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1530: - copymode - PETSc copying options
1532: Output Parameter:
1533: . A - the matrix
1535: Level: intermediate
1537: .seealso: MatHYPRE, PetscCopyMode
1538: */
1539: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1540: {
1541: Mat T;
1542: Mat_HYPRE *hA;
1543: MPI_Comm comm;
1544: PetscInt rstart,rend,cstart,cend,M,N;
1545: PetscBool isseqaij,isseqaijmkl,ismpiaij,isaij,ishyp,isis;
1547: comm = hypre_ParCSRMatrixComm(parcsr);
1548: PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1549: PetscStrcmp(mtype,MATSEQAIJMKL,&isseqaijmkl);
1550: PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1551: PetscStrcmp(mtype,MATAIJ,&isaij);
1552: PetscStrcmp(mtype,MATHYPRE,&ishyp);
1553: PetscStrcmp(mtype,MATIS,&isis);
1554: isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1555: /* TODO */
1557: /* access ParCSRMatrix */
1558: rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1559: rend = hypre_ParCSRMatrixLastRowIndex(parcsr);
1560: cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1561: cend = hypre_ParCSRMatrixLastColDiag(parcsr);
1562: M = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1563: N = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1565: /* fix for empty local rows/columns */
1566: if (rend < rstart) rend = rstart;
1567: if (cend < cstart) cend = cstart;
1569: /* PETSc convention */
1570: rend++;
1571: cend++;
1572: rend = PetscMin(rend,M);
1573: cend = PetscMin(cend,N);
1575: /* create PETSc matrix with MatHYPRE */
1576: MatCreate(comm,&T);
1577: MatSetSizes(T,rend-rstart,cend-cstart,M,N);
1578: MatSetType(T,MATHYPRE);
1579: hA = (Mat_HYPRE*)(T->data);
1581: /* create HYPRE_IJMatrix */
1582: PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij);
1583: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,hA->ij,HYPRE_PARCSR);
1585: // TODO DEV
1586: /* create new ParCSR object if needed */
1587: if (ishyp && copymode == PETSC_COPY_VALUES) {
1588: hypre_ParCSRMatrix *new_parcsr;
1589: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
1590: hypre_CSRMatrix *hdiag,*hoffd,*ndiag,*noffd;
1592: new_parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1593: hdiag = hypre_ParCSRMatrixDiag(parcsr);
1594: hoffd = hypre_ParCSRMatrixOffd(parcsr);
1595: ndiag = hypre_ParCSRMatrixDiag(new_parcsr);
1596: noffd = hypre_ParCSRMatrixOffd(new_parcsr);
1597: PetscArraycpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag));
1598: PetscArraycpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd));
1599: #else
1600: new_parcsr = hypre_ParCSRMatrixClone(parcsr,1);
1601: #endif
1602: parcsr = new_parcsr;
1603: copymode = PETSC_OWN_POINTER;
1604: }
1606: /* set ParCSR object */
1607: hypre_IJMatrixObject(hA->ij) = parcsr;
1608: T->preallocated = PETSC_TRUE;
1610: /* set assembled flag */
1611: hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1612: #if 0
1613: PetscStackCallStandard(HYPRE_IJMatrixInitialize,hA->ij);
1614: #endif
1615: if (ishyp) {
1616: PetscMPIInt myid = 0;
1618: /* make sure we always have row_starts and col_starts available */
1619: if (HYPRE_AssumedPartitionCheck()) {
1620: MPI_Comm_rank(comm,&myid);
1621: }
1622: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1623: if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1624: PetscLayout map;
1626: MatGetLayouts(T,NULL,&map);
1627: PetscLayoutSetUp(map);
1628: hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1629: }
1630: if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1631: PetscLayout map;
1633: MatGetLayouts(T,&map,NULL);
1634: PetscLayoutSetUp(map);
1635: hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1636: }
1637: #endif
1638: /* prevent from freeing the pointer */
1639: if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1640: *A = T;
1641: MatSetOption(*A,MAT_SORTED_FULL,PETSC_TRUE);
1642: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1643: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1644: } else if (isaij) {
1645: if (copymode != PETSC_OWN_POINTER) {
1646: /* prevent from freeing the pointer */
1647: hA->inner_free = PETSC_FALSE;
1648: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1649: MatDestroy(&T);
1650: } else { /* AIJ return type with PETSC_OWN_POINTER */
1651: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1652: *A = T;
1653: }
1654: } else if (isis) {
1655: MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1656: if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1657: MatDestroy(&T);
1658: }
1659: return 0;
1660: }
1662: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1663: {
1664: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1665: HYPRE_Int type;
1668: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type);
1670: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)parcsr);
1671: return 0;
1672: }
1674: /*
1675: MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1677: Not collective
1679: Input Parameters:
1680: + A - the MATHYPRE object
1682: Output Parameter:
1683: . parcsr - the pointer to the hypre_ParCSRMatrix
1685: Level: intermediate
1687: .seealso: MatHYPRE, PetscCopyMode
1688: */
1689: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1690: {
1693: PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1694: return 0;
1695: }
1697: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1698: {
1699: hypre_ParCSRMatrix *parcsr;
1700: hypre_CSRMatrix *ha;
1701: PetscInt rst;
1704: MatGetOwnershipRange(A,&rst,NULL);
1705: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1706: if (missing) *missing = PETSC_FALSE;
1707: if (dd) *dd = -1;
1708: ha = hypre_ParCSRMatrixDiag(parcsr);
1709: if (ha) {
1710: PetscInt size,i;
1711: HYPRE_Int *ii,*jj;
1713: size = hypre_CSRMatrixNumRows(ha);
1714: ii = hypre_CSRMatrixI(ha);
1715: jj = hypre_CSRMatrixJ(ha);
1716: for (i = 0; i < size; i++) {
1717: PetscInt j;
1718: PetscBool found = PETSC_FALSE;
1720: for (j = ii[i]; j < ii[i+1] && !found; j++)
1721: found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE;
1723: if (!found) {
1724: PetscInfo(A,"Matrix is missing local diagonal entry %" PetscInt_FMT "\n",i);
1725: if (missing) *missing = PETSC_TRUE;
1726: if (dd) *dd = i+rst;
1727: return 0;
1728: }
1729: }
1730: if (!size) {
1731: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1732: if (missing) *missing = PETSC_TRUE;
1733: if (dd) *dd = rst;
1734: }
1735: } else {
1736: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1737: if (missing) *missing = PETSC_TRUE;
1738: if (dd) *dd = rst;
1739: }
1740: return 0;
1741: }
1743: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1744: {
1745: hypre_ParCSRMatrix *parcsr;
1746: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1747: hypre_CSRMatrix *ha;
1748: #endif
1749: HYPRE_Complex hs;
1751: PetscHYPREScalarCast(s,&hs);
1752: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1753: #if PETSC_PKG_HYPRE_VERSION_GE(2,19,0)
1754: PetscStackCallStandard(hypre_ParCSRMatrixScale,parcsr,hs);
1755: #else /* diagonal part */
1756: ha = hypre_ParCSRMatrixDiag(parcsr);
1757: if (ha) {
1758: PetscInt size,i;
1759: HYPRE_Int *ii;
1760: HYPRE_Complex *a;
1762: size = hypre_CSRMatrixNumRows(ha);
1763: a = hypre_CSRMatrixData(ha);
1764: ii = hypre_CSRMatrixI(ha);
1765: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1766: }
1767: /* offdiagonal part */
1768: ha = hypre_ParCSRMatrixOffd(parcsr);
1769: if (ha) {
1770: PetscInt size,i;
1771: HYPRE_Int *ii;
1772: HYPRE_Complex *a;
1774: size = hypre_CSRMatrixNumRows(ha);
1775: a = hypre_CSRMatrixData(ha);
1776: ii = hypre_CSRMatrixI(ha);
1777: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1778: }
1779: #endif
1780: return 0;
1781: }
1783: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1784: {
1785: hypre_ParCSRMatrix *parcsr;
1786: HYPRE_Int *lrows;
1787: PetscInt rst,ren,i;
1790: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1791: PetscMalloc1(numRows,&lrows);
1792: MatGetOwnershipRange(A,&rst,&ren);
1793: for (i=0;i<numRows;i++) {
1794: if (rows[i] < rst || rows[i] >= ren)
1795: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1796: lrows[i] = rows[i] - rst;
1797: }
1798: PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,parcsr,numRows,lrows);
1799: PetscFree(lrows);
1800: return 0;
1801: }
1803: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1804: {
1805: if (ha) {
1806: HYPRE_Int *ii, size;
1807: HYPRE_Complex *a;
1809: size = hypre_CSRMatrixNumRows(ha);
1810: a = hypre_CSRMatrixData(ha);
1811: ii = hypre_CSRMatrixI(ha);
1813: if (a) PetscArrayzero(a,ii[size]);
1814: }
1815: return 0;
1816: }
1818: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1819: {
1820: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1822: if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1823: PetscStackCallStandard(HYPRE_IJMatrixSetConstantValues,hA->ij,0.0);
1824: } else {
1825: hypre_ParCSRMatrix *parcsr;
1827: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1828: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1829: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1830: }
1831: return 0;
1832: }
1834: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag)
1835: {
1836: PetscInt ii;
1837: HYPRE_Int *i, *j;
1838: HYPRE_Complex *a;
1840: if (!hA) return 0;
1842: i = hypre_CSRMatrixI(hA);
1843: j = hypre_CSRMatrixJ(hA);
1844: a = hypre_CSRMatrixData(hA);
1846: for (ii = 0; ii < N; ii++) {
1847: HYPRE_Int jj, ibeg, iend, irow;
1849: irow = rows[ii];
1850: ibeg = i[irow];
1851: iend = i[irow+1];
1852: for (jj = ibeg; jj < iend; jj++)
1853: if (j[jj] == irow) a[jj] = diag;
1854: else a[jj] = 0.0;
1855: }
1856: return 0;
1857: }
1859: static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1860: {
1861: hypre_ParCSRMatrix *parcsr;
1862: PetscInt *lrows,len;
1863: HYPRE_Complex hdiag;
1866: PetscHYPREScalarCast(diag,&hdiag);
1867: /* retrieve the internal matrix */
1868: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1869: /* get locally owned rows */
1870: MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1871: /* zero diagonal part */
1872: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag);
1873: /* zero off-diagonal part */
1874: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);
1876: PetscFree(lrows);
1877: return 0;
1878: }
1880: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1881: {
1882: if (mat->nooffprocentries) return 0;
1884: MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1885: return 0;
1886: }
1888: static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1889: {
1890: hypre_ParCSRMatrix *parcsr;
1891: HYPRE_Int hnz;
1893: /* retrieve the internal matrix */
1894: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1895: /* call HYPRE API */
1896: PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v);
1897: if (nz) *nz = (PetscInt)hnz;
1898: return 0;
1899: }
1901: static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1902: {
1903: hypre_ParCSRMatrix *parcsr;
1904: HYPRE_Int hnz;
1906: /* retrieve the internal matrix */
1907: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1908: /* call HYPRE API */
1909: hnz = nz ? (HYPRE_Int)(*nz) : 0;
1910: PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v);
1911: return 0;
1912: }
1914: static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1915: {
1916: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1917: PetscInt i;
1919: if (!m || !n) return 0;
1920: /* Ignore negative row indices
1921: * And negative column indices should be automatically ignored in hypre
1922: * */
1923: for (i=0; i<m; i++) {
1924: if (idxm[i] >= 0) {
1925: HYPRE_Int hn = (HYPRE_Int)n;
1926: PetscStackCallStandard(HYPRE_IJMatrixGetValues,hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n));
1927: }
1928: }
1929: return 0;
1930: }
1932: static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg)
1933: {
1934: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1936: switch (op) {
1937: case MAT_NO_OFF_PROC_ENTRIES:
1938: if (flg) {
1939: PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,hA->ij,0);
1940: }
1941: break;
1942: case MAT_SORTED_FULL:
1943: hA->sorted_full = flg;
1944: break;
1945: default:
1946: break;
1947: }
1948: return 0;
1949: }
1951: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
1952: {
1953: PetscViewerFormat format;
1955: PetscViewerGetFormat(view,&format);
1956: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
1957: if (format != PETSC_VIEWER_NATIVE) {
1958: Mat B;
1959: hypre_ParCSRMatrix *parcsr;
1960: PetscErrorCode (*mview)(Mat,PetscViewer) = NULL;
1962: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1963: MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);
1964: MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);
1966: (*mview)(B,view);
1967: MatDestroy(&B);
1968: } else {
1969: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1970: PetscMPIInt size;
1971: PetscBool isascii;
1972: const char *filename;
1974: /* HYPRE uses only text files */
1975: PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
1977: PetscViewerFileGetName(view,&filename);
1978: PetscStackCallStandard(HYPRE_IJMatrixPrint,hA->ij,filename);
1979: MPI_Comm_size(hA->comm,&size);
1980: if (size > 1) {
1981: PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);
1982: } else {
1983: PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);
1984: }
1985: }
1986: return 0;
1987: }
1989: static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B)
1990: {
1991: hypre_ParCSRMatrix *parcsr = NULL;
1992: PetscCopyMode cpmode;
1994: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1995: if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
1996: parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1997: cpmode = PETSC_OWN_POINTER;
1998: } else {
1999: cpmode = PETSC_COPY_VALUES;
2000: }
2001: MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);
2002: return 0;
2003: }
2005: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2006: {
2007: hypre_ParCSRMatrix *acsr,*bcsr;
2009: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2010: MatHYPREGetParCSR_HYPRE(A,&acsr);
2011: MatHYPREGetParCSR_HYPRE(B,&bcsr);
2012: PetscStackCallStandard(hypre_ParCSRMatrixCopy,acsr,bcsr,1);
2013: MatSetOption(B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2014: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2015: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2016: } else {
2017: MatCopy_Basic(A,B,str);
2018: }
2019: return 0;
2020: }
2022: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2023: {
2024: hypre_ParCSRMatrix *parcsr;
2025: hypre_CSRMatrix *dmat;
2026: HYPRE_Complex *a;
2027: HYPRE_Complex *data = NULL;
2028: HYPRE_Int *diag = NULL;
2029: PetscInt i;
2030: PetscBool cong;
2032: MatHasCongruentLayouts(A,&cong);
2034: if (PetscDefined(USE_DEBUG)) {
2035: PetscBool miss;
2036: MatMissingDiagonal(A,&miss,NULL);
2038: }
2039: MatHYPREGetParCSR_HYPRE(A,&parcsr);
2040: dmat = hypre_ParCSRMatrixDiag(parcsr);
2041: if (dmat) {
2042: /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
2043: VecGetArray(d,(PetscScalar**)&a);
2044: diag = hypre_CSRMatrixI(dmat);
2045: data = hypre_CSRMatrixData(dmat);
2046: for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]];
2047: VecRestoreArray(d,(PetscScalar**)&a);
2048: }
2049: return 0;
2050: }
2052: #include <petscblaslapack.h>
2054: static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str)
2055: {
2056: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2057: {
2058: Mat B;
2059: hypre_ParCSRMatrix *x,*y,*z;
2061: MatHYPREGetParCSR(Y,&y);
2062: MatHYPREGetParCSR(X,&x);
2063: PetscStackCallStandard(hypre_ParCSRMatrixAdd,1.0,y,1.0,x,&z);
2064: MatCreateFromParCSR(z,MATHYPRE,PETSC_OWN_POINTER,&B);
2065: MatHeaderMerge(Y,&B);
2066: }
2067: #else
2068: if (str == SAME_NONZERO_PATTERN) {
2069: hypre_ParCSRMatrix *x,*y;
2070: hypre_CSRMatrix *xloc,*yloc;
2071: PetscInt xnnz,ynnz;
2072: HYPRE_Complex *xarr,*yarr;
2073: PetscBLASInt one=1,bnz;
2075: MatHYPREGetParCSR(Y,&y);
2076: MatHYPREGetParCSR(X,&x);
2078: /* diagonal block */
2079: xloc = hypre_ParCSRMatrixDiag(x);
2080: yloc = hypre_ParCSRMatrixDiag(y);
2081: xnnz = 0;
2082: ynnz = 0;
2083: xarr = NULL;
2084: yarr = NULL;
2085: if (xloc) {
2086: xarr = hypre_CSRMatrixData(xloc);
2087: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2088: }
2089: if (yloc) {
2090: yarr = hypre_CSRMatrixData(yloc);
2091: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2092: }
2094: PetscBLASIntCast(xnnz,&bnz);
2095: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2097: /* off-diagonal block */
2098: xloc = hypre_ParCSRMatrixOffd(x);
2099: yloc = hypre_ParCSRMatrixOffd(y);
2100: xnnz = 0;
2101: ynnz = 0;
2102: xarr = NULL;
2103: yarr = NULL;
2104: if (xloc) {
2105: xarr = hypre_CSRMatrixData(xloc);
2106: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2107: }
2108: if (yloc) {
2109: yarr = hypre_CSRMatrixData(yloc);
2110: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2111: }
2113: PetscBLASIntCast(xnnz,&bnz);
2114: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2115: } else if (str == SUBSET_NONZERO_PATTERN) {
2116: MatAXPY_Basic(Y,a,X,str);
2117: } else {
2118: Mat B;
2120: MatAXPY_Basic_Preallocate(Y,X,&B);
2121: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2122: MatHeaderReplace(Y,&B);
2123: }
2124: #endif
2125: return 0;
2126: }
2128: static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[])
2129: {
2130: MPI_Comm comm;
2131: PetscMPIInt size;
2132: PetscLayout rmap,cmap;
2133: Mat_HYPRE *hmat;
2134: hypre_ParCSRMatrix *parCSR;
2135: hypre_CSRMatrix *diag,*offd;
2136: Mat A,B,cooMat;
2137: PetscScalar *Aa,*Ba;
2138: HYPRE_MemoryLocation hypreMemtype = HYPRE_MEMORY_HOST;
2139: PetscMemType petscMemtype;
2140: MatType matType = MATAIJ; /* default type of cooMat */
2142: /* Build an agent matrix cooMat whose type is either MATAIJ or MATAIJKOKKOS.
2143: It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2144: */
2145: PetscObjectGetComm((PetscObject)mat,&comm);
2146: MPI_Comm_size(comm,&size);
2147: PetscLayoutSetUp(mat->rmap);
2148: PetscLayoutSetUp(mat->cmap);
2149: MatGetLayouts(mat,&rmap,&cmap);
2151: /* I do not know how hypre_ParCSRMatrix stores diagonal elements for non-square matrices, so I just give up now */
2154: #if defined(PETSC_HAVE_DEVICE)
2155: if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
2156: #if defined(PETSC_HAVE_KOKKOS)
2157: matType = MATAIJKOKKOS;
2158: #else
2159: SETERRQ(comm,PETSC_ERR_SUP,"To support MATHYPRE COO assembly on device, we need Kokkos, e.g., --download-kokkos --download-kokkos-kernels");
2160: #endif
2161: }
2162: #endif
2164: /* Do COO preallocation through cooMat */
2165: hmat = (Mat_HYPRE*)mat->data;
2166: MatDestroy(&hmat->cooMat);
2167: MatCreate(comm,&cooMat);
2168: MatSetType(cooMat,matType);
2169: MatSetLayouts(cooMat,rmap,cmap);
2170: MatSetPreallocationCOO(cooMat,coo_n,coo_i,coo_j);
2172: /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */
2173: MatSetOption(mat,MAT_SORTED_FULL,PETSC_TRUE);
2174: MatSetOption(mat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
2175: MatHYPRE_CreateFromMat(cooMat,hmat); /* Create hmat->ij and preallocate it */
2176: MatHYPRE_IJMatrixCopy(cooMat,hmat->ij); /* Copy A's (a,i,j) to hmat->ij. To reuse code. Copying 'a' is not really needed */
2178: mat->preallocated = PETSC_TRUE;
2179: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
2180: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */
2182: /* Alias cooMat's data array to IJMatrix's */
2183: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hmat->ij,(void**)&parCSR);
2184: diag = hypre_ParCSRMatrixDiag(parCSR);
2185: offd = hypre_ParCSRMatrixOffd(parCSR);
2187: hypreMemtype = hypre_CSRMatrixMemoryLocation(diag);
2188: A = (size == 1)? cooMat : ((Mat_MPIAIJ*)cooMat->data)->A;
2189: MatSeqAIJGetCSRAndMemType(A,NULL,NULL,&Aa,&petscMemtype);
2190: PetscAssert((PetscMemTypeHost(petscMemtype) && hypreMemtype == HYPRE_MEMORY_HOST) ||
2191: (PetscMemTypeDevice(petscMemtype) && hypreMemtype == HYPRE_MEMORY_DEVICE),
2192: comm,PETSC_ERR_PLIB,"PETSc and hypre's memory types mismatch");
2194: hmat->diagJ = hypre_CSRMatrixJ(diag);
2195: PetscStackCall("hypre_TFree",hypre_TFree(hypre_CSRMatrixData(diag),hypreMemtype));
2196: hypre_CSRMatrixData(diag) = (HYPRE_Complex*)Aa;
2197: hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */
2199: /* Copy diagonal pointers of A to device to facilitate MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos */
2200: if (hypreMemtype == HYPRE_MEMORY_DEVICE) {
2201: PetscStackCall("hypre_TAlloc",hmat->diag = hypre_TAlloc(PetscInt,rmap->n,hypreMemtype));
2202: MatMarkDiagonal_SeqAIJ(A); /* We need updated diagonal positions */
2203: PetscStackCall("hypre_TMemcpy",hypre_TMemcpy(hmat->diag,((Mat_SeqAIJ*)A->data)->diag,PetscInt,rmap->n,hypreMemtype,HYPRE_MEMORY_HOST));
2204: }
2206: if (size > 1) {
2207: B = ((Mat_MPIAIJ*)cooMat->data)->B;
2208: MatSeqAIJGetCSRAndMemType(B,NULL,NULL,&Ba,&petscMemtype);
2209: hmat->offdJ = hypre_CSRMatrixJ(offd);
2210: PetscStackCall("hypre_TFree",hypre_TFree(hypre_CSRMatrixData(offd),hypreMemtype));
2211: hypre_CSRMatrixData(offd) = (HYPRE_Complex*)Ba;
2212: hypre_CSRMatrixOwnsData(offd) = 0;
2213: }
2215: /* Record cooMat for use in MatSetValuesCOO_HYPRE */
2216: hmat->cooMat = cooMat;
2217: hmat->memType = hypreMemtype;
2218: return 0;
2219: }
2221: static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode)
2222: {
2223: Mat_HYPRE *hmat = (Mat_HYPRE*)mat->data;
2224: PetscMPIInt size;
2225: Mat A;
2228: MPI_Comm_size(hmat->comm,&size);
2229: MatSetValuesCOO(hmat->cooMat,v,imode);
2231: /* Move diagonal elements of the diagonal block to the front of their row, as needed by ParCSRMatrix. So damn hacky */
2232: A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ*)hmat->cooMat->data)->A;
2233: if (hmat->memType == HYPRE_MEMORY_HOST) {
2234: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
2235: PetscInt i,m,*Ai = aij->i,*Adiag = aij->diag;
2236: PetscScalar *Aa = aij->a,tmp;
2238: MatGetSize(A,&m,NULL);
2239: for (i=0; i<m; i++) {
2240: if (Adiag[i] >= Ai[i] && Adiag[i] < Ai[i+1]) { /* Digonal element of this row exists in a[] and j[] */
2241: tmp = Aa[Ai[i]];
2242: Aa[Ai[i]] = Aa[Adiag[i]];
2243: Aa[Adiag[i]] = tmp;
2244: }
2245: }
2246: } else {
2247: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2248: MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos(A,hmat->diag);
2249: #endif
2250: }
2251: return 0;
2252: }
2254: /*MC
2255: MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2256: based on the hypre IJ interface.
2258: Level: intermediate
2260: .seealso: MatCreate()
2261: M*/
2263: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2264: {
2265: Mat_HYPRE *hB;
2267: PetscNewLog(B,&hB);
2269: hB->inner_free = PETSC_TRUE;
2270: hB->available = PETSC_TRUE;
2271: hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */
2272: hB->size = 0;
2273: hB->array = NULL;
2275: B->data = (void*)hB;
2276: B->assembled = PETSC_FALSE;
2278: PetscMemzero(B->ops,sizeof(struct _MatOps));
2279: B->ops->mult = MatMult_HYPRE;
2280: B->ops->multtranspose = MatMultTranspose_HYPRE;
2281: B->ops->multadd = MatMultAdd_HYPRE;
2282: B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
2283: B->ops->setup = MatSetUp_HYPRE;
2284: B->ops->destroy = MatDestroy_HYPRE;
2285: B->ops->assemblyend = MatAssemblyEnd_HYPRE;
2286: B->ops->assemblybegin = MatAssemblyBegin_HYPRE;
2287: B->ops->setvalues = MatSetValues_HYPRE;
2288: B->ops->missingdiagonal = MatMissingDiagonal_HYPRE;
2289: B->ops->scale = MatScale_HYPRE;
2290: B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
2291: B->ops->zeroentries = MatZeroEntries_HYPRE;
2292: B->ops->zerorows = MatZeroRows_HYPRE;
2293: B->ops->getrow = MatGetRow_HYPRE;
2294: B->ops->restorerow = MatRestoreRow_HYPRE;
2295: B->ops->getvalues = MatGetValues_HYPRE;
2296: B->ops->setoption = MatSetOption_HYPRE;
2297: B->ops->duplicate = MatDuplicate_HYPRE;
2298: B->ops->copy = MatCopy_HYPRE;
2299: B->ops->view = MatView_HYPRE;
2300: B->ops->getdiagonal = MatGetDiagonal_HYPRE;
2301: B->ops->axpy = MatAXPY_HYPRE;
2302: B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2303: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2304: B->ops->bindtocpu = MatBindToCPU_HYPRE;
2305: B->boundtocpu = PETSC_FALSE;
2306: #endif
2308: /* build cache for off array entries formed */
2309: MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
2311: PetscCommGetComm(PetscObjectComm((PetscObject)B),&hB->comm);
2312: PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
2313: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
2314: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
2315: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_hypre_C",MatProductSetFromOptions_HYPRE);
2316: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_mpiaij_hypre_C",MatProductSetFromOptions_HYPRE);
2317: PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
2318: PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
2319: PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_HYPRE);
2320: PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_HYPRE);
2321: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2322: #if defined(HYPRE_USING_HIP)
2323: PetscDeviceInitialize(PETSC_DEVICE_HIP);
2324: MatSetVecType(B,VECHIP);
2325: #endif
2326: #if defined(HYPRE_USING_CUDA)
2327: PetscDeviceInitialize(PETSC_DEVICE_CUDA);
2328: MatSetVecType(B,VECCUDA);
2329: #endif
2330: #endif
2331: return 0;
2332: }
2334: static PetscErrorCode hypre_array_destroy(void *ptr)
2335: {
2336: hypre_TFree(ptr,HYPRE_MEMORY_HOST);
2337: return 0;
2338: }