xref: /petsc/src/mat/impls/hypre/mhypre.c (revision fa2f3a98544cffa8147e57031dd6780dd55215d2) 
1 
2 /*
3     Creates hypre ijmatrix from PETSc matrix
4 */
5 
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>
15 
16 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*);
17 
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 
25 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
26 {
27   PetscErrorCode ierr;
28   PetscInt       i,n_d,n_o;
29   const PetscInt *ia_d,*ia_o;
30   PetscBool      done_d=PETSC_FALSE,done_o=PETSC_FALSE;
31   PetscInt       *nnz_d=NULL,*nnz_o=NULL;
32 
33   PetscFunctionBegin;
34   if (A_d) { /* determine number of nonzero entries in local diagonal part */
35     ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);CHKERRQ(ierr);
36     if (done_d) {
37       ierr = PetscMalloc1(n_d,&nnz_d);CHKERRQ(ierr);
38       for (i=0; i<n_d; i++) {
39         nnz_d[i] = ia_d[i+1] - ia_d[i];
40       }
41     }
42     ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);CHKERRQ(ierr);
43   }
44   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
45     ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
46     if (done_o) {
47       ierr = PetscMalloc1(n_o,&nnz_o);CHKERRQ(ierr);
48       for (i=0; i<n_o; i++) {
49         nnz_o[i] = ia_o[i+1] - ia_o[i];
50       }
51     }
52     ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
53   }
54   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
55     if (!done_o) { /* only diagonal part */
56       ierr = PetscMalloc1(n_d,&nnz_o);CHKERRQ(ierr);
57       for (i=0; i<n_d; i++) {
58         nnz_o[i] = 0;
59       }
60     }
61     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
62     ierr = PetscFree(nnz_d);CHKERRQ(ierr);
63     ierr = PetscFree(nnz_o);CHKERRQ(ierr);
64   }
65   PetscFunctionReturn(0);
66 }
67 
68 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
69 {
70   PetscErrorCode ierr;
71   PetscInt       rstart,rend,cstart,cend;
72 
73   PetscFunctionBegin;
74   ierr   = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
75   ierr   = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
76   rstart = A->rmap->rstart;
77   rend   = A->rmap->rend;
78   cstart = A->cmap->rstart;
79   cend   = A->cmap->rend;
80   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
81   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
82   {
83     PetscBool      same;
84     Mat            A_d,A_o;
85     const PetscInt *colmap;
86     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr);
87     if (same) {
88       ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
89       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
90       PetscFunctionReturn(0);
91     }
92     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr);
93     if (same) {
94       ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
95       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
96       PetscFunctionReturn(0);
97     }
98     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr);
99     if (same) {
100       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
101       PetscFunctionReturn(0);
102     }
103     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr);
104     if (same) {
105       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
106       PetscFunctionReturn(0);
107     }
108   }
109   PetscFunctionReturn(0);
110 }
111 
112 static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
113 {
114   PetscErrorCode    ierr;
115   PetscInt          i,rstart,rend,ncols,nr,nc;
116   const PetscScalar *values;
117   const PetscInt    *cols;
118   PetscBool         flg;
119 
120   PetscFunctionBegin;
121   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
122   ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr);
123   if (flg && nr == nc) {
124     ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr);
125     PetscFunctionReturn(0);
126   }
127   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
128   if (flg) {
129     ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr);
130     PetscFunctionReturn(0);
131   }
132 
133   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
134   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
135   for (i=rstart; i<rend; i++) {
136     ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
137     if (ncols) {
138       PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values));
139     }
140     ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
141   }
142   PetscFunctionReturn(0);
143 }
144 
145 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
146 {
147   PetscErrorCode        ierr;
148   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ*)A->data;
149   HYPRE_Int             type;
150   hypre_ParCSRMatrix    *par_matrix;
151   hypre_AuxParCSRMatrix *aux_matrix;
152   hypre_CSRMatrix       *hdiag;
153 
154   PetscFunctionBegin;
155   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
156   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
157   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
158   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
159   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
160   /*
161        this is the Hack part where we monkey directly with the hypre datastructures
162   */
163   ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
164   ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
165   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
166 
167   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
168   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
169   PetscFunctionReturn(0);
170 }
171 
172 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
173 {
174   PetscErrorCode        ierr;
175   Mat_MPIAIJ            *pA = (Mat_MPIAIJ*)A->data;
176   Mat_SeqAIJ            *pdiag,*poffd;
177   PetscInt              i,*garray = pA->garray,*jj,cstart,*pjj;
178   HYPRE_Int             type;
179   hypre_ParCSRMatrix    *par_matrix;
180   hypre_AuxParCSRMatrix *aux_matrix;
181   hypre_CSRMatrix       *hdiag,*hoffd;
182 
183   PetscFunctionBegin;
184   pdiag = (Mat_SeqAIJ*) pA->A->data;
185   poffd = (Mat_SeqAIJ*) pA->B->data;
186   /* cstart is only valid for square MPIAIJ layed out in the usual way */
187   ierr = MatGetOwnershipRange(A,&cstart,NULL);CHKERRQ(ierr);
188 
189   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
190   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
191   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
192   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
193   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
194   hoffd = hypre_ParCSRMatrixOffd(par_matrix);
195 
196   /*
197        this is the Hack part where we monkey directly with the hypre datastructures
198   */
199   ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
200   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
201   jj  = (PetscInt*)hdiag->j;
202   pjj = (PetscInt*)pdiag->j;
203   for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
204   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
205   ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
206   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
207      If we hacked a hypre a bit more we might be able to avoid this step */
208   jj  = (PetscInt*) hoffd->j;
209   pjj = (PetscInt*) poffd->j;
210   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
211   ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));
212 
213   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
214   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
215   PetscFunctionReturn(0);
216 }
217 
218 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
219 {
220   Mat_HYPRE*             mhA = (Mat_HYPRE*)(A->data);
221   Mat                    lA;
222   ISLocalToGlobalMapping rl2g,cl2g;
223   IS                     is;
224   hypre_ParCSRMatrix     *hA;
225   hypre_CSRMatrix        *hdiag,*hoffd;
226   MPI_Comm               comm;
227   PetscScalar            *hdd,*hod,*aa,*data;
228   HYPRE_Int              *col_map_offd,*hdi,*hdj,*hoi,*hoj;
229   PetscInt               *ii,*jj,*iptr,*jptr;
230   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
231   HYPRE_Int              type;
232   PetscErrorCode         ierr;
233 
234   PetscFunctionBegin;
235   comm = PetscObjectComm((PetscObject)A);
236   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
237   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
238   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
239   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
240   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
241   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
242   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
243   hdiag = hypre_ParCSRMatrixDiag(hA);
244   hoffd = hypre_ParCSRMatrixOffd(hA);
245   dr    = hypre_CSRMatrixNumRows(hdiag);
246   dc    = hypre_CSRMatrixNumCols(hdiag);
247   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
248   hdi   = hypre_CSRMatrixI(hdiag);
249   hdj   = hypre_CSRMatrixJ(hdiag);
250   hdd   = hypre_CSRMatrixData(hdiag);
251   oc    = hypre_CSRMatrixNumCols(hoffd);
252   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
253   hoi   = hypre_CSRMatrixI(hoffd);
254   hoj   = hypre_CSRMatrixJ(hoffd);
255   hod   = hypre_CSRMatrixData(hoffd);
256   if (reuse != MAT_REUSE_MATRIX) {
257     PetscInt *aux;
258 
259     /* generate l2g maps for rows and cols */
260     ierr = ISCreateStride(comm,dr,str,1,&is);CHKERRQ(ierr);
261     ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr);
262     ierr = ISDestroy(&is);CHKERRQ(ierr);
263     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
264     ierr = PetscMalloc1(dc+oc,&aux);CHKERRQ(ierr);
265     for (i=0; i<dc; i++) aux[i] = i+stc;
266     for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
267     ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);CHKERRQ(ierr);
268     ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr);
269     ierr = ISDestroy(&is);CHKERRQ(ierr);
270     /* create MATIS object */
271     ierr = MatCreate(comm,B);CHKERRQ(ierr);
272     ierr = MatSetSizes(*B,dr,dc,M,N);CHKERRQ(ierr);
273     ierr = MatSetType(*B,MATIS);CHKERRQ(ierr);
274     ierr = MatSetLocalToGlobalMapping(*B,rl2g,cl2g);CHKERRQ(ierr);
275     ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr);
276     ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr);
277 
278     /* allocate CSR for local matrix */
279     ierr = PetscMalloc1(dr+1,&iptr);CHKERRQ(ierr);
280     ierr = PetscMalloc1(nnz,&jptr);CHKERRQ(ierr);
281     ierr = PetscMalloc1(nnz,&data);CHKERRQ(ierr);
282   } else {
283     PetscInt  nr;
284     PetscBool done;
285     ierr = MatISGetLocalMat(*B,&lA);CHKERRQ(ierr);
286     ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);CHKERRQ(ierr);
287     if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
288     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);
289     ierr = MatSeqAIJGetArray(lA,&data);CHKERRQ(ierr);
290   }
291   /* merge local matrices */
292   ii   = iptr;
293   jj   = jptr;
294   aa   = data;
295   *ii  = *(hdi++) + *(hoi++);
296   for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
297     PetscScalar *aold = aa;
298     PetscInt    *jold = jj,nc = jd+jo;
299     for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
300     for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
301     *(++ii) = *(hdi++) + *(hoi++);
302     ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);CHKERRQ(ierr);
303   }
304   for (; cum<dr; cum++) *(++ii) = nnz;
305   if (reuse != MAT_REUSE_MATRIX) {
306     Mat_SeqAIJ* a;
307 
308     ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr);
309     ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr);
310     /* hack SeqAIJ */
311     a          = (Mat_SeqAIJ*)(lA->data);
312     a->free_a  = PETSC_TRUE;
313     a->free_ij = PETSC_TRUE;
314     ierr = MatDestroy(&lA);CHKERRQ(ierr);
315   }
316   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
317   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
318   if (reuse == MAT_INPLACE_MATRIX) {
319     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
320   }
321   PetscFunctionReturn(0);
322 }
323 
324 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
325 {
326   Mat_HYPRE      *hB;
327   MPI_Comm       comm = PetscObjectComm((PetscObject)A);
328   PetscErrorCode ierr;
329 
330   PetscFunctionBegin;
331   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
332   if (reuse == MAT_REUSE_MATRIX) {
333     /* always destroy the old matrix and create a new memory;
334        hope this does not churn the memory too much. The problem
335        is I do not know if it is possible to put the matrix back to
336        its initial state so that we can directly copy the values
337        the second time through. */
338     hB = (Mat_HYPRE*)((*B)->data);
339     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
340   } else {
341     ierr = MatCreate(comm,B);CHKERRQ(ierr);
342     ierr = MatSetType(*B,MATHYPRE);CHKERRQ(ierr);
343     ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
344     hB   = (Mat_HYPRE*)((*B)->data);
345   }
346   ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr);
347   ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr);
348   (*B)->preallocated = PETSC_TRUE;
349   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
350   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
351   PetscFunctionReturn(0);
352 }
353 
354 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
355 {
356   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
357   hypre_ParCSRMatrix *parcsr;
358   hypre_CSRMatrix    *hdiag,*hoffd;
359   MPI_Comm           comm;
360   PetscScalar        *da,*oa,*aptr;
361   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
362   PetscInt           i,dnnz,onnz,m,n;
363   HYPRE_Int          type;
364   PetscMPIInt        size;
365   PetscErrorCode     ierr;
366 
367   PetscFunctionBegin;
368   comm = PetscObjectComm((PetscObject)A);
369   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
370   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
371   if (reuse == MAT_REUSE_MATRIX) {
372     PetscBool ismpiaij,isseqaij;
373     ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
374     ierr = PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
375     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
376   }
377   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
378 
379   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
380   hdiag = hypre_ParCSRMatrixDiag(parcsr);
381   hoffd = hypre_ParCSRMatrixOffd(parcsr);
382   m     = hypre_CSRMatrixNumRows(hdiag);
383   n     = hypre_CSRMatrixNumCols(hdiag);
384   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
385   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
386   if (reuse == MAT_INITIAL_MATRIX) {
387     ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr);
388     ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr);
389     ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr);
390   } else if (reuse == MAT_REUSE_MATRIX) {
391     PetscInt  nr;
392     PetscBool done;
393     if (size > 1) {
394       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
395 
396       ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
397       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);
398       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);
399       ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr);
400     } else {
401       ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
402       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
403       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);
404       ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr);
405     }
406   } else { /* MAT_INPLACE_MATRIX */
407     dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
408     djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
409     da  = hypre_CSRMatrixData(hdiag);
410   }
411   ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
412   ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr);
413   ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr);
414   iptr = djj;
415   aptr = da;
416   for (i=0; i<m; i++) {
417     PetscInt nc = dii[i+1]-dii[i];
418     ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
419     iptr += nc;
420     aptr += nc;
421   }
422   if (size > 1) {
423     HYPRE_Int *offdj,*coffd;
424 
425     if (reuse == MAT_INITIAL_MATRIX) {
426       ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr);
427       ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr);
428       ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr);
429     } else if (reuse == MAT_REUSE_MATRIX) {
430       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
431       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
432       PetscBool  done;
433 
434       ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr);
435       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);
436       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);
437       ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr);
438     } else { /* MAT_INPLACE_MATRIX */
439       oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
440       ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
441       oa  = hypre_CSRMatrixData(hoffd);
442     }
443     ierr  = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
444     offdj = hypre_CSRMatrixJ(hoffd);
445     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
446     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
447     ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr);
448     iptr = ojj;
449     aptr = oa;
450     for (i=0; i<m; i++) {
451        PetscInt nc = oii[i+1]-oii[i];
452        ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
453        iptr += nc;
454        aptr += nc;
455     }
456     if (reuse == MAT_INITIAL_MATRIX) {
457       Mat_MPIAIJ *b;
458       Mat_SeqAIJ *d,*o;
459 
460       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);CHKERRQ(ierr);
461       /* hack MPIAIJ */
462       b          = (Mat_MPIAIJ*)((*B)->data);
463       d          = (Mat_SeqAIJ*)b->A->data;
464       o          = (Mat_SeqAIJ*)b->B->data;
465       d->free_a  = PETSC_TRUE;
466       d->free_ij = PETSC_TRUE;
467       o->free_a  = PETSC_TRUE;
468       o->free_ij = PETSC_TRUE;
469     } else if (reuse == MAT_INPLACE_MATRIX) {
470       Mat T;
471       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);CHKERRQ(ierr);
472       hypre_CSRMatrixI(hdiag)    = NULL;
473       hypre_CSRMatrixJ(hdiag)    = NULL;
474       hypre_CSRMatrixData(hdiag) = NULL;
475       hypre_CSRMatrixI(hoffd)    = NULL;
476       hypre_CSRMatrixJ(hoffd)    = NULL;
477       hypre_CSRMatrixData(hoffd) = NULL;
478       ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr);
479     }
480   } else {
481     oii  = NULL;
482     ojj  = NULL;
483     oa   = NULL;
484     if (reuse == MAT_INITIAL_MATRIX) {
485       Mat_SeqAIJ* b;
486       ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr);
487       /* hack SeqAIJ */
488       b          = (Mat_SeqAIJ*)((*B)->data);
489       b->free_a  = PETSC_TRUE;
490       b->free_ij = PETSC_TRUE;
491     } else if (reuse == MAT_INPLACE_MATRIX) {
492       Mat T;
493       ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);CHKERRQ(ierr);
494       hypre_CSRMatrixI(hdiag)    = NULL;
495       hypre_CSRMatrixJ(hdiag)    = NULL;
496       hypre_CSRMatrixData(hdiag) = NULL;
497       ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr);
498     }
499   }
500 
501   /* we have to use hypre_Tfree to free the arrays */
502   if (reuse == MAT_INPLACE_MATRIX) {
503     void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
504     const char *names[6] = {"_hypre_csr_dii",
505                             "_hypre_csr_djj",
506                             "_hypre_csr_da",
507                             "_hypre_csr_oii",
508                             "_hypre_csr_ojj",
509                             "_hypre_csr_oa"};
510     for (i=0; i<6; i++) {
511       PetscContainer c;
512 
513       ierr = PetscContainerCreate(comm,&c);CHKERRQ(ierr);
514       ierr = PetscContainerSetPointer(c,ptrs[i]);CHKERRQ(ierr);
515       ierr = PetscContainerSetUserDestroy(c,hypre_array_destroy);CHKERRQ(ierr);
516       ierr = PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);CHKERRQ(ierr);
517       ierr = PetscContainerDestroy(&c);CHKERRQ(ierr);
518     }
519   }
520   PetscFunctionReturn(0);
521 }
522 
523 static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
524 {
525   hypre_ParCSRMatrix *tA;
526   hypre_CSRMatrix    *hdiag,*hoffd;
527   Mat_SeqAIJ         *diag,*offd;
528   PetscInt           *garray,noffd,dnnz,onnz,*row_starts,*col_starts;
529   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
530   PetscBool          ismpiaij,isseqaij;
531   PetscErrorCode     ierr;
532 
533   PetscFunctionBegin;
534   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
535   ierr = PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
536   if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type);
537   if (ismpiaij) {
538     Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);
539 
540     diag   = (Mat_SeqAIJ*)a->A->data;
541     offd   = (Mat_SeqAIJ*)a->B->data;
542     garray = a->garray;
543     noffd  = a->B->cmap->N;
544     dnnz   = diag->nz;
545     onnz   = offd->nz;
546   } else {
547     diag   = (Mat_SeqAIJ*)A->data;
548     offd   = NULL;
549     garray = NULL;
550     noffd  = 0;
551     dnnz   = diag->nz;
552     onnz   = 0;
553   }
554 
555   /* create a temporary ParCSR */
556   if (HYPRE_AssumedPartitionCheck()) {
557     PetscMPIInt myid;
558 
559     ierr       = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr);
560     row_starts = A->rmap->range + myid;
561     col_starts = A->cmap->range + myid;
562   } else {
563     row_starts = A->rmap->range;
564     col_starts = A->cmap->range;
565   }
566   tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_Int*)row_starts,(HYPRE_Int*)col_starts,noffd,dnnz,onnz);
567   hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
568   hypre_ParCSRMatrixSetColStartsOwner(tA,0);
569 
570   /* set diagonal part */
571   hdiag = hypre_ParCSRMatrixDiag(tA);
572   hypre_CSRMatrixI(hdiag)           = (HYPRE_Int*)diag->i;
573   hypre_CSRMatrixJ(hdiag)           = (HYPRE_Int*)diag->j;
574   hypre_CSRMatrixData(hdiag)        = diag->a;
575   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
576   hypre_CSRMatrixSetRownnz(hdiag);
577   hypre_CSRMatrixSetDataOwner(hdiag,0);
578 
579   /* set offdiagonal part */
580   hoffd = hypre_ParCSRMatrixOffd(tA);
581   if (offd) {
582     hypre_CSRMatrixI(hoffd)           = (HYPRE_Int*)offd->i;
583     hypre_CSRMatrixJ(hoffd)           = (HYPRE_Int*)offd->j;
584     hypre_CSRMatrixData(hoffd)        = offd->a;
585     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
586     hypre_CSRMatrixSetRownnz(hoffd);
587     hypre_CSRMatrixSetDataOwner(hoffd,0);
588     hypre_ParCSRMatrixSetNumNonzeros(tA);
589     hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_Int*)garray;
590   }
591   *hA = tA;
592   PetscFunctionReturn(0);
593 }
594 
595 static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
596 {
597   hypre_CSRMatrix    *hdiag,*hoffd;
598 
599   PetscFunctionBegin;
600   hdiag = hypre_ParCSRMatrixDiag(*hA);
601   hoffd = hypre_ParCSRMatrixOffd(*hA);
602   /* set pointers to NULL before destroying tA */
603   hypre_CSRMatrixI(hdiag)           = NULL;
604   hypre_CSRMatrixJ(hdiag)           = NULL;
605   hypre_CSRMatrixData(hdiag)        = NULL;
606   hypre_CSRMatrixI(hoffd)           = NULL;
607   hypre_CSRMatrixJ(hoffd)           = NULL;
608   hypre_CSRMatrixData(hoffd)        = NULL;
609   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
610   hypre_ParCSRMatrixDestroy(*hA);
611   *hA = NULL;
612   PetscFunctionReturn(0);
613 }
614 
615 /* calls RAP from BoomerAMG:
616    the resulting ParCSR will not own the column and row starts
617    It looks like we don't need to have the diagonal entries
618    ordered first in the rows of the diagonal part
619    for boomerAMGBuildCoarseOperator to work */
620 static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
621 {
622   PetscErrorCode ierr;
623   HYPRE_Int      P_owns_col_starts,R_owns_row_starts;
624 
625   PetscFunctionBegin;
626   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
627   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
628   PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP));
629   PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP));
630   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
631   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
632   hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
633   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
634   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
635   PetscFunctionReturn(0);
636 }
637 
638 static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
639 {
640   Mat                B;
641   hypre_ParCSRMatrix *hA,*hP,*hPtAP;
642   PetscErrorCode     ierr;
643 
644   PetscFunctionBegin;
645   ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr);
646   ierr = MatAIJGetParCSR_Private(P,&hP);CHKERRQ(ierr);
647   ierr = MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);CHKERRQ(ierr);
648   ierr = MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);CHKERRQ(ierr);
649   ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr);
650   ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr);
651   ierr = MatAIJRestoreParCSR_Private(P,&hP);CHKERRQ(ierr);
652   PetscFunctionReturn(0);
653 }
654 
655 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat *C)
656 {
657   PetscErrorCode ierr;
658 
659   PetscFunctionBegin;
660   ierr                   = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr);
661   ierr                   = MatSetType(*C,MATAIJ);CHKERRQ(ierr);
662   (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
663   PetscFunctionReturn(0);
664 }
665 
666 static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
667 {
668   Mat                B;
669   Mat_HYPRE          *hP;
670   hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr;
671   HYPRE_Int          type;
672   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
673   PetscBool          ishypre;
674   PetscErrorCode     ierr;
675 
676   PetscFunctionBegin;
677   ierr = PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);CHKERRQ(ierr);
678   if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
679   hP = (Mat_HYPRE*)P->data;
680   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
681   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
682   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
683 
684   ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr);
685   ierr = MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);CHKERRQ(ierr);
686   ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr);
687 
688   /* create temporary matrix and merge to C */
689   ierr = MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);CHKERRQ(ierr);
690   ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr);
691   PetscFunctionReturn(0);
692 }
693 
694 static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
695 {
696   PetscErrorCode ierr;
697 
698   PetscFunctionBegin;
699   if (scall == MAT_INITIAL_MATRIX) {
700     const char *deft = MATAIJ;
701     char       type[256];
702     PetscBool  flg;
703 
704     ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr);
705     ierr = PetscOptionsFList("-matptap_hypre_outtype","Matrix type",NULL,MatList,deft,type,256,&flg);CHKERRQ(ierr);
706     ierr = PetscOptionsEnd();CHKERRQ(ierr);
707     ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
708     ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr);
709     if (flg) {
710       ierr = MatSetType(*C,type);CHKERRQ(ierr);
711     } else {
712       ierr = MatSetType(*C,deft);CHKERRQ(ierr);
713     }
714     (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
715     ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
716   }
717   ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
718   ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr);
719   ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
720   PetscFunctionReturn(0);
721 }
722 
723 static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
724 {
725   Mat                B;
726   hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr;
727   Mat_HYPRE          *hA,*hP;
728   PetscBool          ishypre;
729   HYPRE_Int          type;
730   PetscErrorCode     ierr;
731 
732   PetscFunctionBegin;
733   ierr = PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);CHKERRQ(ierr);
734   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
735   ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr);
736   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
737   hA = (Mat_HYPRE*)A->data;
738   hP = (Mat_HYPRE*)P->data;
739   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
740   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
741   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
742   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
743   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
744   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
745   ierr = MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);CHKERRQ(ierr);
746   ierr = MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);CHKERRQ(ierr);
747   ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr);
748   PetscFunctionReturn(0);
749 }
750 
751 static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
752 {
753   PetscErrorCode ierr;
754 
755   PetscFunctionBegin;
756   if (scall == MAT_INITIAL_MATRIX) {
757     ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
758     ierr                   = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr);
759     ierr                   = MatSetType(*C,MATHYPRE);CHKERRQ(ierr);
760     (*C)->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
761     ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
762   }
763   ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
764   ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr);
765   ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
766   PetscFunctionReturn(0);
767 }
768 
769 /* calls hypre_ParMatmul
770    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
771    hypre_ParMatrixCreate does not duplicate the communicator
772    It looks like we don't need to have the diagonal entries
773    ordered first in the rows of the diagonal part
774    for boomerAMGBuildCoarseOperator to work */
775 static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
776 {
777   PetscFunctionBegin;
778   PetscStackPush("hypre_ParMatmul");
779   *hAB = hypre_ParMatmul(hA,hB);
780   PetscStackPop;
781   PetscFunctionReturn(0);
782 }
783 
784 static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
785 {
786   Mat                D;
787   hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
788   PetscErrorCode     ierr;
789 
790   PetscFunctionBegin;
791   ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr);
792   ierr = MatAIJGetParCSR_Private(B,&hB);CHKERRQ(ierr);
793   ierr = MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);CHKERRQ(ierr);
794   ierr = MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);CHKERRQ(ierr);
795   ierr = MatHeaderMerge(C,&D);CHKERRQ(ierr);
796   ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr);
797   ierr = MatAIJRestoreParCSR_Private(B,&hB);CHKERRQ(ierr);
798   PetscFunctionReturn(0);
799 }
800 
801 PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat *C)
802 {
803   PetscErrorCode ierr;
804 
805   PetscFunctionBegin;
806   ierr                      = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr);
807   ierr                      = MatSetType(*C,MATAIJ);CHKERRQ(ierr);
808   (*C)->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
809   PetscFunctionReturn(0);
810 }
811 
812 static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
813 {
814   Mat                D;
815   hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
816   Mat_HYPRE          *hA,*hB;
817   PetscBool          ishypre;
818   HYPRE_Int          type;
819   PetscErrorCode     ierr;
820 
821   PetscFunctionBegin;
822   ierr = PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);CHKERRQ(ierr);
823   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
824   ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr);
825   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
826   hA = (Mat_HYPRE*)A->data;
827   hB = (Mat_HYPRE*)B->data;
828   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
829   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
830   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
831   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
832   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
833   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
834   ierr = MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);CHKERRQ(ierr);
835   ierr = MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);CHKERRQ(ierr);
836   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
837   ierr = MatHeaderReplace(C,&D);CHKERRQ(ierr);
838   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
839   PetscFunctionReturn(0);
840 }
841 
842 static PetscErrorCode MatMatMult_HYPRE_HYPRE(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
843 {
844   PetscErrorCode ierr;
845 
846   PetscFunctionBegin;
847   if (scall == MAT_INITIAL_MATRIX) {
848     ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr);
849     ierr                      = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr);
850     ierr                      = MatSetType(*C,MATHYPRE);CHKERRQ(ierr);
851     (*C)->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
852     ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr);
853   }
854   ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr);
855   ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr);
856   ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr);
857   PetscFunctionReturn(0);
858 }
859 
860 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
861 {
862   Mat                E;
863   hypre_ParCSRMatrix *hA,*hB,*hC,*hABC;
864   PetscErrorCode     ierr;
865 
866   PetscFunctionBegin;
867   ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr);
868   ierr = MatAIJGetParCSR_Private(B,&hB);CHKERRQ(ierr);
869   ierr = MatAIJGetParCSR_Private(C,&hC);CHKERRQ(ierr);
870   ierr = MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);CHKERRQ(ierr);
871   ierr = MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);CHKERRQ(ierr);
872   ierr = MatHeaderMerge(D,&E);CHKERRQ(ierr);
873   ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr);
874   ierr = MatAIJRestoreParCSR_Private(B,&hB);CHKERRQ(ierr);
875   ierr = MatAIJRestoreParCSR_Private(C,&hC);CHKERRQ(ierr);
876   PetscFunctionReturn(0);
877 }
878 
879 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
880 {
881   PetscErrorCode ierr;
882 
883   PetscFunctionBegin;
884   ierr = MatCreate(PetscObjectComm((PetscObject)A),D);CHKERRQ(ierr);
885   ierr = MatSetType(*D,MATAIJ);CHKERRQ(ierr);
886   PetscFunctionReturn(0);
887 }
888 
889 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
890 {
891   PetscErrorCode ierr;
892 
893   PetscFunctionBegin;
894   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);CHKERRQ(ierr);
895   PetscFunctionReturn(0);
896 }
897 
898 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
899 {
900   PetscErrorCode ierr;
901 
902   PetscFunctionBegin;
903   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);CHKERRQ(ierr);
904   PetscFunctionReturn(0);
905 }
906 
907 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
908 {
909   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
910   hypre_ParCSRMatrix *parcsr;
911   hypre_ParVector    *hx,*hy;
912   PetscScalar        *ax,*ay,*sax,*say;
913   PetscErrorCode     ierr;
914 
915   PetscFunctionBegin;
916   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
917   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
918   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
919   ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
920   ierr = VecGetArray(y,&ay);CHKERRQ(ierr);
921   if (trans) {
922     VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
923     VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
924     hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
925     VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
926     VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
927   } else {
928     VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
929     VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
930     hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
931     VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
932     VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
933   }
934   ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
935   ierr = VecRestoreArray(y,&ay);CHKERRQ(ierr);
936   PetscFunctionReturn(0);
937 }
938 
939 static PetscErrorCode MatDestroy_HYPRE(Mat A)
940 {
941   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
942   PetscErrorCode ierr;
943 
944   PetscFunctionBegin;
945   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
946   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
947   if (hA->ij) {
948     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
949     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
950   }
951   if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); }
952   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr);
953   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr);
954   ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);CHKERRQ(ierr);
955   ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);CHKERRQ(ierr);
956   ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);CHKERRQ(ierr);
957   ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);CHKERRQ(ierr);
958   ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_hypre_C",NULL);CHKERRQ(ierr);
959   ierr = PetscFree(A->data);CHKERRQ(ierr);
960   PetscFunctionReturn(0);
961 }
962 
963 static PetscErrorCode MatSetUp_HYPRE(Mat A)
964 {
965   PetscErrorCode ierr;
966 
967   PetscFunctionBegin;
968   ierr = MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);CHKERRQ(ierr);
969   PetscFunctionReturn(0);
970 }
971 
972 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
973 {
974   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
975   Vec                x,b;
976   PetscErrorCode     ierr;
977 
978   PetscFunctionBegin;
979   if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");
980   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
981   if (hA->x) PetscFunctionReturn(0);
982   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
983   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
984   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr);
985   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr);
986   ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr);
987   ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr);
988   ierr = VecDestroy(&x);CHKERRQ(ierr);
989   ierr = VecDestroy(&b);CHKERRQ(ierr);
990   PetscFunctionReturn(0);
991 }
992 
993 #define MATHYPRE_SCRATCH 2048
994 
995 PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
996 {
997   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
998   PetscScalar        *vals = (PetscScalar *)v;
999   PetscScalar        sscr[MATHYPRE_SCRATCH];
1000   HYPRE_Int          cscr[2][MATHYPRE_SCRATCH];
1001   HYPRE_Int          i,nzc;
1002   PetscErrorCode     ierr;
1003 
1004   PetscFunctionBegin;
1005   for (i=0,nzc=0;i<nc;i++) {
1006     if (cols[i] >= 0) {
1007       cscr[0][nzc  ] = cols[i];
1008       cscr[1][nzc++] = i;
1009     }
1010   }
1011   if (!nzc) PetscFunctionReturn(0);
1012 
1013   if (ins == ADD_VALUES) {
1014     for (i=0;i<nr;i++) {
1015       if (rows[i] >= 0 && nzc) {
1016         PetscInt j;
1017         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1018         PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&nzc,(HYPRE_Int*)(rows+i),cscr[0],sscr));
1019       }
1020       vals += nc;
1021     }
1022   } else { /* INSERT_VALUES */
1023 #if defined(PETSC_USE_DEBUG)
1024     /* Insert values cannot be used to insert offproc entries */
1025     PetscInt rst,ren;
1026     ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
1027     for (i=0;i<nr;i++)
1028       if (rows[i] >= 0 && (rows[i] < rst || rows[i] >= ren)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use INSERT_VALUES for off-proc entries with MatHYPRE. Use ADD_VALUES instead");
1029 #endif
1030     for (i=0;i<nr;i++) {
1031       if (rows[i] >= 0 && nzc) {
1032         PetscInt j;
1033         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1034         PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&nzc,(HYPRE_Int*)(rows+i),cscr[0],sscr));
1035       }
1036       vals += nc;
1037     }
1038   }
1039   PetscFunctionReturn(0);
1040 }
1041 
1042 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1043 {
1044   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1045   HYPRE_Int          *hdnnz,*honnz;
1046   PetscInt           i,rs,re,cs,ce,bs;
1047   PetscMPIInt        size;
1048   PetscErrorCode     ierr;
1049 
1050   PetscFunctionBegin;
1051   ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
1052   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
1053   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
1054   rs   = A->rmap->rstart;
1055   re   = A->rmap->rend;
1056   cs   = A->cmap->rstart;
1057   ce   = A->cmap->rend;
1058   if (!hA->ij) {
1059     PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1060     PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1061   } else {
1062     HYPRE_Int hrs,hre,hcs,hce;
1063     PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1064     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);
1065     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);
1066   }
1067   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1068 
1069   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1070   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;
1071 
1072   if (!dnnz) {
1073     ierr = PetscMalloc1(A->rmap->n,&hdnnz);CHKERRQ(ierr);
1074     for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1075   } else {
1076     hdnnz = (HYPRE_Int*)dnnz;
1077   }
1078   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
1079   if (size > 1) {
1080     if (!onnz) {
1081       ierr = PetscMalloc1(A->rmap->n,&honnz);CHKERRQ(ierr);
1082       for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1083     } else {
1084       honnz = (HYPRE_Int*)onnz;
1085     }
1086     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1087   } else {
1088     honnz = NULL;
1089     PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1090   }
1091   if (!dnnz) {
1092     ierr = PetscFree(hdnnz);CHKERRQ(ierr);
1093   }
1094   if (!onnz && honnz) {
1095     ierr = PetscFree(honnz);CHKERRQ(ierr);
1096   }
1097   A->preallocated = PETSC_TRUE;
1098 
1099   /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0 */
1100   {
1101     hypre_AuxParCSRMatrix *aux_matrix;
1102     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1103     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
1104   }
1105   PetscFunctionReturn(0);
1106 }
1107 
1108 /*@C
1109    MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1110 
1111    Collective on Mat
1112 
1113    Input Parameters:
1114 +  A - the matrix
1115 .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1116           (same value is used for all local rows)
1117 .  dnnz - array containing the number of nonzeros in the various rows of the
1118           DIAGONAL portion of the local submatrix (possibly different for each row)
1119           or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1120           The size of this array is equal to the number of local rows, i.e 'm'.
1121           For matrices that will be factored, you must leave room for (and set)
1122           the diagonal entry even if it is zero.
1123 .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1124           submatrix (same value is used for all local rows).
1125 -  onnz - array containing the number of nonzeros in the various rows of the
1126           OFF-DIAGONAL portion of the local submatrix (possibly different for
1127           each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1128           structure. The size of this array is equal to the number
1129           of local rows, i.e 'm'.
1130 
1131    Notes:
1132     If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored.
1133 
1134    Level: intermediate
1135 
1136 .keywords: matrix, aij, compressed row, sparse, parallel
1137 
1138 .seealso: MatCreate(), MatMPIAIJSetPreallocation, MATHYPRE
1139 @*/
1140 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1141 {
1142   PetscErrorCode ierr;
1143 
1144   PetscFunctionBegin;
1145   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1146   PetscValidType(A,1);
1147   ierr = PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));CHKERRQ(ierr);
1148   PetscFunctionReturn(0);
1149 }
1150 
1151 /*
1152    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix
1153 
1154    Collective
1155 
1156    Input Parameters:
1157 +  vparcsr  - the pointer to the hypre_ParCSRMatrix
1158 .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1159 -  copymode - PETSc copying options
1160 
1161    Output Parameter:
1162 .  A  - the matrix
1163 
1164    Level: intermediate
1165 
1166 .seealso: MatHYPRE, PetscCopyMode
1167 */
1168 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *vparcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1169 {
1170   Mat                   T;
1171   Mat_HYPRE             *hA;
1172   hypre_ParCSRMatrix    *parcsr;
1173   MPI_Comm              comm;
1174   PetscInt              rstart,rend,cstart,cend,M,N;
1175   PetscBool             isseqaij,ismpiaij,isaij,ishyp,isis;
1176   PetscErrorCode        ierr;
1177 
1178   PetscFunctionBegin;
1179   parcsr = (hypre_ParCSRMatrix *)vparcsr;
1180   comm   = hypre_ParCSRMatrixComm(parcsr);
1181   ierr   = PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
1182   ierr   = PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
1183   ierr   = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr);
1184   ierr   = PetscStrcmp(mtype,MATHYPRE,&ishyp);CHKERRQ(ierr);
1185   ierr   = PetscStrcmp(mtype,MATIS,&isis);CHKERRQ(ierr);
1186   isaij  = (PetscBool)(isseqaij || ismpiaij || isaij);
1187   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);
1188   if (ishyp && copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");
1189 
1190   /* access ParCSRMatrix */
1191   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1192   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1193   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1194   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1195   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1196   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1197 
1198   /* fix for empty local rows/columns */
1199   if (rend < rstart) rend = rstart;
1200   if (cend < cstart) cend = cstart;
1201 
1202   /* create PETSc matrix with MatHYPRE */
1203   ierr = MatCreate(comm,&T);CHKERRQ(ierr);
1204   ierr = MatSetSizes(T,rend-rstart+1,cend-cstart+1,M,N);CHKERRQ(ierr);
1205   ierr = MatSetType(T,MATHYPRE);CHKERRQ(ierr);
1206   hA   = (Mat_HYPRE*)(T->data);
1207 
1208   /* create HYPRE_IJMatrix */
1209   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
1210 
1211   /* set ParCSR object */
1212   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1213   hypre_IJMatrixObject(hA->ij) = parcsr;
1214   T->preallocated = PETSC_TRUE;
1215 
1216   /* set assembled flag */
1217   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1218   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1219   if (ishyp) {
1220     PetscMPIInt myid = 0;
1221 
1222     /* make sure we always have row_starts and col_starts available */
1223     if (HYPRE_AssumedPartitionCheck()) {
1224       ierr = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr);
1225     }
1226     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1227       PetscLayout map;
1228 
1229       ierr = MatGetLayouts(T,NULL,&map);CHKERRQ(ierr);
1230       ierr = PetscLayoutSetUp(map);CHKERRQ(ierr);
1231       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1232     }
1233     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1234       PetscLayout map;
1235 
1236       ierr = MatGetLayouts(T,&map,NULL);CHKERRQ(ierr);
1237       ierr = PetscLayoutSetUp(map);CHKERRQ(ierr);
1238       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1239     }
1240     /* prevent from freeing the pointer */
1241     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1242     *A   = T;
1243     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1244     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1245   } else if (isaij) {
1246     if (copymode != PETSC_OWN_POINTER) {
1247       /* prevent from freeing the pointer */
1248       hA->inner_free = PETSC_FALSE;
1249       ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr);
1250       ierr = MatDestroy(&T);CHKERRQ(ierr);
1251     } else { /* AIJ return type with PETSC_OWN_POINTER */
1252       ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr);
1253       *A   = T;
1254     }
1255   } else if (isis) {
1256     ierr = MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr);
1257     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1258     ierr = MatDestroy(&T);CHKERRQ(ierr);
1259   }
1260   PetscFunctionReturn(0);
1261 }
1262 
1263 static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1264 {
1265   Mat_HYPRE*            hA = (Mat_HYPRE*)A->data;
1266   HYPRE_Int             type;
1267   PetscErrorCode        ierr;
1268 
1269   PetscFunctionBegin;
1270   if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1271   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1272   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1273   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1274   PetscFunctionReturn(0);
1275 }
1276 
1277 /*
1278    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1279 
1280    Not collective
1281 
1282    Input Parameters:
1283 +  A  - the MATHYPRE object
1284 
1285    Output Parameter:
1286 .  parcsr  - the pointer to the hypre_ParCSRMatrix
1287 
1288    Level: intermediate
1289 
1290 .seealso: MatHYPRE, PetscCopyMode
1291 */
1292 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1293 {
1294   PetscErrorCode ierr;
1295 
1296   PetscFunctionBegin;
1297   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1298   PetscValidType(A,1);
1299   ierr = PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));CHKERRQ(ierr);
1300   PetscFunctionReturn(0);
1301 }
1302 
1303 static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1304 {
1305   hypre_ParCSRMatrix *parcsr;
1306   hypre_CSRMatrix    *ha;
1307   PetscInt           rst;
1308   PetscErrorCode     ierr;
1309 
1310   PetscFunctionBegin;
1311   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1312   ierr = MatGetOwnershipRange(A,&rst,NULL);CHKERRQ(ierr);
1313   ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr);
1314   if (missing) *missing = PETSC_FALSE;
1315   if (dd) *dd = -1;
1316   ha = hypre_ParCSRMatrixDiag(parcsr);
1317   if (ha) {
1318     PetscInt  size,i;
1319     HYPRE_Int *ii,*jj;
1320 
1321     size = hypre_CSRMatrixNumRows(ha);
1322     ii   = hypre_CSRMatrixI(ha);
1323     jj   = hypre_CSRMatrixJ(ha);
1324     for (i = 0; i < size; i++) {
1325       PetscInt  j;
1326       PetscBool found = PETSC_FALSE;
1327 
1328       for (j = ii[i]; j < ii[i+1] && !found; j++)
1329         found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE;
1330 
1331       if (!found) {
1332         PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1333         if (missing) *missing = PETSC_TRUE;
1334         if (dd) *dd = i+rst;
1335         PetscFunctionReturn(0);
1336       }
1337     }
1338     if (!size) {
1339       PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1340       if (missing) *missing = PETSC_TRUE;
1341       if (dd) *dd = rst;
1342     }
1343   } else {
1344     PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1345     if (missing) *missing = PETSC_TRUE;
1346     if (dd) *dd = rst;
1347   }
1348   PetscFunctionReturn(0);
1349 }
1350 
1351 static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1352 {
1353   hypre_ParCSRMatrix *parcsr;
1354   hypre_CSRMatrix    *ha;
1355   PetscErrorCode     ierr;
1356 
1357   PetscFunctionBegin;
1358   ierr  = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr);
1359   /* diagonal part */
1360   ha = hypre_ParCSRMatrixDiag(parcsr);
1361   if (ha) {
1362     PetscInt    size,i;
1363     HYPRE_Int   *ii;
1364     PetscScalar *a;
1365 
1366     size = hypre_CSRMatrixNumRows(ha);
1367     a    = hypre_CSRMatrixData(ha);
1368     ii   = hypre_CSRMatrixI(ha);
1369     for (i = 0; i < ii[size]; i++) a[i] *= s;
1370   }
1371   /* offdiagonal part */
1372   ha = hypre_ParCSRMatrixOffd(parcsr);
1373   if (ha) {
1374     PetscInt    size,i;
1375     HYPRE_Int   *ii;
1376     PetscScalar *a;
1377 
1378     size = hypre_CSRMatrixNumRows(ha);
1379     a    = hypre_CSRMatrixData(ha);
1380     ii   = hypre_CSRMatrixI(ha);
1381     for (i = 0; i < ii[size]; i++) a[i] *= s;
1382   }
1383   PetscFunctionReturn(0);
1384 }
1385 
1386 static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1387 {
1388   hypre_ParCSRMatrix *parcsr;
1389   HYPRE_Int          *lrows;
1390   PetscInt           rst,ren,i;
1391   PetscErrorCode     ierr;
1392 
1393   PetscFunctionBegin;
1394   if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1395   ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr);
1396   ierr = PetscMalloc1(numRows,&lrows);CHKERRQ(ierr);
1397   ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
1398   for (i=0;i<numRows;i++) {
1399     if (rows[i] < rst || rows[i] >= ren)
1400       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1401     lrows[i] = rows[i] - rst;
1402   }
1403   PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1404   ierr = PetscFree(lrows);CHKERRQ(ierr);
1405   PetscFunctionReturn(0);
1406 }
1407 
1408 /*MC
1409    MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
1410           based on the hypre IJ interface.
1411 
1412    Level: intermediate
1413 
1414 .seealso: MatCreate()
1415 M*/
1416 
1417 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
1418 {
1419   Mat_HYPRE      *hB;
1420   PetscErrorCode ierr;
1421 
1422   PetscFunctionBegin;
1423   ierr           = PetscNewLog(B,&hB);CHKERRQ(ierr);
1424   hB->inner_free = PETSC_TRUE;
1425 
1426   B->data       = (void*)hB;
1427   B->rmap->bs   = 1;
1428   B->assembled  = PETSC_FALSE;
1429 
1430   ierr = PetscMemzero(B->ops,sizeof(struct _MatOps));CHKERRQ(ierr);
1431   B->ops->mult            = MatMult_HYPRE;
1432   B->ops->multtranspose   = MatMultTranspose_HYPRE;
1433   B->ops->setup           = MatSetUp_HYPRE;
1434   B->ops->destroy         = MatDestroy_HYPRE;
1435   B->ops->assemblyend     = MatAssemblyEnd_HYPRE;
1436   B->ops->ptap            = MatPtAP_HYPRE_HYPRE;
1437   B->ops->matmult         = MatMatMult_HYPRE_HYPRE;
1438   B->ops->setvalues       = MatSetValues_HYPRE;
1439   B->ops->missingdiagonal = MatMissingDiagonal_HYPRE;
1440   B->ops->scale           = MatScale_HYPRE;
1441   B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
1442 
1443   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr);
1444   ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr);
1445   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr);
1446   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr);
1447   ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr);
1448   ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr);
1449   ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);CHKERRQ(ierr);
1450   ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);CHKERRQ(ierr);
1451   ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_hypre_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr);
1452   PetscFunctionReturn(0);
1453 }
1454 
1455 static PetscErrorCode hypre_array_destroy(void *ptr)
1456 {
1457    PetscFunctionBegin;
1458    hypre_TFree(ptr,HYPRE_MEMORY_HOST);
1459    PetscFunctionReturn(0);
1460 }
1461