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