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 <../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: static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
23: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
24: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
25: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
26: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,HYPRE_Complex,Vec,HYPRE_Complex,Vec,PetscBool);
27: static PetscErrorCode hypre_array_destroy(void*);
28: static PetscErrorCode MatSetValues_HYPRE(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode ins);
30: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
31: {
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;
39: if (A_d) { /* determine number of nonzero entries in local diagonal part */
40: MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);
41: if (done_d) {
42: PetscMalloc1(n_d,&nnz_d);
43: for (i=0; i<n_d; i++) {
44: nnz_d[i] = ia_d[i+1] - ia_d[i];
45: }
46: }
47: 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: MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
51: if (done_o) {
52: PetscMalloc1(n_o,&nnz_o);
53: for (i=0; i<n_o; i++) {
54: nnz_o[i] = ia_o[i+1] - ia_o[i];
55: }
56: }
57: MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
58: }
59: if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */
60: if (!done_o) { /* only diagonal part */
61: PetscCalloc1(n_d,&nnz_o);
62: }
63: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
64: { /* If we don't do this, the columns of the matrix will be all zeros! */
65: hypre_AuxParCSRMatrix *aux_matrix;
66: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
67: hypre_AuxParCSRMatrixDestroy(aux_matrix);
68: hypre_IJMatrixTranslator(ij) = NULL;
69: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,nnz_d,nnz_o));
70: /* it seems they partially fixed it in 2.19.0 */
71: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
72: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
73: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
74: #endif
75: }
76: #else
77: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,nnz_d,nnz_o));
78: #endif
79: PetscFree(nnz_d);
80: PetscFree(nnz_o);
81: }
82: return(0);
83: }
85: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
86: {
88: PetscInt rstart,rend,cstart,cend;
91: PetscLayoutSetUp(A->rmap);
92: PetscLayoutSetUp(A->cmap);
93: rstart = A->rmap->rstart;
94: rend = A->rmap->rend;
95: cstart = A->cmap->rstart;
96: cend = A->cmap->rend;
97: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
98: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
99: {
100: PetscBool same;
101: Mat A_d,A_o;
102: const PetscInt *colmap;
103: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&same);
104: if (same) {
105: MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);
106: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
107: return(0);
108: }
109: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&same);
110: if (same) {
111: MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);
112: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
113: return(0);
114: }
115: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&same);
116: if (same) {
117: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
118: return(0);
119: }
120: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQBAIJ,&same);
121: if (same) {
122: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
123: return(0);
124: }
125: }
126: return(0);
127: }
129: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
130: {
131: PetscErrorCode ierr;
132: PetscInt i,rstart,rend,ncols,nr,nc;
133: const PetscScalar *values;
134: const PetscInt *cols;
135: PetscBool flg;
138: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
139: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
140: #else
141: PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,(ij,HYPRE_MEMORY_HOST));
142: #endif
143: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
144: MatGetSize(A,&nr,&nc);
145: if (flg && nr == nc) {
146: MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
147: return(0);
148: }
149: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
150: if (flg) {
151: MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
152: return(0);
153: }
155: MatGetOwnershipRange(A,&rstart,&rend);
156: for (i=rstart; i<rend; i++) {
157: MatGetRow(A,i,&ncols,&cols,&values);
158: if (ncols) {
159: HYPRE_Int nc = (HYPRE_Int)ncols;
161: if ((PetscInt)nc != ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",ncols,i);
162: PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,&nc,(HYPRE_BigInt *)&i,(HYPRE_BigInt *)cols,(HYPRE_Complex *)values));
163: }
164: MatRestoreRow(A,i,&ncols,&cols,&values);
165: }
166: return(0);
167: }
169: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
170: {
171: PetscErrorCode ierr;
172: Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data;
173: HYPRE_Int type;
174: hypre_ParCSRMatrix *par_matrix;
175: hypre_AuxParCSRMatrix *aux_matrix;
176: hypre_CSRMatrix *hdiag;
177: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
178: const PetscScalar *pa;
181: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
182: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
183: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
184: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
185: /*
186: this is the Hack part where we monkey directly with the hypre datastructures
187: */
188: if (sameint) {
189: PetscArraycpy(hdiag->i,pdiag->i,A->rmap->n + 1);
190: PetscArraycpy(hdiag->j,pdiag->j,pdiag->nz);
191: } else {
192: PetscInt i;
194: for (i=0;i<A->rmap->n + 1;i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
195: for (i=0;i<pdiag->nz;i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
196: }
198: MatSeqAIJGetArrayRead(A,&pa);
199: PetscArraycpy(hdiag->data,pa,pdiag->nz);
200: MatSeqAIJRestoreArrayRead(A,&pa);
202: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
203: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
204: return(0);
205: }
207: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
208: {
209: PetscErrorCode ierr;
210: Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data;
211: Mat_SeqAIJ *pdiag,*poffd;
212: PetscInt i,*garray = pA->garray,*jj,cstart,*pjj;
213: HYPRE_Int *hjj,type;
214: hypre_ParCSRMatrix *par_matrix;
215: hypre_AuxParCSRMatrix *aux_matrix;
216: hypre_CSRMatrix *hdiag,*hoffd;
217: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
218: const PetscScalar *pa;
221: pdiag = (Mat_SeqAIJ*) pA->A->data;
222: poffd = (Mat_SeqAIJ*) pA->B->data;
223: /* cstart is only valid for square MPIAIJ layed out in the usual way */
224: MatGetOwnershipRange(A,&cstart,NULL);
226: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
227: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
228: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
229: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
230: hoffd = hypre_ParCSRMatrixOffd(par_matrix);
232: /*
233: this is the Hack part where we monkey directly with the hypre datastructures
234: */
235: if (sameint) {
236: PetscArraycpy(hdiag->i,pdiag->i,pA->A->rmap->n + 1);
237: } else {
238: for (i=0; i<pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)(pdiag->i[i]);
239: }
240: /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
241: hjj = hdiag->j;
242: pjj = pdiag->j;
243: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
244: for (i=0; i<pdiag->nz; i++) hjj[i] = pjj[i];
245: #else
246: for (i=0; i<pdiag->nz; i++) hjj[i] = cstart + pjj[i];
247: #endif
248: MatSeqAIJGetArrayRead(pA->A,&pa);
249: PetscArraycpy(hdiag->data,pa,pdiag->nz);
250: MatSeqAIJRestoreArrayRead(pA->A,&pa);
251: if (sameint) {
252: PetscArraycpy(hoffd->i,poffd->i,pA->A->rmap->n + 1);
253: } else {
254: for (i=0; i<pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)(poffd->i[i]);
255: }
257: /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
258: If we hacked a hypre a bit more we might be able to avoid this step */
259: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
260: PetscStackCallStandard(hypre_CSRMatrixBigInitialize,(hoffd));
261: jj = (PetscInt*) hoffd->big_j;
262: #else
263: jj = (PetscInt*) hoffd->j;
264: #endif
265: pjj = poffd->j;
266: for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
268: MatSeqAIJGetArrayRead(pA->B,&pa);
269: PetscArraycpy(hoffd->data,pa,poffd->nz);
270: MatSeqAIJRestoreArrayRead(pA->B,&pa);
272: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
273: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
274: return(0);
275: }
277: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
278: {
279: Mat_HYPRE* mhA = (Mat_HYPRE*)(A->data);
280: Mat lA;
281: ISLocalToGlobalMapping rl2g,cl2g;
282: IS is;
283: hypre_ParCSRMatrix *hA;
284: hypre_CSRMatrix *hdiag,*hoffd;
285: MPI_Comm comm;
286: HYPRE_Complex *hdd,*hod,*aa;
287: PetscScalar *data;
288: HYPRE_BigInt *col_map_offd;
289: HYPRE_Int *hdi,*hdj,*hoi,*hoj;
290: PetscInt *ii,*jj,*iptr,*jptr;
291: PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
292: HYPRE_Int type;
293: PetscErrorCode ierr;
296: comm = PetscObjectComm((PetscObject)A);
297: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
298: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
299: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
300: M = hypre_ParCSRMatrixGlobalNumRows(hA);
301: N = hypre_ParCSRMatrixGlobalNumCols(hA);
302: str = hypre_ParCSRMatrixFirstRowIndex(hA);
303: stc = hypre_ParCSRMatrixFirstColDiag(hA);
304: hdiag = hypre_ParCSRMatrixDiag(hA);
305: hoffd = hypre_ParCSRMatrixOffd(hA);
306: dr = hypre_CSRMatrixNumRows(hdiag);
307: dc = hypre_CSRMatrixNumCols(hdiag);
308: nnz = hypre_CSRMatrixNumNonzeros(hdiag);
309: hdi = hypre_CSRMatrixI(hdiag);
310: hdj = hypre_CSRMatrixJ(hdiag);
311: hdd = hypre_CSRMatrixData(hdiag);
312: oc = hypre_CSRMatrixNumCols(hoffd);
313: nnz += hypre_CSRMatrixNumNonzeros(hoffd);
314: hoi = hypre_CSRMatrixI(hoffd);
315: hoj = hypre_CSRMatrixJ(hoffd);
316: hod = hypre_CSRMatrixData(hoffd);
317: if (reuse != MAT_REUSE_MATRIX) {
318: PetscInt *aux;
320: /* generate l2g maps for rows and cols */
321: ISCreateStride(comm,dr,str,1,&is);
322: ISLocalToGlobalMappingCreateIS(is,&rl2g);
323: ISDestroy(&is);
324: col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
325: PetscMalloc1(dc+oc,&aux);
326: for (i=0; i<dc; i++) aux[i] = i+stc;
327: for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
328: ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
329: ISLocalToGlobalMappingCreateIS(is,&cl2g);
330: ISDestroy(&is);
331: /* create MATIS object */
332: MatCreate(comm,B);
333: MatSetSizes(*B,dr,dc,M,N);
334: MatSetType(*B,MATIS);
335: MatSetLocalToGlobalMapping(*B,rl2g,cl2g);
336: ISLocalToGlobalMappingDestroy(&rl2g);
337: ISLocalToGlobalMappingDestroy(&cl2g);
339: /* allocate CSR for local matrix */
340: PetscMalloc1(dr+1,&iptr);
341: PetscMalloc1(nnz,&jptr);
342: PetscMalloc1(nnz,&data);
343: } else {
344: PetscInt nr;
345: PetscBool done;
346: MatISGetLocalMat(*B,&lA);
347: MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);
348: if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
349: if (iptr[nr] < nnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in local mat! reuse %D requested %D",iptr[nr],nnz);
350: MatSeqAIJGetArray(lA,&data);
351: }
352: /* merge local matrices */
353: ii = iptr;
354: jj = jptr;
355: 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++ */
356: *ii = *(hdi++) + *(hoi++);
357: for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
358: PetscScalar *aold = (PetscScalar*)aa;
359: PetscInt *jold = jj,nc = jd+jo;
360: for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; }
361: for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
362: *(++ii) = *(hdi++) + *(hoi++);
363: PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);
364: }
365: for (; cum<dr; cum++) *(++ii) = nnz;
366: if (reuse != MAT_REUSE_MATRIX) {
367: Mat_SeqAIJ* a;
369: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);
370: MatISSetLocalMat(*B,lA);
371: /* hack SeqAIJ */
372: a = (Mat_SeqAIJ*)(lA->data);
373: a->free_a = PETSC_TRUE;
374: a->free_ij = PETSC_TRUE;
375: MatDestroy(&lA);
376: }
377: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
378: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
379: if (reuse == MAT_INPLACE_MATRIX) {
380: MatHeaderReplace(A,B);
381: }
382: return(0);
383: }
385: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
386: {
387: Mat M = NULL;
388: Mat_HYPRE *hB;
389: MPI_Comm comm = PetscObjectComm((PetscObject)A);
393: if (reuse == MAT_REUSE_MATRIX) {
394: /* always destroy the old matrix and create a new memory;
395: hope this does not churn the memory too much. The problem
396: is I do not know if it is possible to put the matrix back to
397: its initial state so that we can directly copy the values
398: the second time through. */
399: hB = (Mat_HYPRE*)((*B)->data);
400: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
401: } else {
402: MatCreate(comm,&M);
403: MatSetType(M,MATHYPRE);
404: MatSetSizes(M,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
405: hB = (Mat_HYPRE*)(M->data);
406: if (reuse == MAT_INITIAL_MATRIX) *B = M;
407: }
408: MatSetOption(*B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
409: MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
410: MatHYPRE_CreateFromMat(A,hB);
411: MatHYPRE_IJMatrixCopy(A,hB->ij);
412: if (reuse == MAT_INPLACE_MATRIX) {
413: MatHeaderReplace(A,&M);
414: }
415: (*B)->preallocated = PETSC_TRUE;
416: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
417: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
418: return(0);
419: }
421: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
422: {
423: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
424: hypre_ParCSRMatrix *parcsr;
425: hypre_CSRMatrix *hdiag,*hoffd;
426: MPI_Comm comm;
427: PetscScalar *da,*oa,*aptr;
428: PetscInt *dii,*djj,*oii,*ojj,*iptr;
429: PetscInt i,dnnz,onnz,m,n;
430: HYPRE_Int type;
431: PetscMPIInt size;
432: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
433: PetscErrorCode ierr;
436: comm = PetscObjectComm((PetscObject)A);
437: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
438: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
439: if (reuse == MAT_REUSE_MATRIX) {
440: PetscBool ismpiaij,isseqaij;
441: PetscObjectBaseTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
442: PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
443: if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
444: }
445: #if defined(PETSC_HAVE_HYPRE_DEVICE)
446: if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) SETERRQ(comm,PETSC_ERR_SUP,"Not yet implemented");
447: #endif
448: MPI_Comm_size(comm,&size);
450: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
451: hdiag = hypre_ParCSRMatrixDiag(parcsr);
452: hoffd = hypre_ParCSRMatrixOffd(parcsr);
453: m = hypre_CSRMatrixNumRows(hdiag);
454: n = hypre_CSRMatrixNumCols(hdiag);
455: dnnz = hypre_CSRMatrixNumNonzeros(hdiag);
456: onnz = hypre_CSRMatrixNumNonzeros(hoffd);
457: if (reuse == MAT_INITIAL_MATRIX) {
458: PetscMalloc1(m+1,&dii);
459: PetscMalloc1(dnnz,&djj);
460: PetscMalloc1(dnnz,&da);
461: } else if (reuse == MAT_REUSE_MATRIX) {
462: PetscInt nr;
463: PetscBool done;
464: if (size > 1) {
465: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
467: MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
468: if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows in diag part! %D != %D",nr,m);
469: if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in diag part! reuse %D hypre %D",dii[nr],dnnz);
470: MatSeqAIJGetArray(b->A,&da);
471: } else {
472: MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
473: if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
474: if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros! reuse %D hypre %D",dii[nr],dnnz);
475: MatSeqAIJGetArray(*B,&da);
476: }
477: } else { /* MAT_INPLACE_MATRIX */
478: if (!sameint) {
479: PetscMalloc1(m+1,&dii);
480: PetscMalloc1(dnnz,&djj);
481: } else {
482: dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
483: djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
484: }
485: da = (PetscScalar*)hypre_CSRMatrixData(hdiag);
486: }
488: if (!sameint) {
489: for (i=0;i<m+1;i++) dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]);
490: for (i=0;i<dnnz;i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]);
491: } else {
492: PetscArraycpy(dii,hypre_CSRMatrixI(hdiag),m+1);
493: PetscArraycpy(djj,hypre_CSRMatrixJ(hdiag),dnnz);
494: }
495: PetscArraycpy(da,hypre_CSRMatrixData(hdiag),dnnz);
496: iptr = djj;
497: aptr = da;
498: for (i=0; i<m; i++) {
499: PetscInt nc = dii[i+1]-dii[i];
500: PetscSortIntWithScalarArray(nc,iptr,aptr);
501: iptr += nc;
502: aptr += nc;
503: }
504: if (size > 1) {
505: HYPRE_BigInt *coffd;
506: HYPRE_Int *offdj;
508: if (reuse == MAT_INITIAL_MATRIX) {
509: PetscMalloc1(m+1,&oii);
510: PetscMalloc1(onnz,&ojj);
511: PetscMalloc1(onnz,&oa);
512: } else if (reuse == MAT_REUSE_MATRIX) {
513: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
514: PetscInt nr,hr = hypre_CSRMatrixNumRows(hoffd);
515: PetscBool done;
517: MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
518: if (nr != hr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows in offdiag part! %D != %D",nr,hr);
519: if (oii[nr] < onnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in offdiag part! reuse %D hypre %D",oii[nr],onnz);
520: MatSeqAIJGetArray(b->B,&oa);
521: } else { /* MAT_INPLACE_MATRIX */
522: if (!sameint) {
523: PetscMalloc1(m+1,&oii);
524: PetscMalloc1(onnz,&ojj);
525: } else {
526: oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
527: ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
528: }
529: oa = (PetscScalar*)hypre_CSRMatrixData(hoffd);
530: }
531: if (!sameint) {
532: for (i=0;i<m+1;i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]);
533: } else {
534: PetscArraycpy(oii,hypre_CSRMatrixI(hoffd),m+1);
535: }
536: offdj = hypre_CSRMatrixJ(hoffd);
537: coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
538: for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
539: PetscArraycpy(oa,hypre_CSRMatrixData(hoffd),onnz);
540: iptr = ojj;
541: aptr = oa;
542: for (i=0; i<m; i++) {
543: PetscInt nc = oii[i+1]-oii[i];
544: PetscSortIntWithScalarArray(nc,iptr,aptr);
545: iptr += nc;
546: aptr += nc;
547: }
548: if (reuse == MAT_INITIAL_MATRIX) {
549: Mat_MPIAIJ *b;
550: Mat_SeqAIJ *d,*o;
552: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
553: /* hack MPIAIJ */
554: b = (Mat_MPIAIJ*)((*B)->data);
555: d = (Mat_SeqAIJ*)b->A->data;
556: o = (Mat_SeqAIJ*)b->B->data;
557: d->free_a = PETSC_TRUE;
558: d->free_ij = PETSC_TRUE;
559: o->free_a = PETSC_TRUE;
560: o->free_ij = PETSC_TRUE;
561: } else if (reuse == MAT_INPLACE_MATRIX) {
562: Mat T;
564: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
565: if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
566: hypre_CSRMatrixI(hdiag) = NULL;
567: hypre_CSRMatrixJ(hdiag) = NULL;
568: hypre_CSRMatrixI(hoffd) = NULL;
569: hypre_CSRMatrixJ(hoffd) = NULL;
570: } else { /* Hack MPIAIJ -> free ij but not a */
571: Mat_MPIAIJ *b = (Mat_MPIAIJ*)(T->data);
572: Mat_SeqAIJ *d = (Mat_SeqAIJ*)(b->A->data);
573: Mat_SeqAIJ *o = (Mat_SeqAIJ*)(b->B->data);
575: d->free_ij = PETSC_TRUE;
576: o->free_ij = PETSC_TRUE;
577: }
578: hypre_CSRMatrixData(hdiag) = NULL;
579: hypre_CSRMatrixData(hoffd) = NULL;
580: MatHeaderReplace(A,&T);
581: }
582: } else {
583: oii = NULL;
584: ojj = NULL;
585: oa = NULL;
586: if (reuse == MAT_INITIAL_MATRIX) {
587: Mat_SeqAIJ* b;
589: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
590: /* hack SeqAIJ */
591: b = (Mat_SeqAIJ*)((*B)->data);
592: b->free_a = PETSC_TRUE;
593: b->free_ij = PETSC_TRUE;
594: } else if (reuse == MAT_INPLACE_MATRIX) {
595: Mat T;
597: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
598: if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
599: hypre_CSRMatrixI(hdiag) = NULL;
600: hypre_CSRMatrixJ(hdiag) = NULL;
601: } else { /* free ij but not a */
602: Mat_SeqAIJ* b = (Mat_SeqAIJ*)(T->data);
604: b->free_ij = PETSC_TRUE;
605: }
606: hypre_CSRMatrixData(hdiag) = NULL;
607: MatHeaderReplace(A,&T);
608: }
609: }
611: /* we have to use hypre_Tfree to free the HYPRE arrays
612: that PETSc now onws */
613: if (reuse == MAT_INPLACE_MATRIX) {
614: PetscInt nh;
615: void *ptrs[6] = {da,oa,dii,djj,oii,ojj};
616: const char *names[6] = {"_hypre_csr_da",
617: "_hypre_csr_oa",
618: "_hypre_csr_dii",
619: "_hypre_csr_djj",
620: "_hypre_csr_oii",
621: "_hypre_csr_ojj"};
622: nh = sameint ? 6 : 2;
623: for (i=0; i<nh; i++) {
624: PetscContainer c;
626: PetscContainerCreate(comm,&c);
627: PetscContainerSetPointer(c,ptrs[i]);
628: PetscContainerSetUserDestroy(c,hypre_array_destroy);
629: PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
630: PetscContainerDestroy(&c);
631: }
632: }
633: return(0);
634: }
636: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
637: {
638: hypre_ParCSRMatrix *tA;
639: hypre_CSRMatrix *hdiag,*hoffd;
640: Mat_SeqAIJ *diag,*offd;
641: PetscInt *garray,i,noffd,dnnz,onnz,*row_starts,*col_starts;
642: MPI_Comm comm = PetscObjectComm((PetscObject)A);
643: PetscBool ismpiaij,isseqaij;
644: PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
645: PetscErrorCode ierr;
646: HYPRE_Int *hdi = NULL,*hdj = NULL,*hoi = NULL,*hoj = NULL;
647: PetscInt *pdi,*pdj,*poi,*poj;
648: #if defined(PETSC_HAVE_HYPRE_DEVICE)
649: PetscBool iscuda = PETSC_FALSE;
650: #endif
653: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
654: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
655: if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type_name);
656: if (ismpiaij) {
657: Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);
659: diag = (Mat_SeqAIJ*)a->A->data;
660: offd = (Mat_SeqAIJ*)a->B->data;
661: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE) && defined(HYPRE_USING_CUDA)
662: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJCUSPARSE,&iscuda);
663: if (iscuda && !A->boundtocpu) {
664: sameint = PETSC_TRUE;
665: MatSeqAIJCUSPARSEGetIJ(a->A,PETSC_FALSE,(const HYPRE_Int**)&hdi,(const HYPRE_Int**)&hdj);
666: MatSeqAIJCUSPARSEGetIJ(a->B,PETSC_FALSE,(const HYPRE_Int**)&hoi,(const HYPRE_Int**)&hoj);
667: } else {
668: #else
669: {
670: #endif
671: pdi = diag->i;
672: pdj = diag->j;
673: poi = offd->i;
674: poj = offd->j;
675: if (sameint) {
676: hdi = (HYPRE_Int*)pdi;
677: hdj = (HYPRE_Int*)pdj;
678: hoi = (HYPRE_Int*)poi;
679: hoj = (HYPRE_Int*)poj;
680: }
681: }
682: garray = a->garray;
683: noffd = a->B->cmap->N;
684: dnnz = diag->nz;
685: onnz = offd->nz;
686: } else {
687: diag = (Mat_SeqAIJ*)A->data;
688: offd = NULL;
689: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE)
690: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJCUSPARSE,&iscuda);
691: if (iscuda && !A->boundtocpu) {
692: sameint = PETSC_TRUE;
693: MatSeqAIJCUSPARSEGetIJ(A,PETSC_FALSE,(const HYPRE_Int**)&hdi,(const HYPRE_Int**)&hdj);
694: } else {
695: #else
696: {
697: #endif
698: pdi = diag->i;
699: pdj = diag->j;
700: if (sameint) {
701: hdi = (HYPRE_Int*)pdi;
702: hdj = (HYPRE_Int*)pdj;
703: }
704: }
705: garray = NULL;
706: noffd = 0;
707: dnnz = diag->nz;
708: onnz = 0;
709: }
711: /* create a temporary ParCSR */
712: if (HYPRE_AssumedPartitionCheck()) {
713: PetscMPIInt myid;
715: MPI_Comm_rank(comm,&myid);
716: row_starts = A->rmap->range + myid;
717: col_starts = A->cmap->range + myid;
718: } else {
719: row_starts = A->rmap->range;
720: col_starts = A->cmap->range;
721: }
722: tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_BigInt*)row_starts,(HYPRE_BigInt*)col_starts,noffd,dnnz,onnz);
723: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
724: hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
725: hypre_ParCSRMatrixSetColStartsOwner(tA,0);
726: #endif
728: /* set diagonal part */
729: hdiag = hypre_ParCSRMatrixDiag(tA);
730: if (!sameint) { /* malloc CSR pointers */
731: PetscMalloc2(A->rmap->n+1,&hdi,dnnz,&hdj);
732: for (i = 0; i < A->rmap->n+1; i++) hdi[i] = (HYPRE_Int)(pdi[i]);
733: for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)(pdj[i]);
734: }
735: hypre_CSRMatrixI(hdiag) = hdi;
736: hypre_CSRMatrixJ(hdiag) = hdj;
737: hypre_CSRMatrixData(hdiag) = (HYPRE_Complex*)diag->a;
738: hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
739: hypre_CSRMatrixSetRownnz(hdiag);
740: hypre_CSRMatrixSetDataOwner(hdiag,0);
742: /* set offdiagonal part */
743: hoffd = hypre_ParCSRMatrixOffd(tA);
744: if (offd) {
745: if (!sameint) { /* malloc CSR pointers */
746: PetscMalloc2(A->rmap->n+1,&hoi,onnz,&hoj);
747: for (i = 0; i < A->rmap->n+1; i++) hoi[i] = (HYPRE_Int)(poi[i]);
748: for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)(poj[i]);
749: }
750: hypre_CSRMatrixI(hoffd) = hoi;
751: hypre_CSRMatrixJ(hoffd) = hoj;
752: hypre_CSRMatrixData(hoffd) = (HYPRE_Complex*)offd->a;
753: hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
754: hypre_CSRMatrixSetRownnz(hoffd);
755: hypre_CSRMatrixSetDataOwner(hoffd,0);
756: }
757: #if defined(PETSC_HAVE_HYPRE_DEVICE)
758: PetscStackCallStandard(hypre_ParCSRMatrixInitialize_v2,(tA,iscuda ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST));
759: #else
760: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
761: PetscStackCallStandard(hypre_ParCSRMatrixInitialize,(tA));
762: #else
763: PetscStackCallStandard(hypre_ParCSRMatrixInitialize_v2,(tA,HYPRE_MEMORY_HOST));
764: #endif
765: #endif
766: hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA),HYPRE_MEMORY_HOST);
767: hypre_ParCSRMatrixSetNumNonzeros(tA);
768: hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt*)garray;
769: if (!hypre_ParCSRMatrixCommPkg(tA)) PetscStackCallStandard(hypre_MatvecCommPkgCreate,(tA));
770: *hA = tA;
771: return(0);
772: }
774: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
775: {
776: hypre_CSRMatrix *hdiag,*hoffd;
777: PetscBool ismpiaij,sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
778: PetscErrorCode ierr;
779: #if defined(PETSC_HAVE_HYPRE_DEVICE)
780: PetscBool iscuda = PETSC_FALSE;
781: #endif
784: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
785: #if defined(PETSC_HAVE_HYPRE_DEVICE)
786: PetscObjectTypeCompareAny((PetscObject)A,&iscuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
787: if (iscuda) sameint = PETSC_TRUE;
788: #endif
789: hdiag = hypre_ParCSRMatrixDiag(*hA);
790: hoffd = hypre_ParCSRMatrixOffd(*hA);
791: /* free temporary memory allocated by PETSc
792: set pointers to NULL before destroying tA */
793: if (!sameint) {
794: HYPRE_Int *hi,*hj;
796: hi = hypre_CSRMatrixI(hdiag);
797: hj = hypre_CSRMatrixJ(hdiag);
798: PetscFree2(hi,hj);
799: if (ismpiaij) {
800: hi = hypre_CSRMatrixI(hoffd);
801: hj = hypre_CSRMatrixJ(hoffd);
802: PetscFree2(hi,hj);
803: }
804: }
805: hypre_CSRMatrixI(hdiag) = NULL;
806: hypre_CSRMatrixJ(hdiag) = NULL;
807: hypre_CSRMatrixData(hdiag) = NULL;
808: if (ismpiaij) {
809: hypre_CSRMatrixI(hoffd) = NULL;
810: hypre_CSRMatrixJ(hoffd) = NULL;
811: hypre_CSRMatrixData(hoffd) = NULL;
812: }
813: hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
814: hypre_ParCSRMatrixDestroy(*hA);
815: *hA = NULL;
816: return(0);
817: }
819: /* calls RAP from BoomerAMG:
820: the resulting ParCSR will not own the column and row starts
821: It looks like we don't need to have the diagonal entries ordered first */
822: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
823: {
824: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
825: HYPRE_Int P_owns_col_starts,R_owns_row_starts;
826: #endif
829: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
830: P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
831: R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
832: #endif
833: /* can be replaced by version test later */
834: #if defined(PETSC_HAVE_HYPRE_DEVICE)
835: PetscStackPush("hypre_ParCSRMatrixRAP");
836: *hRAP = hypre_ParCSRMatrixRAP(hR,hA,hP);
837: PetscStackPop;
838: #else
839: PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP));
840: PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP));
841: #endif
842: /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
843: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
844: hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
845: hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
846: if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
847: if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
848: #endif
849: return(0);
850: }
852: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
853: {
854: Mat B;
855: hypre_ParCSRMatrix *hA,*hP,*hPtAP = NULL;
856: PetscErrorCode ierr;
857: 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: {
878: MatSetType(C,MATAIJ);
879: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
880: C->ops->productnumeric = MatProductNumeric_PtAP;
881: return(0);
882: }
884: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
885: {
886: Mat B;
887: Mat_HYPRE *hP;
888: hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr = NULL;
889: HYPRE_Int type;
890: MPI_Comm comm = PetscObjectComm((PetscObject)A);
891: PetscBool ishypre;
892: PetscErrorCode ierr;
895: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
896: if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
897: hP = (Mat_HYPRE*)P->data;
898: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
899: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
900: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
902: MatAIJGetParCSR_Private(A,&hA);
903: MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
904: MatAIJRestoreParCSR_Private(A,&hA);
906: /* create temporary matrix and merge to C */
907: MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
908: MatHeaderMerge(C,&B);
909: return(0);
910: }
912: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
913: {
914: Mat B;
915: hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr = NULL;
916: Mat_HYPRE *hA,*hP;
917: PetscBool ishypre;
918: HYPRE_Int type;
919: PetscErrorCode ierr;
922: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
923: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
924: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
925: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
926: hA = (Mat_HYPRE*)A->data;
927: hP = (Mat_HYPRE*)P->data;
928: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
929: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
930: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
931: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
932: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
933: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
934: MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
935: MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
936: MatHeaderMerge(C,&B);
937: return(0);
938: }
940: /* calls hypre_ParMatmul
941: hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
942: hypre_ParMatrixCreate does not duplicate the communicator
943: It looks like we don't need to have the diagonal entries ordered first */
944: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
945: {
947: /* can be replaced by version test later */
948: #if defined(PETSC_HAVE_HYPRE_DEVICE)
949: PetscStackPush("hypre_ParCSRMatMat");
950: *hAB = hypre_ParCSRMatMat(hA,hB);
951: #else
952: PetscStackPush("hypre_ParMatmul");
953: *hAB = hypre_ParMatmul(hA,hB);
954: #endif
955: PetscStackPop;
956: return(0);
957: }
959: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
960: {
961: Mat D;
962: hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
963: PetscErrorCode ierr;
964: Mat_Product *product=C->product;
967: MatAIJGetParCSR_Private(A,&hA);
968: MatAIJGetParCSR_Private(B,&hB);
969: MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
970: MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);
972: MatHeaderMerge(C,&D);
973: C->product = product;
975: MatAIJRestoreParCSR_Private(A,&hA);
976: MatAIJRestoreParCSR_Private(B,&hB);
977: return(0);
978: }
980: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat C)
981: {
985: MatSetType(C,MATAIJ);
986: C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
987: C->ops->productnumeric = MatProductNumeric_AB;
988: return(0);
989: }
991: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
992: {
993: Mat D;
994: hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
995: Mat_HYPRE *hA,*hB;
996: PetscBool ishypre;
997: HYPRE_Int type;
998: PetscErrorCode ierr;
999: Mat_Product *product;
1002: PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
1003: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
1004: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1005: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
1006: hA = (Mat_HYPRE*)A->data;
1007: hB = (Mat_HYPRE*)B->data;
1008: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1009: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
1010: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
1011: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
1012: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
1013: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
1014: MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
1015: MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);
1017: /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1018: product = C->product; /* save it from MatHeaderReplace() */
1019: C->product = NULL;
1020: MatHeaderReplace(C,&D);
1021: C->product = product;
1022: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1023: C->ops->productnumeric = MatProductNumeric_AB;
1024: return(0);
1025: }
1027: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
1028: {
1029: Mat E;
1030: hypre_ParCSRMatrix *hA,*hB,*hC,*hABC = NULL;
1031: PetscErrorCode ierr;
1034: MatAIJGetParCSR_Private(A,&hA);
1035: MatAIJGetParCSR_Private(B,&hB);
1036: MatAIJGetParCSR_Private(C,&hC);
1037: MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
1038: MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
1039: MatHeaderMerge(D,&E);
1040: MatAIJRestoreParCSR_Private(A,&hA);
1041: MatAIJRestoreParCSR_Private(B,&hB);
1042: MatAIJRestoreParCSR_Private(C,&hC);
1043: return(0);
1044: }
1046: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
1047: {
1051: MatSetType(D,MATAIJ);
1052: return(0);
1053: }
1055: /* ---------------------------------------------------- */
1056: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1057: {
1059: C->ops->productnumeric = MatProductNumeric_AB;
1060: return(0);
1061: }
1063: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1064: {
1066: Mat_Product *product = C->product;
1067: PetscBool Ahypre;
1070: PetscObjectTypeCompare((PetscObject)product->A,MATHYPRE,&Ahypre);
1071: if (Ahypre) { /* A is a Hypre matrix */
1072: MatSetType(C,MATHYPRE);
1073: C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1074: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1075: return(0);
1076: }
1077: return(0);
1078: }
1080: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1081: {
1083: C->ops->productnumeric = MatProductNumeric_PtAP;
1084: return(0);
1085: }
1087: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1088: {
1090: Mat_Product *product = C->product;
1091: PetscBool flg;
1092: PetscInt type = 0;
1093: const char *outTypes[4] = {"aij","seqaij","mpiaij","hypre"};
1094: PetscInt ntype = 4;
1095: Mat A = product->A;
1096: PetscBool Ahypre;
1099: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&Ahypre);
1100: if (Ahypre) { /* A is a Hypre matrix */
1101: MatSetType(C,MATHYPRE);
1102: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1103: C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
1104: return(0);
1105: }
1107: /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1108: /* Get runtime option */
1109: if (product->api_user) {
1110: PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatPtAP_HYPRE","Mat");
1111: PetscOptionsEList("-matptap_hypre_outtype","MatPtAP outtype","MatPtAP outtype",outTypes,ntype,outTypes[type],&type,&flg);
1112: PetscOptionsEnd();
1113: } else {
1114: PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_PtAP_HYPRE","Mat");
1115: PetscOptionsEList("-matproduct_ptap_hypre_outtype","MatProduct_PtAP outtype","MatProduct_PtAP",outTypes,ntype,outTypes[type],&type,&flg);
1116: PetscOptionsEnd();
1117: }
1119: if (type == 0 || type == 1 || type == 2) {
1120: MatSetType(C,MATAIJ);
1121: } else if (type == 3) {
1122: MatSetType(C,MATHYPRE);
1123: } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatPtAP outtype is not supported");
1124: C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1125: C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
1126: return(0);
1127: }
1129: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1130: {
1132: Mat_Product *product = C->product;
1135: switch (product->type) {
1136: case MATPRODUCT_AB:
1137: MatProductSetFromOptions_HYPRE_AB(C);
1138: break;
1139: case MATPRODUCT_PtAP:
1140: MatProductSetFromOptions_HYPRE_PtAP(C);
1141: break;
1142: default:
1143: break;
1144: }
1145: return(0);
1146: }
1148: /* -------------------------------------------------------- */
1150: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1151: {
1155: MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);
1156: return(0);
1157: }
1159: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1160: {
1164: MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);
1165: return(0);
1166: }
1168: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1169: {
1173: if (y != z) {
1174: VecCopy(y,z);
1175: }
1176: MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);
1177: return(0);
1178: }
1180: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1181: {
1185: if (y != z) {
1186: VecCopy(y,z);
1187: }
1188: MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);
1189: return(0);
1190: }
1192: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1193: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1194: {
1195: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1196: hypre_ParCSRMatrix *parcsr;
1197: hypre_ParVector *hx,*hy;
1198: PetscErrorCode ierr;
1201: if (trans) {
1202: VecHYPRE_IJVectorPushVecRead(hA->b,x);
1203: if (b != 0.0) { VecHYPRE_IJVectorPushVec(hA->x,y); }
1204: else { VecHYPRE_IJVectorPushVecWrite(hA->x,y); }
1205: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b->ij,(void**)&hx));
1206: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x->ij,(void**)&hy));
1207: } else {
1208: VecHYPRE_IJVectorPushVecRead(hA->x,x);
1209: if (b != 0.0) { VecHYPRE_IJVectorPushVec(hA->b,y); }
1210: else { VecHYPRE_IJVectorPushVecWrite(hA->b,y); }
1211: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x->ij,(void**)&hx));
1212: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b->ij,(void**)&hy));
1213: }
1214: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1215: if (trans) {
1216: PetscStackCallStandard(hypre_ParCSRMatrixMatvecT,(a,parcsr,hx,b,hy));
1217: } else {
1218: PetscStackCallStandard(hypre_ParCSRMatrixMatvec,(a,parcsr,hx,b,hy));
1219: }
1220: VecHYPRE_IJVectorPopVec(hA->x);
1221: VecHYPRE_IJVectorPopVec(hA->b);
1222: return(0);
1223: }
1225: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1226: {
1227: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1231: VecHYPRE_IJVectorDestroy(&hA->x);
1232: VecHYPRE_IJVectorDestroy(&hA->b);
1233: if (hA->ij) {
1234: if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1235: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
1236: }
1237: if (hA->comm) {MPI_Comm_free(&hA->comm);}
1239: MatStashDestroy_Private(&A->stash);
1241: PetscFree(hA->array);
1243: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
1244: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
1245: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_hypre_C",NULL);
1246: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_mpiaij_hypre_C",NULL);
1247: PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
1248: PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
1249: PetscFree(A->data);
1250: return(0);
1251: }
1253: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1254: {
1258: MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
1259: return(0);
1260: }
1262: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1263: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1264: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1265: {
1266: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1267: HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
1268: PetscErrorCode ierr;
1271: A->boundtocpu = bind;
1272: if (hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1273: hypre_ParCSRMatrix *parcsr;
1274: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1275: PetscStackCallStandard(hypre_ParCSRMatrixMigrate,(parcsr, hmem));
1276: }
1277: if (hA->x) { VecHYPRE_IJBindToCPU(hA->x,bind); }
1278: if (hA->b) { VecHYPRE_IJBindToCPU(hA->b,bind); }
1279: return(0);
1280: }
1281: #endif
1283: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1284: {
1285: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1286: PetscMPIInt n;
1287: PetscInt i,j,rstart,ncols,flg;
1288: PetscInt *row,*col;
1289: PetscScalar *val;
1290: PetscErrorCode ierr;
1293: if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");
1295: if (!A->nooffprocentries) {
1296: while (1) {
1297: MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
1298: if (!flg) break;
1300: for (i=0; i<n;) {
1301: /* Now identify the consecutive vals belonging to the same row */
1302: for (j=i,rstart=row[j]; j<n; j++) {
1303: if (row[j] != rstart) break;
1304: }
1305: if (j < n) ncols = j-i;
1306: else ncols = n-i;
1307: /* Now assemble all these values with a single function call */
1308: MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);
1310: i = j;
1311: }
1312: }
1313: MatStashScatterEnd_Private(&A->stash);
1314: }
1316: PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1317: /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1318: /* 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 */
1319: if (!hA->sorted_full) {
1320: hypre_AuxParCSRMatrix *aux_matrix;
1322: /* call destroy just to make sure we do not leak anything */
1323: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1324: PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,(aux_matrix));
1325: hypre_IJMatrixTranslator(hA->ij) = NULL;
1327: /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1328: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1329: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1330: if (aux_matrix) {
1331: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1332: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1333: PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,(aux_matrix));
1334: #else
1335: PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize_v2,(aux_matrix,HYPRE_MEMORY_HOST));
1336: #endif
1337: }
1338: }
1339: {
1340: hypre_ParCSRMatrix *parcsr;
1342: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1343: if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscStackCallStandard(hypre_MatvecCommPkgCreate,(parcsr));
1344: }
1345: if (!hA->x) { VecHYPRE_IJVectorCreate(A->cmap,&hA->x); }
1346: if (!hA->b) { VecHYPRE_IJVectorCreate(A->rmap,&hA->b); }
1347: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1348: MatBindToCPU_HYPRE(A,A->boundtocpu);
1349: #endif
1350: return(0);
1351: }
1353: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1354: {
1355: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1356: PetscErrorCode ierr;
1359: if (!hA->available) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Temporary space is in use");
1361: if (hA->size >= size) {
1362: *array = hA->array;
1363: } else {
1364: PetscFree(hA->array);
1365: hA->size = size;
1366: PetscMalloc(hA->size,&hA->array);
1367: *array = hA->array;
1368: }
1370: hA->available = PETSC_FALSE;
1371: return(0);
1372: }
1374: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1375: {
1376: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1379: *array = NULL;
1380: hA->available = PETSC_TRUE;
1381: return(0);
1382: }
1384: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1385: {
1386: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1387: PetscScalar *vals = (PetscScalar *)v;
1388: HYPRE_Complex *sscr;
1389: PetscInt *cscr[2];
1390: PetscInt i,nzc;
1391: void *array = NULL;
1395: MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(HYPRE_Complex)*nc*nr,&array);
1396: cscr[0] = (PetscInt*)array;
1397: cscr[1] = ((PetscInt*)array)+nc;
1398: sscr = (HYPRE_Complex*)(((PetscInt*)array)+nc*2);
1399: for (i=0,nzc=0;i<nc;i++) {
1400: if (cols[i] >= 0) {
1401: cscr[0][nzc ] = cols[i];
1402: cscr[1][nzc++] = i;
1403: }
1404: }
1405: if (!nzc) {
1406: MatRestoreArray_HYPRE(A,&array);
1407: return(0);
1408: }
1410: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1411: if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1412: hypre_ParCSRMatrix *parcsr;
1414: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1415: PetscStackCallStandard(hypre_ParCSRMatrixMigrate,(parcsr, HYPRE_MEMORY_HOST));
1416: }
1417: #endif
1419: if (ins == ADD_VALUES) {
1420: for (i=0;i<nr;i++) {
1421: if (rows[i] >= 0) {
1422: PetscInt j;
1423: HYPRE_Int hnc = (HYPRE_Int)nzc;
1425: if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1426: for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1427: PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1428: }
1429: vals += nc;
1430: }
1431: } else { /* INSERT_VALUES */
1432: PetscInt rst,ren;
1434: MatGetOwnershipRange(A,&rst,&ren);
1435: for (i=0;i<nr;i++) {
1436: if (rows[i] >= 0) {
1437: PetscInt j;
1438: HYPRE_Int hnc = (HYPRE_Int)nzc;
1440: if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1441: for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1442: /* nonlocal values */
1443: if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],(PetscScalar*)sscr,PETSC_FALSE); }
1444: /* local values */
1445: else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1446: }
1447: vals += nc;
1448: }
1449: }
1451: MatRestoreArray_HYPRE(A,&array);
1452: return(0);
1453: }
1455: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1456: {
1457: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1458: HYPRE_Int *hdnnz,*honnz;
1459: PetscInt i,rs,re,cs,ce,bs;
1460: PetscMPIInt size;
1464: MatGetBlockSize(A,&bs);
1465: PetscLayoutSetUp(A->rmap);
1466: PetscLayoutSetUp(A->cmap);
1467: rs = A->rmap->rstart;
1468: re = A->rmap->rend;
1469: cs = A->cmap->rstart;
1470: ce = A->cmap->rend;
1471: if (!hA->ij) {
1472: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1473: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1474: } else {
1475: HYPRE_BigInt hrs,hre,hcs,hce;
1476: PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1477: if (hre-hrs+1 != re -rs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local rows: IJMatrix [%D,%D), PETSc [%D,%D)",hrs,hre+1,rs,re);
1478: if (hce-hcs+1 != ce -cs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local cols: IJMatrix [%D,%D), PETSc [%D,%D)",hcs,hce+1,cs,ce);
1479: }
1480: if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1481: if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;
1483: if (!dnnz) {
1484: PetscMalloc1(A->rmap->n,&hdnnz);
1485: for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1486: } else {
1487: hdnnz = (HYPRE_Int*)dnnz;
1488: }
1489: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1490: if (size > 1) {
1491: hypre_AuxParCSRMatrix *aux_matrix;
1492: if (!onnz) {
1493: PetscMalloc1(A->rmap->n,&honnz);
1494: for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1495: } else honnz = (HYPRE_Int*)onnz;
1496: /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1497: they assume the user will input the entire row values, properly sorted
1498: In PETSc, we don't make such an assumption and set this flag to 1,
1499: unless the option MAT_SORTED_FULL is set to true.
1500: Also, to avoid possible memory leaks, we destroy and recreate the translator
1501: This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1502: the IJ matrix for us */
1503: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1504: hypre_AuxParCSRMatrixDestroy(aux_matrix);
1505: hypre_IJMatrixTranslator(hA->ij) = NULL;
1506: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1507: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1508: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full;
1509: } else {
1510: honnz = NULL;
1511: PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1512: }
1514: /* reset assembled flag and call the initialize method */
1515: hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1516: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1517: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1518: #else
1519: PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,(hA->ij,HYPRE_MEMORY_HOST));
1520: #endif
1521: if (!dnnz) {
1522: PetscFree(hdnnz);
1523: }
1524: if (!onnz && honnz) {
1525: PetscFree(honnz);
1526: }
1527: /* Match AIJ logic */
1528: A->preallocated = PETSC_TRUE;
1529: A->assembled = PETSC_FALSE;
1530: return(0);
1531: }
1533: /*@C
1534: MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1536: Collective on Mat
1538: Input Parameters:
1539: + A - the matrix
1540: . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix
1541: (same value is used for all local rows)
1542: . dnnz - array containing the number of nonzeros in the various rows of the
1543: DIAGONAL portion of the local submatrix (possibly different for each row)
1544: or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1545: The size of this array is equal to the number of local rows, i.e 'm'.
1546: For matrices that will be factored, you must leave room for (and set)
1547: the diagonal entry even if it is zero.
1548: . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1549: submatrix (same value is used for all local rows).
1550: - onnz - array containing the number of nonzeros in the various rows of the
1551: OFF-DIAGONAL portion of the local submatrix (possibly different for
1552: each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1553: structure. The size of this array is equal to the number
1554: of local rows, i.e 'm'.
1556: Notes:
1557: If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored.
1559: Level: intermediate
1561: .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE
1562: @*/
1563: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1564: {
1570: PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1571: return(0);
1572: }
1574: /*
1575: MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix
1577: Collective
1579: Input Parameters:
1580: + parcsr - the pointer to the hypre_ParCSRMatrix
1581: . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1582: - copymode - PETSc copying options
1584: Output Parameter:
1585: . A - the matrix
1587: Level: intermediate
1589: .seealso: MatHYPRE, PetscCopyMode
1590: */
1591: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1592: {
1593: Mat T;
1594: Mat_HYPRE *hA;
1595: MPI_Comm comm;
1596: PetscInt rstart,rend,cstart,cend,M,N;
1597: PetscBool isseqaij,isseqaijmkl,ismpiaij,isaij,ishyp,isis;
1601: comm = hypre_ParCSRMatrixComm(parcsr);
1602: PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1603: PetscStrcmp(mtype,MATSEQAIJMKL,&isseqaijmkl);
1604: PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1605: PetscStrcmp(mtype,MATAIJ,&isaij);
1606: PetscStrcmp(mtype,MATHYPRE,&ishyp);
1607: PetscStrcmp(mtype,MATIS,&isis);
1608: isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1609: /* TODO */
1610: if (!isaij && !ishyp && !isis) SETERRQ7(comm,PETSC_ERR_SUP,"Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s",mtype,MATAIJ,MATSEQAIJ,MATSEQAIJMKL,MATMPIAIJ,MATIS,MATHYPRE);
1611: /* access ParCSRMatrix */
1612: rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1613: rend = hypre_ParCSRMatrixLastRowIndex(parcsr);
1614: cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1615: cend = hypre_ParCSRMatrixLastColDiag(parcsr);
1616: M = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1617: N = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1619: /* fix for empty local rows/columns */
1620: if (rend < rstart) rend = rstart;
1621: if (cend < cstart) cend = cstart;
1623: /* PETSc convention */
1624: rend++;
1625: cend++;
1626: rend = PetscMin(rend,M);
1627: cend = PetscMin(cend,N);
1629: /* create PETSc matrix with MatHYPRE */
1630: MatCreate(comm,&T);
1631: MatSetSizes(T,rend-rstart,cend-cstart,M,N);
1632: MatSetType(T,MATHYPRE);
1633: hA = (Mat_HYPRE*)(T->data);
1635: /* create HYPRE_IJMatrix */
1636: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
1637: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1639: // TODO DEV
1640: /* create new ParCSR object if needed */
1641: if (ishyp && copymode == PETSC_COPY_VALUES) {
1642: hypre_ParCSRMatrix *new_parcsr;
1643: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
1644: hypre_CSRMatrix *hdiag,*hoffd,*ndiag,*noffd;
1646: new_parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1647: hdiag = hypre_ParCSRMatrixDiag(parcsr);
1648: hoffd = hypre_ParCSRMatrixOffd(parcsr);
1649: ndiag = hypre_ParCSRMatrixDiag(new_parcsr);
1650: noffd = hypre_ParCSRMatrixOffd(new_parcsr);
1651: PetscArraycpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag));
1652: PetscArraycpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd));
1653: #else
1654: new_parcsr = hypre_ParCSRMatrixClone(parcsr,1);
1655: #endif
1656: parcsr = new_parcsr;
1657: copymode = PETSC_OWN_POINTER;
1658: }
1660: /* set ParCSR object */
1661: hypre_IJMatrixObject(hA->ij) = parcsr;
1662: T->preallocated = PETSC_TRUE;
1664: /* set assembled flag */
1665: hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1666: #if 0
1667: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1668: #endif
1669: if (ishyp) {
1670: PetscMPIInt myid = 0;
1672: /* make sure we always have row_starts and col_starts available */
1673: if (HYPRE_AssumedPartitionCheck()) {
1674: MPI_Comm_rank(comm,&myid);
1675: }
1676: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1677: if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1678: PetscLayout map;
1680: MatGetLayouts(T,NULL,&map);
1681: PetscLayoutSetUp(map);
1682: hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1683: }
1684: if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1685: PetscLayout map;
1687: MatGetLayouts(T,&map,NULL);
1688: PetscLayoutSetUp(map);
1689: hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1690: }
1691: #endif
1692: /* prevent from freeing the pointer */
1693: if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1694: *A = T;
1695: MatSetOption(*A,MAT_SORTED_FULL,PETSC_TRUE);
1696: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1697: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1698: } else if (isaij) {
1699: if (copymode != PETSC_OWN_POINTER) {
1700: /* prevent from freeing the pointer */
1701: hA->inner_free = PETSC_FALSE;
1702: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1703: MatDestroy(&T);
1704: } else { /* AIJ return type with PETSC_OWN_POINTER */
1705: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1706: *A = T;
1707: }
1708: } else if (isis) {
1709: MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1710: if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1711: MatDestroy(&T);
1712: }
1713: return(0);
1714: }
1716: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1717: {
1718: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1719: HYPRE_Int type;
1722: if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1723: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1724: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1725: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1726: return(0);
1727: }
1729: /*
1730: MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1732: Not collective
1734: Input Parameters:
1735: + A - the MATHYPRE object
1737: Output Parameter:
1738: . parcsr - the pointer to the hypre_ParCSRMatrix
1740: Level: intermediate
1742: .seealso: MatHYPRE, PetscCopyMode
1743: */
1744: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1745: {
1751: PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1752: return(0);
1753: }
1755: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1756: {
1757: hypre_ParCSRMatrix *parcsr;
1758: hypre_CSRMatrix *ha;
1759: PetscInt rst;
1760: PetscErrorCode ierr;
1763: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1764: MatGetOwnershipRange(A,&rst,NULL);
1765: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1766: if (missing) *missing = PETSC_FALSE;
1767: if (dd) *dd = -1;
1768: ha = hypre_ParCSRMatrixDiag(parcsr);
1769: if (ha) {
1770: PetscInt size,i;
1771: HYPRE_Int *ii,*jj;
1773: size = hypre_CSRMatrixNumRows(ha);
1774: ii = hypre_CSRMatrixI(ha);
1775: jj = hypre_CSRMatrixJ(ha);
1776: for (i = 0; i < size; i++) {
1777: PetscInt j;
1778: PetscBool found = PETSC_FALSE;
1780: for (j = ii[i]; j < ii[i+1] && !found; j++)
1781: found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE;
1783: if (!found) {
1784: PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1785: if (missing) *missing = PETSC_TRUE;
1786: if (dd) *dd = i+rst;
1787: return(0);
1788: }
1789: }
1790: if (!size) {
1791: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1792: if (missing) *missing = PETSC_TRUE;
1793: if (dd) *dd = rst;
1794: }
1795: } else {
1796: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1797: if (missing) *missing = PETSC_TRUE;
1798: if (dd) *dd = rst;
1799: }
1800: return(0);
1801: }
1803: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1804: {
1805: hypre_ParCSRMatrix *parcsr;
1806: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1807: hypre_CSRMatrix *ha;
1808: #endif
1809: PetscErrorCode ierr;
1810: HYPRE_Complex hs;
1813: PetscHYPREScalarCast(s,&hs);
1814: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1815: #if PETSC_PKG_HYPRE_VERSION_GE(2,19,0)
1816: PetscStackCallStandard(hypre_ParCSRMatrixScale,(parcsr,hs));
1817: #else /* diagonal part */
1818: ha = hypre_ParCSRMatrixDiag(parcsr);
1819: if (ha) {
1820: PetscInt size,i;
1821: HYPRE_Int *ii;
1822: HYPRE_Complex *a;
1824: size = hypre_CSRMatrixNumRows(ha);
1825: a = hypre_CSRMatrixData(ha);
1826: ii = hypre_CSRMatrixI(ha);
1827: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1828: }
1829: /* offdiagonal part */
1830: ha = hypre_ParCSRMatrixOffd(parcsr);
1831: if (ha) {
1832: PetscInt size,i;
1833: HYPRE_Int *ii;
1834: HYPRE_Complex *a;
1836: size = hypre_CSRMatrixNumRows(ha);
1837: a = hypre_CSRMatrixData(ha);
1838: ii = hypre_CSRMatrixI(ha);
1839: for (i = 0; i < ii[size]; i++) a[i] *= hs;
1840: }
1841: #endif
1842: return(0);
1843: }
1845: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1846: {
1847: hypre_ParCSRMatrix *parcsr;
1848: HYPRE_Int *lrows;
1849: PetscInt rst,ren,i;
1850: PetscErrorCode ierr;
1853: if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1854: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1855: PetscMalloc1(numRows,&lrows);
1856: MatGetOwnershipRange(A,&rst,&ren);
1857: for (i=0;i<numRows;i++) {
1858: if (rows[i] < rst || rows[i] >= ren)
1859: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1860: lrows[i] = rows[i] - rst;
1861: }
1862: PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1863: PetscFree(lrows);
1864: return(0);
1865: }
1867: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1868: {
1869: PetscErrorCode ierr;
1872: if (ha) {
1873: HYPRE_Int *ii, size;
1874: HYPRE_Complex *a;
1876: size = hypre_CSRMatrixNumRows(ha);
1877: a = hypre_CSRMatrixData(ha);
1878: ii = hypre_CSRMatrixI(ha);
1880: if (a) {PetscArrayzero(a,ii[size]);}
1881: }
1882: return(0);
1883: }
1885: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1886: {
1887: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1890: if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1891: PetscStackCallStandard(HYPRE_IJMatrixSetConstantValues,(hA->ij,0.0));
1892: } else {
1893: hypre_ParCSRMatrix *parcsr;
1894: PetscErrorCode ierr;
1896: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1897: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1898: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1899: }
1900: return(0);
1901: }
1903: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag)
1904: {
1905: PetscInt ii;
1906: HYPRE_Int *i, *j;
1907: HYPRE_Complex *a;
1910: if (!hA) return(0);
1912: i = hypre_CSRMatrixI(hA);
1913: j = hypre_CSRMatrixJ(hA);
1914: a = hypre_CSRMatrixData(hA);
1916: for (ii = 0; ii < N; ii++) {
1917: HYPRE_Int jj, ibeg, iend, irow;
1919: irow = rows[ii];
1920: ibeg = i[irow];
1921: iend = i[irow+1];
1922: for (jj = ibeg; jj < iend; jj++)
1923: if (j[jj] == irow) a[jj] = diag;
1924: else a[jj] = 0.0;
1925: }
1926: return(0);
1927: }
1929: static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1930: {
1931: hypre_ParCSRMatrix *parcsr;
1932: PetscInt *lrows,len;
1933: HYPRE_Complex hdiag;
1934: PetscErrorCode ierr;
1937: if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
1938: PetscHYPREScalarCast(diag,&hdiag);
1939: /* retrieve the internal matrix */
1940: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1941: /* get locally owned rows */
1942: MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1943: /* zero diagonal part */
1944: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag);
1945: /* zero off-diagonal part */
1946: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);
1948: PetscFree(lrows);
1949: return(0);
1950: }
1952: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1953: {
1957: if (mat->nooffprocentries) return(0);
1959: MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1960: return(0);
1961: }
1963: static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1964: {
1965: hypre_ParCSRMatrix *parcsr;
1966: HYPRE_Int hnz;
1967: PetscErrorCode ierr;
1970: /* retrieve the internal matrix */
1971: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1972: /* call HYPRE API */
1973: PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1974: if (nz) *nz = (PetscInt)hnz;
1975: return(0);
1976: }
1978: static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1979: {
1980: hypre_ParCSRMatrix *parcsr;
1981: HYPRE_Int hnz;
1982: PetscErrorCode ierr;
1985: /* retrieve the internal matrix */
1986: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1987: /* call HYPRE API */
1988: hnz = nz ? (HYPRE_Int)(*nz) : 0;
1989: PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1990: return(0);
1991: }
1993: static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1994: {
1995: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1996: PetscInt i;
1999: if (!m || !n) return(0);
2000: /* Ignore negative row indices
2001: * And negative column indices should be automatically ignored in hypre
2002: * */
2003: for (i=0; i<m; i++) {
2004: if (idxm[i] >= 0) {
2005: HYPRE_Int hn = (HYPRE_Int)n;
2006: PetscStackCallStandard(HYPRE_IJMatrixGetValues,(hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n)));
2007: }
2008: }
2009: return(0);
2010: }
2012: static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg)
2013: {
2014: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
2017: switch (op) {
2018: case MAT_NO_OFF_PROC_ENTRIES:
2019: if (flg) {
2020: PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,(hA->ij,0));
2021: }
2022: break;
2023: case MAT_SORTED_FULL:
2024: hA->sorted_full = flg;
2025: break;
2026: default:
2027: break;
2028: }
2029: return(0);
2030: }
2032: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2033: {
2034: PetscErrorCode ierr;
2035: PetscViewerFormat format;
2038: PetscViewerGetFormat(view,&format);
2039: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
2040: if (format != PETSC_VIEWER_NATIVE) {
2041: Mat B;
2042: hypre_ParCSRMatrix *parcsr;
2043: PetscErrorCode (*mview)(Mat,PetscViewer) = NULL;
2045: MatHYPREGetParCSR_HYPRE(A,&parcsr);
2046: MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);
2047: MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);
2048: if (!mview) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation");
2049: (*mview)(B,view);
2050: MatDestroy(&B);
2051: } else {
2052: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
2053: PetscMPIInt size;
2054: PetscBool isascii;
2055: const char *filename;
2057: /* HYPRE uses only text files */
2058: PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
2059: if (!isascii) SETERRQ1(PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name);
2060: PetscViewerFileGetName(view,&filename);
2061: PetscStackCallStandard(HYPRE_IJMatrixPrint,(hA->ij,filename));
2062: MPI_Comm_size(hA->comm,&size);
2063: if (size > 1) {
2064: PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);
2065: } else {
2066: PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);
2067: }
2068: }
2069: return(0);
2070: }
2072: static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B)
2073: {
2074: hypre_ParCSRMatrix *parcsr = NULL;
2075: PetscErrorCode ierr;
2076: PetscCopyMode cpmode;
2079: MatHYPREGetParCSR_HYPRE(A,&parcsr);
2080: if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2081: parcsr = hypre_ParCSRMatrixClone(parcsr,0);
2082: cpmode = PETSC_OWN_POINTER;
2083: } else {
2084: cpmode = PETSC_COPY_VALUES;
2085: }
2086: MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);
2087: return(0);
2088: }
2090: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2091: {
2092: hypre_ParCSRMatrix *acsr,*bcsr;
2093: PetscErrorCode ierr;
2096: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2097: MatHYPREGetParCSR_HYPRE(A,&acsr);
2098: MatHYPREGetParCSR_HYPRE(B,&bcsr);
2099: PetscStackCallStandard(hypre_ParCSRMatrixCopy,(acsr,bcsr,1));
2100: MatSetOption(B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2101: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2102: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2103: } else {
2104: MatCopy_Basic(A,B,str);
2105: }
2106: return(0);
2107: }
2109: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2110: {
2111: hypre_ParCSRMatrix *parcsr;
2112: hypre_CSRMatrix *dmat;
2113: HYPRE_Complex *a;
2114: HYPRE_Complex *data = NULL;
2115: HYPRE_Int *diag = NULL;
2116: PetscInt i;
2117: PetscBool cong;
2118: PetscErrorCode ierr;
2121: MatHasCongruentLayouts(A,&cong);
2122: if (!cong) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns");
2123: if (PetscDefined(USE_DEBUG)) {
2124: PetscBool miss;
2125: MatMissingDiagonal(A,&miss,NULL);
2126: if (miss && A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing");
2127: }
2128: MatHYPREGetParCSR_HYPRE(A,&parcsr);
2129: dmat = hypre_ParCSRMatrixDiag(parcsr);
2130: if (dmat) {
2131: /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
2132: VecGetArray(d,(PetscScalar**)&a);
2133: diag = hypre_CSRMatrixI(dmat);
2134: data = hypre_CSRMatrixData(dmat);
2135: for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]];
2136: VecRestoreArray(d,(PetscScalar**)&a);
2137: }
2138: return(0);
2139: }
2141: #include <petscblaslapack.h>
2143: static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str)
2144: {
2148: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2149: {
2150: Mat B;
2151: hypre_ParCSRMatrix *x,*y,*z;
2153: MatHYPREGetParCSR(Y,&y);
2154: MatHYPREGetParCSR(X,&x);
2155: PetscStackCallStandard(hypre_ParCSRMatrixAdd,(1.0,y,1.0,x,&z));
2156: MatCreateFromParCSR(z,MATHYPRE,PETSC_OWN_POINTER,&B);
2157: MatHeaderMerge(Y,&B);
2158: }
2159: #else
2160: if (str == SAME_NONZERO_PATTERN) {
2161: hypre_ParCSRMatrix *x,*y;
2162: hypre_CSRMatrix *xloc,*yloc;
2163: PetscInt xnnz,ynnz;
2164: HYPRE_Complex *xarr,*yarr;
2165: PetscBLASInt one=1,bnz;
2167: MatHYPREGetParCSR(Y,&y);
2168: MatHYPREGetParCSR(X,&x);
2170: /* diagonal block */
2171: xloc = hypre_ParCSRMatrixDiag(x);
2172: yloc = hypre_ParCSRMatrixDiag(y);
2173: xnnz = 0;
2174: ynnz = 0;
2175: xarr = NULL;
2176: yarr = NULL;
2177: if (xloc) {
2178: xarr = hypre_CSRMatrixData(xloc);
2179: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2180: }
2181: if (yloc) {
2182: yarr = hypre_CSRMatrixData(yloc);
2183: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2184: }
2185: if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %D != %D",xnnz,ynnz);
2186: PetscBLASIntCast(xnnz,&bnz);
2187: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2189: /* off-diagonal block */
2190: xloc = hypre_ParCSRMatrixOffd(x);
2191: yloc = hypre_ParCSRMatrixOffd(y);
2192: xnnz = 0;
2193: ynnz = 0;
2194: xarr = NULL;
2195: yarr = NULL;
2196: if (xloc) {
2197: xarr = hypre_CSRMatrixData(xloc);
2198: xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2199: }
2200: if (yloc) {
2201: yarr = hypre_CSRMatrixData(yloc);
2202: ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2203: }
2204: if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %D != %D",xnnz,ynnz);
2205: PetscBLASIntCast(xnnz,&bnz);
2206: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2207: } else if (str == SUBSET_NONZERO_PATTERN) {
2208: MatAXPY_Basic(Y,a,X,str);
2209: } else {
2210: Mat B;
2212: MatAXPY_Basic_Preallocate(Y,X,&B);
2213: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2214: MatHeaderReplace(Y,&B);
2215: }
2216: #endif
2217: return(0);
2218: }
2220: /*MC
2221: MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2222: based on the hypre IJ interface.
2224: Level: intermediate
2226: .seealso: MatCreate()
2227: M*/
2229: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2230: {
2231: Mat_HYPRE *hB;
2235: PetscNewLog(B,&hB);
2237: hB->inner_free = PETSC_TRUE;
2238: hB->available = PETSC_TRUE;
2239: hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */
2240: hB->size = 0;
2241: hB->array = NULL;
2243: B->data = (void*)hB;
2244: B->rmap->bs = 1;
2245: B->assembled = PETSC_FALSE;
2247: PetscMemzero(B->ops,sizeof(struct _MatOps));
2248: B->ops->mult = MatMult_HYPRE;
2249: B->ops->multtranspose = MatMultTranspose_HYPRE;
2250: B->ops->multadd = MatMultAdd_HYPRE;
2251: B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
2252: B->ops->setup = MatSetUp_HYPRE;
2253: B->ops->destroy = MatDestroy_HYPRE;
2254: B->ops->assemblyend = MatAssemblyEnd_HYPRE;
2255: B->ops->assemblybegin = MatAssemblyBegin_HYPRE;
2256: B->ops->setvalues = MatSetValues_HYPRE;
2257: B->ops->missingdiagonal = MatMissingDiagonal_HYPRE;
2258: B->ops->scale = MatScale_HYPRE;
2259: B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
2260: B->ops->zeroentries = MatZeroEntries_HYPRE;
2261: B->ops->zerorows = MatZeroRows_HYPRE;
2262: B->ops->getrow = MatGetRow_HYPRE;
2263: B->ops->restorerow = MatRestoreRow_HYPRE;
2264: B->ops->getvalues = MatGetValues_HYPRE;
2265: B->ops->setoption = MatSetOption_HYPRE;
2266: B->ops->duplicate = MatDuplicate_HYPRE;
2267: B->ops->copy = MatCopy_HYPRE;
2268: B->ops->view = MatView_HYPRE;
2269: B->ops->getdiagonal = MatGetDiagonal_HYPRE;
2270: B->ops->axpy = MatAXPY_HYPRE;
2271: B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2272: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2273: B->ops->bindtocpu = MatBindToCPU_HYPRE;
2274: B->boundtocpu = PETSC_FALSE;
2275: #endif
2277: /* build cache for off array entries formed */
2278: MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
2280: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
2281: PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
2282: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
2283: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
2284: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_hypre_C",MatProductSetFromOptions_HYPRE);
2285: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_mpiaij_hypre_C",MatProductSetFromOptions_HYPRE);
2286: PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
2287: PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
2288: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2289: #if defined(HYPRE_USING_HIP)
2290: PetscHIPInitializeCheck();
2291: MatSetVecType(B,VECHIP);
2292: #endif
2293: #if defined(HYPRE_USING_CUDA)
2294: PetscCUDAInitializeCheck();
2295: MatSetVecType(B,VECCUDA);
2296: #endif
2297: #endif
2298: return(0);
2299: }
2301: static PetscErrorCode hypre_array_destroy(void *ptr)
2302: {
2304: hypre_TFree(ptr,HYPRE_MEMORY_HOST);
2305: return(0);
2306: }