Actual source code: mhypre.c
petsc-3.10.5 2019-03-28
2: /*
3: Creates hypre ijmatrix from PETSc matrix
4: */
6: #include <petscmathypre.h>
7: #include <petsc/private/matimpl.h>
8: #include <../src/mat/impls/hypre/mhypre.h>
9: #include <../src/mat/impls/aij/mpi/mpiaij.h>
10: #include <../src/vec/vec/impls/hypre/vhyp.h>
11: #include <HYPRE.h>
12: #include <HYPRE_utilities.h>
13: #include <_hypre_parcsr_ls.h>
14: #include <_hypre_sstruct_ls.h>
16: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*);
18: static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
19: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
20: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
21: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
22: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,Vec,Vec,PetscBool);
23: static PetscErrorCode hypre_array_destroy(void*);
24: PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode ins);
26: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
27: {
29: PetscInt i,n_d,n_o;
30: const PetscInt *ia_d,*ia_o;
31: PetscBool done_d=PETSC_FALSE,done_o=PETSC_FALSE;
32: PetscInt *nnz_d=NULL,*nnz_o=NULL;
35: if (A_d) { /* determine number of nonzero entries in local diagonal part */
36: MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);
37: if (done_d) {
38: PetscMalloc1(n_d,&nnz_d);
39: for (i=0; i<n_d; i++) {
40: nnz_d[i] = ia_d[i+1] - ia_d[i];
41: }
42: }
43: MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);
44: }
45: if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
46: MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
47: if (done_o) {
48: PetscMalloc1(n_o,&nnz_o);
49: for (i=0; i<n_o; i++) {
50: nnz_o[i] = ia_o[i+1] - ia_o[i];
51: }
52: }
53: MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
54: }
55: if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */
56: if (!done_o) { /* only diagonal part */
57: PetscMalloc1(n_d,&nnz_o);
58: for (i=0; i<n_d; i++) {
59: nnz_o[i] = 0;
60: }
61: }
62: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
63: PetscFree(nnz_d);
64: PetscFree(nnz_o);
65: }
66: return(0);
67: }
69: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
70: {
72: PetscInt rstart,rend,cstart,cend;
75: PetscLayoutSetUp(A->rmap);
76: PetscLayoutSetUp(A->cmap);
77: rstart = A->rmap->rstart;
78: rend = A->rmap->rend;
79: cstart = A->cmap->rstart;
80: cend = A->cmap->rend;
81: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
82: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
83: {
84: PetscBool same;
85: Mat A_d,A_o;
86: const PetscInt *colmap;
87: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);
88: if (same) {
89: MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);
90: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
91: return(0);
92: }
93: PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);
94: if (same) {
95: MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);
96: MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
97: return(0);
98: }
99: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);
100: if (same) {
101: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
102: return(0);
103: }
104: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);
105: if (same) {
106: MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
107: return(0);
108: }
109: }
110: return(0);
111: }
113: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
114: {
115: PetscErrorCode ierr;
116: PetscInt i,rstart,rend,ncols,nr,nc;
117: const PetscScalar *values;
118: const PetscInt *cols;
119: PetscBool flg;
122: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
123: MatGetSize(A,&nr,&nc);
124: if (flg && nr == nc) {
125: MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
126: return(0);
127: }
128: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
129: if (flg) {
130: MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
131: return(0);
132: }
134: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
135: MatGetOwnershipRange(A,&rstart,&rend);
136: for (i=rstart; i<rend; i++) {
137: MatGetRow(A,i,&ncols,&cols,&values);
138: if (ncols) {
139: PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values));
140: }
141: MatRestoreRow(A,i,&ncols,&cols,&values);
142: }
143: return(0);
144: }
146: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
147: {
148: PetscErrorCode ierr;
149: Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data;
150: HYPRE_Int type;
151: hypre_ParCSRMatrix *par_matrix;
152: hypre_AuxParCSRMatrix *aux_matrix;
153: hypre_CSRMatrix *hdiag;
156: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
157: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
158: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
159: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
160: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
161: /*
162: this is the Hack part where we monkey directly with the hypre datastructures
163: */
164: PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
165: PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
166: PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
168: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
169: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
170: return(0);
171: }
173: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
174: {
175: PetscErrorCode ierr;
176: Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data;
177: Mat_SeqAIJ *pdiag,*poffd;
178: PetscInt i,*garray = pA->garray,*jj,cstart,*pjj;
179: HYPRE_Int type;
180: hypre_ParCSRMatrix *par_matrix;
181: hypre_AuxParCSRMatrix *aux_matrix;
182: hypre_CSRMatrix *hdiag,*hoffd;
185: pdiag = (Mat_SeqAIJ*) pA->A->data;
186: poffd = (Mat_SeqAIJ*) pA->B->data;
187: /* cstart is only valid for square MPIAIJ layed out in the usual way */
188: MatGetOwnershipRange(A,&cstart,NULL);
190: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
191: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
192: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
193: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
194: hdiag = hypre_ParCSRMatrixDiag(par_matrix);
195: hoffd = hypre_ParCSRMatrixOffd(par_matrix);
197: /*
198: this is the Hack part where we monkey directly with the hypre datastructures
199: */
200: PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
201: /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
202: jj = (PetscInt*)hdiag->j;
203: pjj = (PetscInt*)pdiag->j;
204: for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
205: PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
206: PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
207: /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
208: If we hacked a hypre a bit more we might be able to avoid this step */
209: jj = (PetscInt*) hoffd->j;
210: pjj = (PetscInt*) poffd->j;
211: for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
212: PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));
214: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
215: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
216: return(0);
217: }
219: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
220: {
221: Mat_HYPRE* mhA = (Mat_HYPRE*)(A->data);
222: Mat lA;
223: ISLocalToGlobalMapping rl2g,cl2g;
224: IS is;
225: hypre_ParCSRMatrix *hA;
226: hypre_CSRMatrix *hdiag,*hoffd;
227: MPI_Comm comm;
228: PetscScalar *hdd,*hod,*aa,*data;
229: HYPRE_Int *col_map_offd,*hdi,*hdj,*hoi,*hoj;
230: PetscInt *ii,*jj,*iptr,*jptr;
231: PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
232: HYPRE_Int type;
233: PetscErrorCode ierr;
236: comm = PetscObjectComm((PetscObject)A);
237: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
238: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
239: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
240: M = hypre_ParCSRMatrixGlobalNumRows(hA);
241: N = hypre_ParCSRMatrixGlobalNumCols(hA);
242: str = hypre_ParCSRMatrixFirstRowIndex(hA);
243: stc = hypre_ParCSRMatrixFirstColDiag(hA);
244: hdiag = hypre_ParCSRMatrixDiag(hA);
245: hoffd = hypre_ParCSRMatrixOffd(hA);
246: dr = hypre_CSRMatrixNumRows(hdiag);
247: dc = hypre_CSRMatrixNumCols(hdiag);
248: nnz = hypre_CSRMatrixNumNonzeros(hdiag);
249: hdi = hypre_CSRMatrixI(hdiag);
250: hdj = hypre_CSRMatrixJ(hdiag);
251: hdd = hypre_CSRMatrixData(hdiag);
252: oc = hypre_CSRMatrixNumCols(hoffd);
253: nnz += hypre_CSRMatrixNumNonzeros(hoffd);
254: hoi = hypre_CSRMatrixI(hoffd);
255: hoj = hypre_CSRMatrixJ(hoffd);
256: hod = hypre_CSRMatrixData(hoffd);
257: if (reuse != MAT_REUSE_MATRIX) {
258: PetscInt *aux;
260: /* generate l2g maps for rows and cols */
261: ISCreateStride(comm,dr,str,1,&is);
262: ISLocalToGlobalMappingCreateIS(is,&rl2g);
263: ISDestroy(&is);
264: col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
265: PetscMalloc1(dc+oc,&aux);
266: for (i=0; i<dc; i++) aux[i] = i+stc;
267: for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
268: ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
269: ISLocalToGlobalMappingCreateIS(is,&cl2g);
270: ISDestroy(&is);
271: /* create MATIS object */
272: MatCreate(comm,B);
273: MatSetSizes(*B,dr,dc,M,N);
274: MatSetType(*B,MATIS);
275: MatSetLocalToGlobalMapping(*B,rl2g,cl2g);
276: ISLocalToGlobalMappingDestroy(&rl2g);
277: ISLocalToGlobalMappingDestroy(&cl2g);
279: /* allocate CSR for local matrix */
280: PetscMalloc1(dr+1,&iptr);
281: PetscMalloc1(nnz,&jptr);
282: PetscMalloc1(nnz,&data);
283: } else {
284: PetscInt nr;
285: PetscBool done;
286: MatISGetLocalMat(*B,&lA);
287: MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);
288: if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
289: 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);
290: MatSeqAIJGetArray(lA,&data);
291: }
292: /* merge local matrices */
293: ii = iptr;
294: jj = jptr;
295: aa = data;
296: *ii = *(hdi++) + *(hoi++);
297: for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
298: PetscScalar *aold = aa;
299: PetscInt *jold = jj,nc = jd+jo;
300: for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; }
301: for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
302: *(++ii) = *(hdi++) + *(hoi++);
303: PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);
304: }
305: for (; cum<dr; cum++) *(++ii) = nnz;
306: if (reuse != MAT_REUSE_MATRIX) {
307: Mat_SeqAIJ* a;
309: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);
310: MatISSetLocalMat(*B,lA);
311: /* hack SeqAIJ */
312: a = (Mat_SeqAIJ*)(lA->data);
313: a->free_a = PETSC_TRUE;
314: a->free_ij = PETSC_TRUE;
315: MatDestroy(&lA);
316: }
317: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
318: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
319: if (reuse == MAT_INPLACE_MATRIX) {
320: MatHeaderReplace(A,B);
321: }
322: return(0);
323: }
325: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
326: {
327: Mat_HYPRE *hB;
328: MPI_Comm comm = PetscObjectComm((PetscObject)A);
332: if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
333: if (reuse == MAT_REUSE_MATRIX) {
334: /* always destroy the old matrix and create a new memory;
335: hope this does not churn the memory too much. The problem
336: is I do not know if it is possible to put the matrix back to
337: its initial state so that we can directly copy the values
338: the second time through. */
339: hB = (Mat_HYPRE*)((*B)->data);
340: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
341: } else {
342: MatCreate(comm,B);
343: MatSetType(*B,MATHYPRE);
344: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
345: hB = (Mat_HYPRE*)((*B)->data);
346: }
347: MatHYPRE_CreateFromMat(A,hB);
348: MatHYPRE_IJMatrixCopy(A,hB->ij);
349: (*B)->preallocated = PETSC_TRUE;
350: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
351: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
352: return(0);
353: }
355: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
356: {
357: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
358: hypre_ParCSRMatrix *parcsr;
359: hypre_CSRMatrix *hdiag,*hoffd;
360: MPI_Comm comm;
361: PetscScalar *da,*oa,*aptr;
362: PetscInt *dii,*djj,*oii,*ojj,*iptr;
363: PetscInt i,dnnz,onnz,m,n;
364: HYPRE_Int type;
365: PetscMPIInt size;
366: PetscErrorCode ierr;
369: comm = PetscObjectComm((PetscObject)A);
370: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
371: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
372: if (reuse == MAT_REUSE_MATRIX) {
373: PetscBool ismpiaij,isseqaij;
374: PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
375: PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
376: if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
377: }
378: MPI_Comm_size(comm,&size);
380: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
381: hdiag = hypre_ParCSRMatrixDiag(parcsr);
382: hoffd = hypre_ParCSRMatrixOffd(parcsr);
383: m = hypre_CSRMatrixNumRows(hdiag);
384: n = hypre_CSRMatrixNumCols(hdiag);
385: dnnz = hypre_CSRMatrixNumNonzeros(hdiag);
386: onnz = hypre_CSRMatrixNumNonzeros(hoffd);
387: if (reuse == MAT_INITIAL_MATRIX) {
388: PetscMalloc1(m+1,&dii);
389: PetscMalloc1(dnnz,&djj);
390: PetscMalloc1(dnnz,&da);
391: } else if (reuse == MAT_REUSE_MATRIX) {
392: PetscInt nr;
393: PetscBool done;
394: if (size > 1) {
395: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
397: MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
398: 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);
399: 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);
400: MatSeqAIJGetArray(b->A,&da);
401: } else {
402: MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
403: if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
404: 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);
405: MatSeqAIJGetArray(*B,&da);
406: }
407: } else { /* MAT_INPLACE_MATRIX */
408: dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
409: djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
410: da = hypre_CSRMatrixData(hdiag);
411: }
412: PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));
413: PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));
414: PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));
415: iptr = djj;
416: aptr = da;
417: for (i=0; i<m; i++) {
418: PetscInt nc = dii[i+1]-dii[i];
419: PetscSortIntWithScalarArray(nc,iptr,aptr);
420: iptr += nc;
421: aptr += nc;
422: }
423: if (size > 1) {
424: HYPRE_Int *offdj,*coffd;
426: if (reuse == MAT_INITIAL_MATRIX) {
427: PetscMalloc1(m+1,&oii);
428: PetscMalloc1(onnz,&ojj);
429: PetscMalloc1(onnz,&oa);
430: } else if (reuse == MAT_REUSE_MATRIX) {
431: Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
432: PetscInt nr,hr = hypre_CSRMatrixNumRows(hoffd);
433: PetscBool done;
435: MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
436: 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);
437: 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);
438: MatSeqAIJGetArray(b->B,&oa);
439: } else { /* MAT_INPLACE_MATRIX */
440: oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
441: ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
442: oa = hypre_CSRMatrixData(hoffd);
443: }
444: PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));
445: offdj = hypre_CSRMatrixJ(hoffd);
446: coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
447: for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
448: PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));
449: iptr = ojj;
450: aptr = oa;
451: for (i=0; i<m; i++) {
452: PetscInt nc = oii[i+1]-oii[i];
453: PetscSortIntWithScalarArray(nc,iptr,aptr);
454: iptr += nc;
455: aptr += nc;
456: }
457: if (reuse == MAT_INITIAL_MATRIX) {
458: Mat_MPIAIJ *b;
459: Mat_SeqAIJ *d,*o;
461: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
462: /* hack MPIAIJ */
463: b = (Mat_MPIAIJ*)((*B)->data);
464: d = (Mat_SeqAIJ*)b->A->data;
465: o = (Mat_SeqAIJ*)b->B->data;
466: d->free_a = PETSC_TRUE;
467: d->free_ij = PETSC_TRUE;
468: o->free_a = PETSC_TRUE;
469: o->free_ij = PETSC_TRUE;
470: } else if (reuse == MAT_INPLACE_MATRIX) {
471: Mat T;
472: MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
473: hypre_CSRMatrixI(hdiag) = NULL;
474: hypre_CSRMatrixJ(hdiag) = NULL;
475: hypre_CSRMatrixData(hdiag) = NULL;
476: hypre_CSRMatrixI(hoffd) = NULL;
477: hypre_CSRMatrixJ(hoffd) = NULL;
478: hypre_CSRMatrixData(hoffd) = NULL;
479: MatHeaderReplace(A,&T);
480: }
481: } else {
482: oii = NULL;
483: ojj = NULL;
484: oa = NULL;
485: if (reuse == MAT_INITIAL_MATRIX) {
486: Mat_SeqAIJ* b;
487: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
488: /* hack SeqAIJ */
489: b = (Mat_SeqAIJ*)((*B)->data);
490: b->free_a = PETSC_TRUE;
491: b->free_ij = PETSC_TRUE;
492: } else if (reuse == MAT_INPLACE_MATRIX) {
493: Mat T;
494: MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
495: hypre_CSRMatrixI(hdiag) = NULL;
496: hypre_CSRMatrixJ(hdiag) = NULL;
497: hypre_CSRMatrixData(hdiag) = NULL;
498: MatHeaderReplace(A,&T);
499: }
500: }
502: /* we have to use hypre_Tfree to free the arrays */
503: if (reuse == MAT_INPLACE_MATRIX) {
504: void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
505: const char *names[6] = {"_hypre_csr_dii",
506: "_hypre_csr_djj",
507: "_hypre_csr_da",
508: "_hypre_csr_oii",
509: "_hypre_csr_ojj",
510: "_hypre_csr_oa"};
511: for (i=0; i<6; i++) {
512: PetscContainer c;
514: PetscContainerCreate(comm,&c);
515: PetscContainerSetPointer(c,ptrs[i]);
516: PetscContainerSetUserDestroy(c,hypre_array_destroy);
517: PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
518: PetscContainerDestroy(&c);
519: }
520: }
521: return(0);
522: }
524: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
525: {
526: hypre_ParCSRMatrix *tA;
527: hypre_CSRMatrix *hdiag,*hoffd;
528: Mat_SeqAIJ *diag,*offd;
529: PetscInt *garray,noffd,dnnz,onnz,*row_starts,*col_starts;
530: MPI_Comm comm = PetscObjectComm((PetscObject)A);
531: PetscBool ismpiaij,isseqaij;
532: PetscErrorCode ierr;
535: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
536: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
537: if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type);
538: if (ismpiaij) {
539: Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);
541: diag = (Mat_SeqAIJ*)a->A->data;
542: offd = (Mat_SeqAIJ*)a->B->data;
543: garray = a->garray;
544: noffd = a->B->cmap->N;
545: dnnz = diag->nz;
546: onnz = offd->nz;
547: } else {
548: diag = (Mat_SeqAIJ*)A->data;
549: offd = NULL;
550: garray = NULL;
551: noffd = 0;
552: dnnz = diag->nz;
553: onnz = 0;
554: }
556: /* create a temporary ParCSR */
557: if (HYPRE_AssumedPartitionCheck()) {
558: PetscMPIInt myid;
560: MPI_Comm_rank(comm,&myid);
561: row_starts = A->rmap->range + myid;
562: col_starts = A->cmap->range + myid;
563: } else {
564: row_starts = A->rmap->range;
565: col_starts = A->cmap->range;
566: }
567: tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_Int*)row_starts,(HYPRE_Int*)col_starts,noffd,dnnz,onnz);
568: hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
569: hypre_ParCSRMatrixSetColStartsOwner(tA,0);
571: /* set diagonal part */
572: hdiag = hypre_ParCSRMatrixDiag(tA);
573: hypre_CSRMatrixI(hdiag) = (HYPRE_Int*)diag->i;
574: hypre_CSRMatrixJ(hdiag) = (HYPRE_Int*)diag->j;
575: hypre_CSRMatrixData(hdiag) = diag->a;
576: hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
577: hypre_CSRMatrixSetRownnz(hdiag);
578: hypre_CSRMatrixSetDataOwner(hdiag,0);
580: /* set offdiagonal part */
581: hoffd = hypre_ParCSRMatrixOffd(tA);
582: if (offd) {
583: hypre_CSRMatrixI(hoffd) = (HYPRE_Int*)offd->i;
584: hypre_CSRMatrixJ(hoffd) = (HYPRE_Int*)offd->j;
585: hypre_CSRMatrixData(hoffd) = offd->a;
586: hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
587: hypre_CSRMatrixSetRownnz(hoffd);
588: hypre_CSRMatrixSetDataOwner(hoffd,0);
589: hypre_ParCSRMatrixSetNumNonzeros(tA);
590: hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_Int*)garray;
591: }
592: *hA = tA;
593: return(0);
594: }
596: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
597: {
598: hypre_CSRMatrix *hdiag,*hoffd;
601: hdiag = hypre_ParCSRMatrixDiag(*hA);
602: hoffd = hypre_ParCSRMatrixOffd(*hA);
603: /* set pointers to NULL before destroying tA */
604: hypre_CSRMatrixI(hdiag) = NULL;
605: hypre_CSRMatrixJ(hdiag) = NULL;
606: hypre_CSRMatrixData(hdiag) = NULL;
607: hypre_CSRMatrixI(hoffd) = NULL;
608: hypre_CSRMatrixJ(hoffd) = NULL;
609: hypre_CSRMatrixData(hoffd) = NULL;
610: hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
611: hypre_ParCSRMatrixDestroy(*hA);
612: *hA = NULL;
613: return(0);
614: }
616: /* calls RAP from BoomerAMG:
617: the resulting ParCSR will not own the column and row starts
618: It looks like we don't need to have the diagonal entries
619: ordered first in the rows of the diagonal part
620: for boomerAMGBuildCoarseOperator to work */
621: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
622: {
624: HYPRE_Int P_owns_col_starts,R_owns_row_starts;
627: P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
628: R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
629: PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP));
630: PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP));
631: /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
632: hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
633: hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
634: if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
635: if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
636: return(0);
637: }
639: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
640: {
641: Mat B;
642: hypre_ParCSRMatrix *hA,*hP,*hPtAP;
643: PetscErrorCode ierr;
646: MatAIJGetParCSR_Private(A,&hA);
647: MatAIJGetParCSR_Private(P,&hP);
648: MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);
649: MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);
650: MatHeaderMerge(C,&B);
651: MatAIJRestoreParCSR_Private(A,&hA);
652: MatAIJRestoreParCSR_Private(P,&hP);
653: return(0);
654: }
656: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat *C)
657: {
661: MatCreate(PetscObjectComm((PetscObject)A),C);
662: MatSetType(*C,MATAIJ);
663: (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
664: return(0);
665: }
667: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
668: {
669: Mat B;
670: Mat_HYPRE *hP;
671: hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr;
672: HYPRE_Int type;
673: MPI_Comm comm = PetscObjectComm((PetscObject)A);
674: PetscBool ishypre;
675: PetscErrorCode ierr;
678: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
679: if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
680: hP = (Mat_HYPRE*)P->data;
681: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
682: if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
683: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
685: MatAIJGetParCSR_Private(A,&hA);
686: MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
687: MatAIJRestoreParCSR_Private(A,&hA);
689: /* create temporary matrix and merge to C */
690: MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
691: MatHeaderMerge(C,&B);
692: return(0);
693: }
695: static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
696: {
700: if (scall == MAT_INITIAL_MATRIX) {
701: const char *deft = MATAIJ;
702: char type[256];
703: PetscBool flg;
705: PetscObjectOptionsBegin((PetscObject)A);
706: PetscOptionsFList("-matptap_hypre_outtype","Matrix type",NULL,MatList,deft,type,256,&flg);
707: PetscOptionsEnd();
708: PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
709: MatCreate(PetscObjectComm((PetscObject)A),C);
710: if (flg) {
711: MatSetType(*C,type);
712: } else {
713: MatSetType(*C,deft);
714: }
715: (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
716: PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
717: }
718: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
719: (*(*C)->ops->ptapnumeric)(A,P,*C);
720: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
721: return(0);
722: }
724: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
725: {
726: Mat B;
727: hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr;
728: Mat_HYPRE *hA,*hP;
729: PetscBool ishypre;
730: HYPRE_Int type;
731: PetscErrorCode ierr;
734: PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
735: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
736: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
737: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
738: hA = (Mat_HYPRE*)A->data;
739: hP = (Mat_HYPRE*)P->data;
740: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
741: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
742: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
743: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
744: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
745: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
746: MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
747: MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
748: MatHeaderMerge(C,&B);
749: return(0);
750: }
752: static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
753: {
757: if (scall == MAT_INITIAL_MATRIX) {
758: PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
759: MatCreate(PetscObjectComm((PetscObject)A),C);
760: MatSetType(*C,MATHYPRE);
761: (*C)->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
762: PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
763: }
764: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
765: (*(*C)->ops->ptapnumeric)(A,P,*C);
766: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
767: return(0);
768: }
770: /* calls hypre_ParMatmul
771: hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
772: hypre_ParMatrixCreate does not duplicate the communicator
773: It looks like we don't need to have the diagonal entries
774: ordered first in the rows of the diagonal part
775: for boomerAMGBuildCoarseOperator to work */
776: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
777: {
779: PetscStackPush("hypre_ParMatmul");
780: *hAB = hypre_ParMatmul(hA,hB);
781: PetscStackPop;
782: return(0);
783: }
785: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
786: {
787: Mat D;
788: hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
789: PetscErrorCode ierr;
792: MatAIJGetParCSR_Private(A,&hA);
793: MatAIJGetParCSR_Private(B,&hB);
794: MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
795: MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);
796: MatHeaderMerge(C,&D);
797: MatAIJRestoreParCSR_Private(A,&hA);
798: MatAIJRestoreParCSR_Private(B,&hB);
799: return(0);
800: }
802: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat *C)
803: {
807: MatCreate(PetscObjectComm((PetscObject)A),C);
808: MatSetType(*C,MATAIJ);
809: (*C)->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
810: return(0);
811: }
813: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
814: {
815: Mat D;
816: hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
817: Mat_HYPRE *hA,*hB;
818: PetscBool ishypre;
819: HYPRE_Int type;
820: PetscErrorCode ierr;
823: PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
824: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
825: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
826: if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
827: hA = (Mat_HYPRE*)A->data;
828: hB = (Mat_HYPRE*)B->data;
829: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
830: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
831: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
832: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
833: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
834: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
835: MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
836: MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);
837: /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
838: MatHeaderReplace(C,&D);
839: C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
840: return(0);
841: }
843: static PetscErrorCode MatMatMult_HYPRE_HYPRE(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
844: {
848: if (scall == MAT_INITIAL_MATRIX) {
849: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
850: MatCreate(PetscObjectComm((PetscObject)A),C);
851: MatSetType(*C,MATHYPRE);
852: (*C)->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
853: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
854: }
855: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
856: (*(*C)->ops->matmultnumeric)(A,B,*C);
857: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
858: return(0);
859: }
861: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
862: {
863: Mat E;
864: hypre_ParCSRMatrix *hA,*hB,*hC,*hABC;
865: PetscErrorCode ierr;
868: MatAIJGetParCSR_Private(A,&hA);
869: MatAIJGetParCSR_Private(B,&hB);
870: MatAIJGetParCSR_Private(C,&hC);
871: MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
872: MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
873: MatHeaderMerge(D,&E);
874: MatAIJRestoreParCSR_Private(A,&hA);
875: MatAIJRestoreParCSR_Private(B,&hB);
876: MatAIJRestoreParCSR_Private(C,&hC);
877: return(0);
878: }
880: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
881: {
885: MatCreate(PetscObjectComm((PetscObject)A),D);
886: MatSetType(*D,MATAIJ);
887: return(0);
888: }
890: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
891: {
895: MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);
896: return(0);
897: }
899: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
900: {
904: MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);
905: return(0);
906: }
908: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
909: {
910: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
911: hypre_ParCSRMatrix *parcsr;
912: hypre_ParVector *hx,*hy;
913: PetscScalar *ax,*ay,*sax,*say;
914: PetscErrorCode ierr;
917: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
918: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
919: PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
920: VecGetArrayRead(x,(const PetscScalar**)&ax);
921: VecGetArray(y,&ay);
922: if (trans) {
923: VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
924: VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
925: hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
926: VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
927: VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
928: } else {
929: VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
930: VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
931: hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
932: VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
933: VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
934: }
935: VecRestoreArrayRead(x,(const PetscScalar**)&ax);
936: VecRestoreArray(y,&ay);
937: return(0);
938: }
940: static PetscErrorCode MatDestroy_HYPRE(Mat A)
941: {
942: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
946: if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
947: if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
948: if (hA->ij) {
949: if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
950: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
951: }
952: if (hA->comm) { MPI_Comm_free(&hA->comm); }
954: MatStashDestroy_Private(&A->stash);
956: PetscFree(hA->array);
958: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
959: PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
960: PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);
961: PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);
962: PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
963: PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
964: PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_hypre_C",NULL);
965: PetscFree(A->data);
966: return(0);
967: }
969: static PetscErrorCode MatSetUp_HYPRE(Mat A)
970: {
974: MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
975: return(0);
976: }
979: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
980: {
981: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
982: Vec x,b;
983: PetscMPIInt n;
984: PetscInt i,j,rstart,ncols,flg;
985: PetscInt *row,*col;
986: PetscScalar *val;
987: PetscErrorCode ierr;
990: if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");
992: if (!A->nooffprocentries) {
993: while (1) {
994: MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
995: if (!flg) break;
997: for (i=0; i<n; ) {
998: /* Now identify the consecutive vals belonging to the same row */
999: for (j=i,rstart=row[j]; j<n; j++) {
1000: if (row[j] != rstart) break;
1001: }
1002: if (j < n) ncols = j-i;
1003: else ncols = n-i;
1004: /* Now assemble all these values with a single function call */
1005: MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);
1007: i = j;
1008: }
1009: }
1010: MatStashScatterEnd_Private(&A->stash);
1011: }
1013: PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1014: if (hA->x) return(0);
1015: PetscLayoutSetUp(A->rmap);
1016: PetscLayoutSetUp(A->cmap);
1017: VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);
1018: VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);
1019: VecHYPRE_IJVectorCreate(x,&hA->x);
1020: VecHYPRE_IJVectorCreate(b,&hA->b);
1021: VecDestroy(&x);
1022: VecDestroy(&b);
1023: return(0);
1024: }
1026: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1027: {
1028: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1029: PetscErrorCode ierr;
1032: if (!hA->available) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Temporary space is in use");
1034: if (hA->size >= size) *array = hA->array;
1035: else {
1036: PetscFree(hA->array);
1037: hA->size = size;
1038: PetscMalloc(hA->size,&hA->array);
1039: *array = hA->array;
1040: }
1042: hA->available = PETSC_FALSE;
1043: return(0);
1044: }
1046: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1047: {
1048: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1051: *array = NULL;
1052: hA->available = PETSC_TRUE;
1053: return(0);
1054: }
1057: PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1058: {
1059: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1060: PetscScalar *vals = (PetscScalar *)v;
1061: PetscScalar *sscr;
1062: PetscInt *cscr[2];
1063: PetscInt i,nzc;
1064: void *array = NULL;
1065: PetscErrorCode ierr;
1068: MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(PetscScalar)*nc*nr,&array);
1069: cscr[0] = (PetscInt*)array;
1070: cscr[1] = ((PetscInt*)array)+nc;
1071: sscr = (PetscScalar*)(((PetscInt*)array)+nc*2);
1072: for (i=0,nzc=0;i<nc;i++) {
1073: if (cols[i] >= 0) {
1074: cscr[0][nzc ] = cols[i];
1075: cscr[1][nzc++] = i;
1076: }
1077: }
1078: if (!nzc) {
1079: MatRestoreArray_HYPRE(A,&array);
1080: return(0);
1081: }
1083: if (ins == ADD_VALUES) {
1084: for (i=0;i<nr;i++) {
1085: if (rows[i] >= 0 && nzc) {
1086: PetscInt j;
1087: for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1088: PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1089: }
1090: vals += nc;
1091: }
1092: } else { /* INSERT_VALUES */
1094: PetscInt rst,ren;
1095: MatGetOwnershipRange(A,&rst,&ren);
1097: for (i=0;i<nr;i++) {
1098: if (rows[i] >= 0 && nzc) {
1099: PetscInt j;
1100: for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1101: /* nonlocal values */
1102: if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],sscr, PETSC_FALSE); }
1103: /* local values */
1104: else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1105: }
1106: vals += nc;
1107: }
1108: }
1110: MatRestoreArray_HYPRE(A,&array);
1111: return(0);
1112: }
1114: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1115: {
1116: Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1117: HYPRE_Int *hdnnz,*honnz;
1118: PetscInt i,rs,re,cs,ce,bs;
1119: PetscMPIInt size;
1120: PetscErrorCode ierr;
1123: MatGetBlockSize(A,&bs);
1124: PetscLayoutSetUp(A->rmap);
1125: PetscLayoutSetUp(A->cmap);
1126: rs = A->rmap->rstart;
1127: re = A->rmap->rend;
1128: cs = A->cmap->rstart;
1129: ce = A->cmap->rend;
1130: if (!hA->ij) {
1131: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1132: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1133: } else {
1134: HYPRE_Int hrs,hre,hcs,hce;
1135: PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1136: 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);
1137: 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);
1138: }
1139: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1141: if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1142: if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;
1144: if (!dnnz) {
1145: PetscMalloc1(A->rmap->n,&hdnnz);
1146: for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1147: } else {
1148: hdnnz = (HYPRE_Int*)dnnz;
1149: }
1150: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1151: if (size > 1) {
1152: if (!onnz) {
1153: PetscMalloc1(A->rmap->n,&honnz);
1154: for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1155: } else {
1156: honnz = (HYPRE_Int*)onnz;
1157: }
1158: PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1159: } else {
1160: honnz = NULL;
1161: PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1162: }
1163: if (!dnnz) {
1164: PetscFree(hdnnz);
1165: }
1166: if (!onnz && honnz) {
1167: PetscFree(honnz);
1168: }
1169: A->preallocated = PETSC_TRUE;
1171: /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0 */
1172: {
1173: hypre_AuxParCSRMatrix *aux_matrix;
1174: aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1175: hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
1176: }
1177: return(0);
1178: }
1180: /*@C
1181: MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1183: Collective on Mat
1185: Input Parameters:
1186: + A - the matrix
1187: . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix
1188: (same value is used for all local rows)
1189: . dnnz - array containing the number of nonzeros in the various rows of the
1190: DIAGONAL portion of the local submatrix (possibly different for each row)
1191: or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1192: The size of this array is equal to the number of local rows, i.e 'm'.
1193: For matrices that will be factored, you must leave room for (and set)
1194: the diagonal entry even if it is zero.
1195: . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1196: submatrix (same value is used for all local rows).
1197: - onnz - array containing the number of nonzeros in the various rows of the
1198: OFF-DIAGONAL portion of the local submatrix (possibly different for
1199: each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1200: structure. The size of this array is equal to the number
1201: of local rows, i.e 'm'.
1203: Notes:
1204: If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored.
1206: Level: intermediate
1208: .keywords: matrix, aij, compressed row, sparse, parallel
1210: .seealso: MatCreate(), MatMPIAIJSetPreallocation, MATHYPRE
1211: @*/
1212: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1213: {
1219: PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1220: return(0);
1221: }
1223: /*
1224: MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix
1226: Collective
1228: Input Parameters:
1229: + vparcsr - the pointer to the hypre_ParCSRMatrix
1230: . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1231: - copymode - PETSc copying options
1233: Output Parameter:
1234: . A - the matrix
1236: Level: intermediate
1238: .seealso: MatHYPRE, PetscCopyMode
1239: */
1240: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *vparcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1241: {
1242: Mat T;
1243: Mat_HYPRE *hA;
1244: hypre_ParCSRMatrix *parcsr;
1245: MPI_Comm comm;
1246: PetscInt rstart,rend,cstart,cend,M,N;
1247: PetscBool isseqaij,ismpiaij,isaij,ishyp,isis;
1248: PetscErrorCode ierr;
1251: parcsr = (hypre_ParCSRMatrix *)vparcsr;
1252: comm = hypre_ParCSRMatrixComm(parcsr);
1253: PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1254: PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1255: PetscStrcmp(mtype,MATAIJ,&isaij);
1256: PetscStrcmp(mtype,MATHYPRE,&ishyp);
1257: PetscStrcmp(mtype,MATIS,&isis);
1258: isaij = (PetscBool)(isseqaij || ismpiaij || isaij);
1259: if (!isaij && !ishyp && !isis) SETERRQ6(comm,PETSC_ERR_SUP,"Unsupported MatType %s! Supported types are %s, %s, %s, %s, and %s",mtype,MATAIJ,MATSEQAIJ,MATMPIAIJ,MATIS,MATHYPRE);
1260: if (ishyp && copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");
1262: /* access ParCSRMatrix */
1263: rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1264: rend = hypre_ParCSRMatrixLastRowIndex(parcsr);
1265: cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1266: cend = hypre_ParCSRMatrixLastColDiag(parcsr);
1267: M = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1268: N = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1270: /* fix for empty local rows/columns */
1271: if (rend < rstart) rend = rstart;
1272: if (cend < cstart) cend = cstart;
1274: /* create PETSc matrix with MatHYPRE */
1275: MatCreate(comm,&T);
1276: MatSetSizes(T,rend-rstart+1,cend-cstart+1,M,N);
1277: MatSetType(T,MATHYPRE);
1278: hA = (Mat_HYPRE*)(T->data);
1280: /* create HYPRE_IJMatrix */
1281: PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
1283: /* set ParCSR object */
1284: PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1285: hypre_IJMatrixObject(hA->ij) = parcsr;
1286: T->preallocated = PETSC_TRUE;
1288: /* set assembled flag */
1289: hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1290: PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1291: if (ishyp) {
1292: PetscMPIInt myid = 0;
1294: /* make sure we always have row_starts and col_starts available */
1295: if (HYPRE_AssumedPartitionCheck()) {
1296: MPI_Comm_rank(comm,&myid);
1297: }
1298: if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1299: PetscLayout map;
1301: MatGetLayouts(T,NULL,&map);
1302: PetscLayoutSetUp(map);
1303: hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1304: }
1305: if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1306: PetscLayout map;
1308: MatGetLayouts(T,&map,NULL);
1309: PetscLayoutSetUp(map);
1310: hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1311: }
1312: /* prevent from freeing the pointer */
1313: if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1314: *A = T;
1315: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1316: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1317: } else if (isaij) {
1318: if (copymode != PETSC_OWN_POINTER) {
1319: /* prevent from freeing the pointer */
1320: hA->inner_free = PETSC_FALSE;
1321: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1322: MatDestroy(&T);
1323: } else { /* AIJ return type with PETSC_OWN_POINTER */
1324: MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1325: *A = T;
1326: }
1327: } else if (isis) {
1328: MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1329: if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1330: MatDestroy(&T);
1331: }
1332: return(0);
1333: }
1335: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1336: {
1337: Mat_HYPRE* hA = (Mat_HYPRE*)A->data;
1338: HYPRE_Int type;
1339: PetscErrorCode ierr;
1342: if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1343: PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1344: if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1345: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1346: return(0);
1347: }
1349: /*
1350: MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1352: Not collective
1354: Input Parameters:
1355: + A - the MATHYPRE object
1357: Output Parameter:
1358: . parcsr - the pointer to the hypre_ParCSRMatrix
1360: Level: intermediate
1362: .seealso: MatHYPRE, PetscCopyMode
1363: */
1364: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1365: {
1371: PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1372: return(0);
1373: }
1375: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1376: {
1377: hypre_ParCSRMatrix *parcsr;
1378: hypre_CSRMatrix *ha;
1379: PetscInt rst;
1380: PetscErrorCode ierr;
1383: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1384: MatGetOwnershipRange(A,&rst,NULL);
1385: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1386: if (missing) *missing = PETSC_FALSE;
1387: if (dd) *dd = -1;
1388: ha = hypre_ParCSRMatrixDiag(parcsr);
1389: if (ha) {
1390: PetscInt size,i;
1391: HYPRE_Int *ii,*jj;
1393: size = hypre_CSRMatrixNumRows(ha);
1394: ii = hypre_CSRMatrixI(ha);
1395: jj = hypre_CSRMatrixJ(ha);
1396: for (i = 0; i < size; i++) {
1397: PetscInt j;
1398: PetscBool found = PETSC_FALSE;
1400: for (j = ii[i]; j < ii[i+1] && !found; j++)
1401: found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE;
1403: if (!found) {
1404: PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1405: if (missing) *missing = PETSC_TRUE;
1406: if (dd) *dd = i+rst;
1407: return(0);
1408: }
1409: }
1410: if (!size) {
1411: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1412: if (missing) *missing = PETSC_TRUE;
1413: if (dd) *dd = rst;
1414: }
1415: } else {
1416: PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1417: if (missing) *missing = PETSC_TRUE;
1418: if (dd) *dd = rst;
1419: }
1420: return(0);
1421: }
1423: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1424: {
1425: hypre_ParCSRMatrix *parcsr;
1426: hypre_CSRMatrix *ha;
1427: PetscErrorCode ierr;
1430: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1431: /* diagonal part */
1432: ha = hypre_ParCSRMatrixDiag(parcsr);
1433: if (ha) {
1434: PetscInt size,i;
1435: HYPRE_Int *ii;
1436: PetscScalar *a;
1438: size = hypre_CSRMatrixNumRows(ha);
1439: a = hypre_CSRMatrixData(ha);
1440: ii = hypre_CSRMatrixI(ha);
1441: for (i = 0; i < ii[size]; i++) a[i] *= s;
1442: }
1443: /* offdiagonal part */
1444: ha = hypre_ParCSRMatrixOffd(parcsr);
1445: if (ha) {
1446: PetscInt size,i;
1447: HYPRE_Int *ii;
1448: PetscScalar *a;
1450: size = hypre_CSRMatrixNumRows(ha);
1451: a = hypre_CSRMatrixData(ha);
1452: ii = hypre_CSRMatrixI(ha);
1453: for (i = 0; i < ii[size]; i++) a[i] *= s;
1454: }
1455: return(0);
1456: }
1458: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1459: {
1460: hypre_ParCSRMatrix *parcsr;
1461: HYPRE_Int *lrows;
1462: PetscInt rst,ren,i;
1463: PetscErrorCode ierr;
1466: if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1467: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1468: PetscMalloc1(numRows,&lrows);
1469: MatGetOwnershipRange(A,&rst,&ren);
1470: for (i=0;i<numRows;i++) {
1471: if (rows[i] < rst || rows[i] >= ren)
1472: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1473: lrows[i] = rows[i] - rst;
1474: }
1475: PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1476: PetscFree(lrows);
1477: return(0);
1478: }
1480: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1481: {
1482: PetscErrorCode ierr;
1485: if (ha) {
1486: HYPRE_Int *ii, size;
1487: HYPRE_Complex *a;
1489: size = hypre_CSRMatrixNumRows(ha);
1490: a = hypre_CSRMatrixData(ha);
1491: ii = hypre_CSRMatrixI(ha);
1493: if (a) { PetscMemzero(a,(ii[size])*sizeof(HYPRE_Complex)); }
1494: }
1495: return(0);
1496: }
1498: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1499: {
1500: hypre_ParCSRMatrix *parcsr;
1501: PetscErrorCode ierr;
1504: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1505: /* diagonal part */
1506: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1507: /* off-diagonal part */
1508: MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1509: return(0);
1510: }
1513: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],PetscScalar diag)
1514: {
1515: PetscInt ii, jj, ibeg, iend, irow;
1516: PetscInt *i, *j;
1517: PetscScalar *a;
1521: if (!hA) return(0);
1523: i = (PetscInt*) hypre_CSRMatrixI(hA);
1524: j = (PetscInt*) hypre_CSRMatrixJ(hA);
1525: a = hypre_CSRMatrixData(hA);
1527: for (ii = 0; ii < N; ii++) {
1528: irow = rows[ii];
1529: ibeg = i[irow];
1530: iend = i[irow+1];
1531: for (jj = ibeg; jj < iend; jj++)
1532: if (j[jj] == irow) a[jj] = diag;
1533: else a[jj] = 0.0;
1534: }
1536: return(0);
1537: }
1539: PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1540: {
1541: hypre_ParCSRMatrix *parcsr;
1542: PetscInt *lrows,len;
1543: PetscErrorCode ierr;
1546: if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size \n");
1547: /* retrieve the internal matrix */
1548: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1549: /* get locally owned rows */
1550: MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1551: /* zero diagonal part */
1552: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,diag);
1553: /* zero off-diagonal part */
1554: MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);
1556: PetscFree(lrows);
1557: return(0);
1558: }
1561: PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1562: {
1566: if (mat->nooffprocentries) return(0);
1568: MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1569: return(0);
1570: }
1572: PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1573: {
1574: hypre_ParCSRMatrix *parcsr;
1575: PetscErrorCode ierr;
1578: /* retrieve the internal matrix */
1579: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1580: /* call HYPRE API */
1581: PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1582: return(0);
1583: }
1585: PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1586: {
1587: hypre_ParCSRMatrix *parcsr;
1588: PetscErrorCode ierr;
1591: /* retrieve the internal matrix */
1592: MatHYPREGetParCSR_HYPRE(A,&parcsr);
1593: /* call HYPRE API */
1594: PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1595: return(0);
1596: }
1600: PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1601: {
1602: HYPRE_IJMatrix *hIJ = (HYPRE_IJMatrix*)A->data;
1603: PetscErrorCode ierr;
1604: PetscInt i;
1606: if (!m || !n) return(0);
1608: /* Ignore negative row indices
1609: * And negative column indices should be automatically ignored in hypre
1610: * */
1611: for (i=0; i<m; i++)
1612: if (idxm[i] >= 0) PetscStackCallStandard(HYPRE_IJMatrixGetValues,(*hIJ,1,(HYPRE_Int*)&n,(HYPRE_Int*)&idxm[i],(HYPRE_Int*)idxn,&v[i*n]));
1614: return(0);
1615: }
1618: /*MC
1619: MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
1620: based on the hypre IJ interface.
1622: Level: intermediate
1624: .seealso: MatCreate()
1625: M*/
1627: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
1628: {
1629: Mat_HYPRE *hB;
1633: PetscNewLog(B,&hB);
1634: hB->inner_free = PETSC_TRUE;
1635: hB->available = PETSC_TRUE;
1636: hB->size = 0;
1637: hB->array = NULL;
1639: B->data = (void*)hB;
1640: B->rmap->bs = 1;
1641: B->assembled = PETSC_FALSE;
1643: PetscMemzero(B->ops,sizeof(struct _MatOps));
1644: B->ops->mult = MatMult_HYPRE;
1645: B->ops->multtranspose = MatMultTranspose_HYPRE;
1646: B->ops->setup = MatSetUp_HYPRE;
1647: B->ops->destroy = MatDestroy_HYPRE;
1649: /* build cache for off array entries formed */
1650: MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
1651: B->ops->assemblyend = MatAssemblyEnd_HYPRE;
1652: B->ops->assemblybegin = MatAssemblyBegin_HYPRE;
1654: B->ops->ptap = MatPtAP_HYPRE_HYPRE;
1655: B->ops->matmult = MatMatMult_HYPRE_HYPRE;
1656: B->ops->setvalues = MatSetValues_HYPRE;
1657: B->ops->missingdiagonal = MatMissingDiagonal_HYPRE;
1658: B->ops->scale = MatScale_HYPRE;
1659: B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
1660: B->ops->zeroentries = MatZeroEntries_HYPRE;
1661: B->ops->zerorows = MatZeroRows_HYPRE;
1662: B->ops->getrow = MatGetRow_HYPRE;
1663: B->ops->restorerow = MatRestoreRow_HYPRE;
1664: B->ops->getvalues = MatGetValues_HYPRE;
1666: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
1667: PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
1668: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
1669: PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
1670: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);
1671: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);
1672: PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
1673: PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
1674: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_hypre_C",MatPtAP_IS_XAIJ);
1675: return(0);
1676: }
1678: static PetscErrorCode hypre_array_destroy(void *ptr)
1679: {
1681: hypre_TFree(ptr,HYPRE_MEMORY_HOST);
1682: return(0);
1683: }