xref: /petsc/src/mat/impls/hypre/mhypre.c (revision 41073d710092a75ce18d397101144c768a8ad705) 
1 
2 /*
3     Creates hypre ijmatrix from PETSc matrix
4 */
5 
6 /*MC
7    MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
8           based on the hypre IJ interface.
9 
10    Level: intermediate
11 
12 .seealso: MatCreate()
13 M*/
14 
15 #include <petscsys.h>
16 #include <petsc/private/matimpl.h>
17 #include <../src/mat/impls/hypre/mhypre.h>
18 #include <../src/mat/impls/aij/mpi/mpiaij.h>
19 #include <../src/vec/vec/impls/hypre/vhyp.h>
20 #include <HYPRE.h>
21 #include <HYPRE_utilities.h>
22 #include <_hypre_parcsr_ls.h>
23 
24 static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
25 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
26 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
27 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
28 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,Vec,Vec,PetscBool);
29 static PetscErrorCode hypre_array_destroy(void*);
30 
31 #undef __FUNCT__
32 #define __FUNCT__ "MatHYPRE_IJMatrixPreallocate"
33 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
34 {
35   PetscErrorCode ierr;
36   PetscInt       i,n_d,n_o;
37   const PetscInt *ia_d,*ia_o;
38   PetscBool      done_d=PETSC_FALSE,done_o=PETSC_FALSE;
39   PetscInt       *nnz_d=NULL,*nnz_o=NULL;
40 
41   PetscFunctionBegin;
42   if (A_d) { /* determine number of nonzero entries in local diagonal part */
43     ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);CHKERRQ(ierr);
44     if (done_d) {
45       ierr = PetscMalloc1(n_d,&nnz_d);CHKERRQ(ierr);
46       for (i=0; i<n_d; i++) {
47         nnz_d[i] = ia_d[i+1] - ia_d[i];
48       }
49     }
50     ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);CHKERRQ(ierr);
51   }
52   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
53     ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
54     if (done_o) {
55       ierr = PetscMalloc1(n_o,&nnz_o);CHKERRQ(ierr);
56       for (i=0; i<n_o; i++) {
57         nnz_o[i] = ia_o[i+1] - ia_o[i];
58       }
59     }
60     ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
61   }
62   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
63     if (!done_o) { /* only diagonal part */
64       ierr = PetscMalloc1(n_d,&nnz_o);CHKERRQ(ierr);
65       for (i=0; i<n_d; i++) {
66         nnz_o[i] = 0;
67       }
68     }
69     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
70     ierr = PetscFree(nnz_d);CHKERRQ(ierr);
71     ierr = PetscFree(nnz_o);CHKERRQ(ierr);
72   }
73   PetscFunctionReturn(0);
74 }
75 
76 #undef __FUNCT__
77 #define __FUNCT__ "MatHYPRE_CreateFromMat"
78 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
79 {
80   PetscErrorCode ierr;
81   PetscInt       rstart,rend,cstart,cend;
82 
83   PetscFunctionBegin;
84   ierr   = MatSetUp(A);CHKERRQ(ierr);
85   rstart = A->rmap->rstart;
86   rend   = A->rmap->rend;
87   cstart = A->cmap->rstart;
88   cend   = A->cmap->rend;
89   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
90   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
91   {
92     PetscBool      same;
93     Mat            A_d,A_o;
94     const PetscInt *colmap;
95     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr);
96     if (same) {
97       ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
98       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
99       PetscFunctionReturn(0);
100     }
101     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr);
102     if (same) {
103       ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
104       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
105       PetscFunctionReturn(0);
106     }
107     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr);
108     if (same) {
109       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
110       PetscFunctionReturn(0);
111     }
112     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr);
113     if (same) {
114       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
115       PetscFunctionReturn(0);
116     }
117   }
118   PetscFunctionReturn(0);
119 }
120 
121 #undef __FUNCT__
122 #define __FUNCT__ "MatHYPRE_IJMatrixCopy"
123 static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
124 {
125   PetscErrorCode    ierr;
126   PetscInt          i,rstart,rend,ncols,nr,nc;
127   const PetscScalar *values;
128   const PetscInt    *cols;
129   PetscBool         flg;
130 
131   PetscFunctionBegin;
132   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
133   ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr);
134   if (flg && nr == nc) {
135     ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr);
136     PetscFunctionReturn(0);
137   }
138   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
139   if (flg) {
140     ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr);
141     PetscFunctionReturn(0);
142   }
143 
144   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
145   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
146   for (i=rstart; i<rend; i++) {
147     ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
148     PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values));
149     ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
150   }
151   PetscFunctionReturn(0);
152 }
153 
154 #undef __FUNCT__
155 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_SeqAIJ"
156 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
157 {
158   PetscErrorCode        ierr;
159   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ*)A->data;
160   HYPRE_Int             type;
161   hypre_ParCSRMatrix    *par_matrix;
162   hypre_AuxParCSRMatrix *aux_matrix;
163   hypre_CSRMatrix       *hdiag;
164 
165   PetscFunctionBegin;
166   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
167   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
168   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
169   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
170   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
171   /*
172        this is the Hack part where we monkey directly with the hypre datastructures
173   */
174   ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
175   ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
176   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
177 
178   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
179   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
180   PetscFunctionReturn(0);
181 }
182 
183 #undef __FUNCT__
184 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_MPIAIJ"
185 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
186 {
187   PetscErrorCode        ierr;
188   Mat_MPIAIJ            *pA = (Mat_MPIAIJ*)A->data;
189   Mat_SeqAIJ            *pdiag,*poffd;
190   PetscInt              i,*garray = pA->garray,*jj,cstart,*pjj;
191   HYPRE_Int             type;
192   hypre_ParCSRMatrix    *par_matrix;
193   hypre_AuxParCSRMatrix *aux_matrix;
194   hypre_CSRMatrix       *hdiag,*hoffd;
195 
196   PetscFunctionBegin;
197   pdiag = (Mat_SeqAIJ*) pA->A->data;
198   poffd = (Mat_SeqAIJ*) pA->B->data;
199   /* cstart is only valid for square MPIAIJ layed out in the usual way */
200   ierr = MatGetOwnershipRange(A,&cstart,NULL);CHKERRQ(ierr);
201 
202   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
203   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
204   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
205   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
206   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
207   hoffd = hypre_ParCSRMatrixOffd(par_matrix);
208 
209   /*
210        this is the Hack part where we monkey directly with the hypre datastructures
211   */
212   ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
213   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
214   jj  = (PetscInt*)hdiag->j;
215   pjj = (PetscInt*)pdiag->j;
216   for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
217   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
218   ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
219   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
220      If we hacked a hypre a bit more we might be able to avoid this step */
221   jj  = (PetscInt*) hoffd->j;
222   pjj = (PetscInt*) poffd->j;
223   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
224   ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));
225 
226   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
227   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
228   PetscFunctionReturn(0);
229 }
230 
231 #undef __FUNCT__
232 #define __FUNCT__ "MatConvert_HYPRE_IS"
233 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
234 {
235   Mat_HYPRE*             mhA = (Mat_HYPRE*)(A->data);
236   Mat                    lA;
237   ISLocalToGlobalMapping rl2g,cl2g;
238   IS                     is;
239   hypre_ParCSRMatrix     *hA;
240   hypre_CSRMatrix        *hdiag,*hoffd;
241   MPI_Comm               comm;
242   PetscScalar            *hdd,*hod,*aa,*data;
243   PetscInt               *col_map_offd,*hdi,*hdj,*hoi,*hoj;
244   PetscInt               *ii,*jj,*iptr,*jptr;
245   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
246   HYPRE_Int              type;
247   PetscErrorCode         ierr;
248 
249   PetscFunctionBegin;
250   comm = PetscObjectComm((PetscObject)A);
251   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
252   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
253   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
254   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
255   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
256   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
257   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
258   hdiag = hypre_ParCSRMatrixDiag(hA);
259   hoffd = hypre_ParCSRMatrixOffd(hA);
260   dr    = hypre_CSRMatrixNumRows(hdiag);
261   dc    = hypre_CSRMatrixNumCols(hdiag);
262   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
263   hdi   = hypre_CSRMatrixI(hdiag);
264   hdj   = hypre_CSRMatrixJ(hdiag);
265   hdd   = hypre_CSRMatrixData(hdiag);
266   oc    = hypre_CSRMatrixNumCols(hoffd);
267   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
268   hoi   = hypre_CSRMatrixI(hoffd);
269   hoj   = hypre_CSRMatrixJ(hoffd);
270   hod   = hypre_CSRMatrixData(hoffd);
271   if (reuse != MAT_REUSE_MATRIX) {
272     PetscInt *aux;
273 
274     /* generate l2g maps for rows and cols */
275     ierr = ISCreateStride(comm,dr,str,1,&is);CHKERRQ(ierr);
276     ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr);
277     ierr = ISDestroy(&is);CHKERRQ(ierr);
278     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
279     ierr = PetscMalloc1(dc+oc,&aux);CHKERRQ(ierr);
280     for (i=0; i<dc; i++) aux[i] = i+stc;
281     for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
282     ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);CHKERRQ(ierr);
283     ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr);
284     ierr = ISDestroy(&is);CHKERRQ(ierr);
285     /* create MATIS object */
286     ierr = MatCreate(comm,B);CHKERRQ(ierr);
287     ierr = MatSetSizes(*B,dr,dc,M,N);CHKERRQ(ierr);
288     ierr = MatSetType(*B,MATIS);CHKERRQ(ierr);
289     ierr = MatSetLocalToGlobalMapping(*B,rl2g,cl2g);CHKERRQ(ierr);
290     ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr);
291     ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr);
292 
293     /* allocate CSR for local matrix */
294     ierr = PetscMalloc1(dr+1,&iptr);CHKERRQ(ierr);
295     ierr = PetscMalloc1(nnz,&jptr);CHKERRQ(ierr);
296     ierr = PetscMalloc1(nnz,&data);CHKERRQ(ierr);
297   } else {
298     PetscInt  nr;
299     PetscBool done;
300     ierr = MatISGetLocalMat(*B,&lA);CHKERRQ(ierr);
301     ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);CHKERRQ(ierr);
302     if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
303     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);
304     ierr = MatSeqAIJGetArray(lA,&data);CHKERRQ(ierr);
305   }
306   /* merge local matrices */
307   ii   = iptr;
308   jj   = jptr;
309   aa   = data;
310   *ii  = *(hdi++) + *(hoi++);
311   for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
312     PetscScalar *aold = aa;
313     PetscInt    *jold = jj,nc = jd+jo;
314     for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
315     for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
316     *(++ii) = *(hdi++) + *(hoi++);
317     ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);CHKERRQ(ierr);
318   }
319   for (; cum<dr; cum++) *(++ii) = nnz;
320   if (reuse != MAT_REUSE_MATRIX) {
321     Mat_SeqAIJ* a;
322 
323     ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr);
324     ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr);
325     /* hack SeqAIJ */
326     a          = (Mat_SeqAIJ*)(lA->data);
327     a->free_a  = PETSC_TRUE;
328     a->free_ij = PETSC_TRUE;
329     ierr = MatDestroy(&lA);CHKERRQ(ierr);
330   }
331   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
332   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
333   if (reuse == MAT_INPLACE_MATRIX) {
334     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
335   }
336   PetscFunctionReturn(0);
337 }
338 
339 #undef __FUNCT__
340 #define __FUNCT__ "MatConvert_AIJ_HYPRE"
341 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
342 {
343   Mat_HYPRE      *hB;
344   MPI_Comm       comm = PetscObjectComm((PetscObject)A);
345   PetscErrorCode ierr;
346 
347   PetscFunctionBegin;
348   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
349   if (reuse == MAT_REUSE_MATRIX) {
350     /* always destroy the old matrix and create a new memory;
351        hope this does not churn the memory too much. The problem
352        is I do not know if it is possible to put the matrix back to
353        its initial state so that we can directly copy the values
354        the second time through. */
355     hB = (Mat_HYPRE*)((*B)->data);
356     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
357   } else {
358     ierr = MatCreate(comm,B);CHKERRQ(ierr);
359     ierr = MatSetType(*B,MATHYPRE);CHKERRQ(ierr);
360     ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
361     ierr = MatSetUp(*B);CHKERRQ(ierr);
362     hB   = (Mat_HYPRE*)((*B)->data);
363   }
364   ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr);
365   ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr);
366   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
367   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
368   PetscFunctionReturn(0);
369 }
370 
371 #undef __FUNCT__
372 #define __FUNCT__ "MatConvert_HYPRE_AIJ"
373 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
374 {
375   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
376   hypre_ParCSRMatrix *parcsr;
377   hypre_CSRMatrix    *hdiag,*hoffd;
378   MPI_Comm           comm;
379   PetscScalar        *da,*oa,*aptr;
380   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
381   PetscInt           i,dnnz,onnz,m,n;
382   HYPRE_Int          type;
383   PetscMPIInt        size;
384   PetscErrorCode     ierr;
385 
386   PetscFunctionBegin;
387   comm = PetscObjectComm((PetscObject)A);
388   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
389   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
390   if (reuse == MAT_REUSE_MATRIX) {
391     PetscBool ismpiaij,isseqaij;
392     ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
393     ierr = PetscObjectTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
394     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
395   }
396   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
397 
398   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
399   hdiag = hypre_ParCSRMatrixDiag(parcsr);
400   hoffd = hypre_ParCSRMatrixOffd(parcsr);
401   m     = hypre_CSRMatrixNumRows(hdiag);
402   n     = hypre_CSRMatrixNumCols(hdiag);
403   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
404   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
405   if (reuse == MAT_INITIAL_MATRIX) {
406     ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr);
407     ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr);
408     ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr);
409   } else if (reuse == MAT_REUSE_MATRIX) {
410     PetscInt  nr;
411     PetscBool done;
412     if (size > 1) {
413       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
414 
415       ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
416       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);
417       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);
418       ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr);
419     } else {
420       ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
421       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
422       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);
423       ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr);
424     }
425   } else { /* MAT_INPLACE_MATRIX */
426     dii = hypre_CSRMatrixI(hdiag);
427     djj = hypre_CSRMatrixJ(hdiag);
428     da  = hypre_CSRMatrixData(hdiag);
429   }
430   ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
431   ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr);
432   ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr);
433   iptr = djj;
434   aptr = da;
435   for (i=0; i<m; i++) {
436     PetscInt nc = dii[i+1]-dii[i];
437     ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
438     iptr += nc;
439     aptr += nc;
440   }
441   if (size > 1) {
442     PetscInt *offdj,*coffd;
443 
444     if (reuse == MAT_INITIAL_MATRIX) {
445       ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr);
446       ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr);
447       ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr);
448     } else if (reuse == MAT_REUSE_MATRIX) {
449       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
450       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
451       PetscBool  done;
452 
453       ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr);
454       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);
455       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);
456       ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr);
457     } else { /* MAT_INPLACE_MATRIX */
458       oii = hypre_CSRMatrixI(hoffd);
459       ojj = hypre_CSRMatrixJ(hoffd);
460       oa  = hypre_CSRMatrixData(hoffd);
461     }
462     ierr  = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
463     offdj = hypre_CSRMatrixJ(hoffd);
464     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
465     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
466     ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr);
467     iptr = ojj;
468     aptr = oa;
469     for (i=0; i<m; i++) {
470        PetscInt nc = oii[i+1]-oii[i];
471        ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
472        iptr += nc;
473        aptr += nc;
474     }
475     if (reuse == MAT_INITIAL_MATRIX) {
476       Mat_MPIAIJ *b;
477       Mat_SeqAIJ *d,*o;
478 
479       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);CHKERRQ(ierr);
480       /* hack MPIAIJ */
481       b          = (Mat_MPIAIJ*)((*B)->data);
482       d          = (Mat_SeqAIJ*)b->A->data;
483       o          = (Mat_SeqAIJ*)b->B->data;
484       d->free_a  = PETSC_TRUE;
485       d->free_ij = PETSC_TRUE;
486       o->free_a  = PETSC_TRUE;
487       o->free_ij = PETSC_TRUE;
488     } else if (reuse == MAT_INPLACE_MATRIX) {
489       Mat T;
490       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);CHKERRQ(ierr);
491       hypre_CSRMatrixI(hdiag)    = NULL;
492       hypre_CSRMatrixJ(hdiag)    = NULL;
493       hypre_CSRMatrixData(hdiag) = NULL;
494       hypre_CSRMatrixI(hoffd)    = NULL;
495       hypre_CSRMatrixJ(hoffd)    = NULL;
496       hypre_CSRMatrixData(hoffd) = NULL;
497       ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr);
498     }
499   } else {
500     oii  = NULL;
501     ojj  = NULL;
502     oa   = NULL;
503     if (reuse == MAT_INITIAL_MATRIX) {
504       Mat_SeqAIJ* b;
505       ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr);
506       /* hack SeqAIJ */
507       b          = (Mat_SeqAIJ*)((*B)->data);
508       b->free_a  = PETSC_TRUE;
509       b->free_ij = PETSC_TRUE;
510     } else if (reuse == MAT_INPLACE_MATRIX) {
511       Mat T;
512       ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);CHKERRQ(ierr);
513       hypre_CSRMatrixI(hdiag)    = NULL;
514       hypre_CSRMatrixJ(hdiag)    = NULL;
515       hypre_CSRMatrixData(hdiag) = NULL;
516       ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr);
517     }
518   }
519 
520   /* we have to use hypre_Tfree to free the arrays */
521   if (reuse == MAT_INPLACE_MATRIX) {
522     void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
523     const char *names[6] = {"_hypre_csr_dii",
524                             "_hypre_csr_djj",
525                             "_hypre_csr_da",
526                             "_hypre_csr_oii",
527                             "_hypre_csr_ojj",
528                             "_hypre_csr_oa"};
529     for (i=0; i<6; i++) {
530       PetscContainer c;
531 
532       ierr = PetscContainerCreate(comm,&c);CHKERRQ(ierr);
533       ierr = PetscContainerSetPointer(c,ptrs[i]);CHKERRQ(ierr);
534       ierr = PetscContainerSetUserDestroy(c,hypre_array_destroy);CHKERRQ(ierr);
535       ierr = PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);CHKERRQ(ierr);
536       ierr = PetscContainerDestroy(&c);CHKERRQ(ierr);
537     }
538   }
539   PetscFunctionReturn(0);
540 }
541 
542 #undef __FUNCT__
543 #define __FUNCT__ "MatPtAP_AIJ_HYPRE"
544 static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
545 {
546   Mat_HYPRE          *hP = (Mat_HYPRE*)P->data;
547   hypre_ParCSRMatrix *tA,*Pparcsr,*ptapparcsr;
548   hypre_CSRMatrix    *hdiag,*hoffd;
549   Mat_SeqAIJ         *diag,*offd;
550   PetscInt           *garray,noffd,dnnz,onnz,*row_starts,*col_starts;
551   HYPRE_Int          type,P_owns_col_starts;
552   PetscBool          ismpiaij,isseqaij;
553   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
554   char               mtype[256];
555   PetscErrorCode     ierr;
556 
557   PetscFunctionBegin;
558   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
559   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
560   if (scall == MAT_REUSE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_REUSE_MATRIX");
561   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
562   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
563   if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type);
564 
565   /* It looks like we don't need to have the diagonal entries
566      ordered first in the rows of the diagonal part
567      for boomerAMGBuildCoarseOperator to work */
568   if (ismpiaij) {
569     Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);
570 
571     diag   = (Mat_SeqAIJ*)a->A->data;
572     offd   = (Mat_SeqAIJ*)a->B->data;
573     garray = a->garray;
574     noffd  = a->B->cmap->N;
575     dnnz   = diag->nz;
576     onnz   = offd->nz;
577   } else {
578     diag    = (Mat_SeqAIJ*)A->data;
579     offd    = NULL;
580     garray  = NULL;
581     noffd   = 0;
582     dnnz    = diag->nz;
583     onnz    = 0;
584   }
585 
586   /* create a temporary ParCSR */
587   if (HYPRE_AssumedPartitionCheck()) {
588    PetscMPIInt myid;
589 
590    ierr       = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr);
591    row_starts = A->rmap->range + myid;
592    col_starts = A->cmap->range + myid;
593   } else {
594    row_starts = A->rmap->range;
595    col_starts = A->cmap->range;
596   }
597   tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,row_starts,col_starts,noffd,dnnz,onnz);
598   hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
599   hypre_ParCSRMatrixSetColStartsOwner(tA,0);
600 
601   /* set diagonal part */
602   hdiag = hypre_ParCSRMatrixDiag(tA);
603   hypre_CSRMatrixI(hdiag)           = diag->i;
604   hypre_CSRMatrixJ(hdiag)           = diag->j;
605   hypre_CSRMatrixData(hdiag)        = diag->a;
606   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
607   hypre_CSRMatrixSetRownnz(hdiag);
608   hypre_CSRMatrixSetDataOwner(hdiag,0);
609 
610   /* set offdiagonal part */
611   hoffd = hypre_ParCSRMatrixOffd(tA);
612   if (offd) {
613     hypre_CSRMatrixI(hoffd)           = offd->i;
614     hypre_CSRMatrixJ(hoffd)           = offd->j;
615     hypre_CSRMatrixData(hoffd)        = offd->a;
616     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
617     hypre_CSRMatrixSetRownnz(hoffd);
618     hypre_CSRMatrixSetDataOwner(hoffd,0);
619     hypre_ParCSRMatrixSetNumNonzeros(tA);
620     hypre_ParCSRMatrixColMapOffd(tA) = garray;
621   }
622 
623   /* call RAP from BoomerAMG */
624   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
625      from Pparcsr (even if it does not own them)! */
626   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
627   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(Pparcsr);
628   ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
629   PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(Pparcsr,tA,Pparcsr,&ptapparcsr));
630   ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
631   hypre_ParCSRMatrixSetRowStartsOwner(ptapparcsr,0);
632   hypre_ParCSRMatrixSetColStartsOwner(ptapparcsr,0);
633   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(Pparcsr, 1);
634 
635   /* set pointers to NULL before destroying tA */
636   hypre_CSRMatrixI(hdiag)          = NULL;
637   hypre_CSRMatrixJ(hdiag)          = NULL;
638   hypre_CSRMatrixData(hdiag)       = NULL;
639   hypre_CSRMatrixI(hoffd)          = NULL;
640   hypre_CSRMatrixJ(hoffd)          = NULL;
641   hypre_CSRMatrixData(hoffd)       = NULL;
642   hypre_ParCSRMatrixColMapOffd(tA) = NULL;
643   hypre_ParCSRMatrixDestroy(tA);
644 
645   /* create C depending on mtype */
646   sprintf(mtype,MATAIJ);
647   ierr = PetscOptionsGetString(((PetscObject)A)->options,((PetscObject)A)->prefix,"-matptap_hypre_outtype",mtype,256,NULL);CHKERRQ(ierr);
648   ierr = MatCreateFromParCSR(ptapparcsr,mtype,PETSC_OWN_POINTER,C);CHKERRQ(ierr);
649   PetscFunctionReturn(0);
650 }
651 
652 #undef __FUNCT__
653 #define __FUNCT__ "MatPtAP_HYPRE_HYPRE"
654 static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
655 {
656   hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr;
657   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data, *hP = (Mat_HYPRE*)P->data;
658   HYPRE_Int          type,P_owns_col_starts;
659   PetscErrorCode     ierr;
660 
661   PetscFunctionBegin;
662   if (scall == MAT_REUSE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Unsupported MAT_REUSE_MATRIX");
663   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
664   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
665   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
666   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
667   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
668   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
669 
670   /* call RAP from BoomerAMG */
671   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
672      from Pparcsr (even if it does not own them)! */
673   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(Pparcsr);
674   ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
675   PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr));
676   ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
677   hypre_ParCSRMatrixSetRowStartsOwner(ptapparcsr,0);
678   hypre_ParCSRMatrixSetColStartsOwner(ptapparcsr,0);
679   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(Pparcsr, 1);
680 
681   /* create MatHYPRE */
682   ierr = MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,C);CHKERRQ(ierr);
683   PetscFunctionReturn(0);
684 }
685 
686 #undef __FUNCT__
687 #define __FUNCT__ "MatMultTranspose_HYPRE"
688 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
689 {
690   PetscErrorCode ierr;
691 
692   PetscFunctionBegin;
693   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);CHKERRQ(ierr);
694   PetscFunctionReturn(0);
695 }
696 
697 #undef __FUNCT__
698 #define __FUNCT__ "MatMult_HYPRE"
699 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
700 {
701   PetscErrorCode ierr;
702 
703   PetscFunctionBegin;
704   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);CHKERRQ(ierr);
705   PetscFunctionReturn(0);
706 }
707 
708 #undef __FUNCT__
709 #define __FUNCT__ "MatHYPRE_MultKernel_Private"
710 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
711 {
712   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
713   hypre_ParCSRMatrix *parcsr;
714   hypre_ParVector    *hx,*hy;
715   PetscScalar        *ax,*ay,*sax,*say;
716   PetscErrorCode     ierr;
717 
718   PetscFunctionBegin;
719   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
720   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
721   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
722   ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
723   ierr = VecGetArray(y,&ay);CHKERRQ(ierr);
724   if (trans) {
725     VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
726     VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
727     hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
728     VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
729     VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
730   } else {
731     VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
732     VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
733     hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
734     VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
735     VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
736   }
737   ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
738   ierr = VecRestoreArray(y,&ay);CHKERRQ(ierr);
739   PetscFunctionReturn(0);
740 }
741 
742 #undef __FUNCT__
743 #define __FUNCT__ "MatDestroy_HYPRE"
744 static PetscErrorCode MatDestroy_HYPRE(Mat A)
745 {
746   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
747   PetscErrorCode ierr;
748 
749   PetscFunctionBegin;
750   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
751   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
752   if (hA->ij) {
753     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
754     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
755   }
756   if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); }
757   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr);
758   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr);
759   ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);CHKERRQ(ierr);
760   ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);CHKERRQ(ierr);
761   ierr = PetscFree(A->data);CHKERRQ(ierr);
762   PetscFunctionReturn(0);
763 }
764 
765 #undef __FUNCT__
766 #define __FUNCT__ "MatSetUp_HYPRE"
767 static PetscErrorCode MatSetUp_HYPRE(Mat A)
768 {
769   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
770   Vec                x,b;
771   PetscErrorCode     ierr;
772 
773   PetscFunctionBegin;
774   if (hA->x) PetscFunctionReturn(0);
775   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
776   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
777   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr);
778   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr);
779   ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr);
780   ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr);
781   ierr = VecDestroy(&x);CHKERRQ(ierr);
782   ierr = VecDestroy(&b);CHKERRQ(ierr);
783   PetscFunctionReturn(0);
784 }
785 
786 #undef __FUNCT__
787 #define __FUNCT__ "MatAssemblyEnd_HYPRE"
788 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
789 {
790   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
791   PetscErrorCode     ierr;
792 
793   PetscFunctionBegin;
794   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
795   PetscFunctionReturn(0);
796 }
797 
798 /*
799    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix
800 
801    Collective
802 
803    Input Parameters:
804 +  vparcsr  - the pointer to the hypre_ParCSRMatrix
805 .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
806 -  copymode - PETSc copying options
807 
808    Output Parameter:
809 .  A  - the matrix
810 
811    Level: intermediate
812 
813 .seealso: MatHYPRE, PetscCopyMode
814 */
815 #undef __FUNCT__
816 #define __FUNCT__ "MatCreateFromParCSR"
817 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(void *vparcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
818 {
819   Mat                   T;
820   Mat_HYPRE             *hA;
821   hypre_ParCSRMatrix    *parcsr;
822   MPI_Comm              comm;
823   PetscInt              rstart,rend,cstart,cend,M,N;
824   PetscBool             isseqaij,ismpiaij,isaij,ishyp,isis;
825   PetscErrorCode        ierr;
826 
827   PetscFunctionBegin;
828   parcsr = (hypre_ParCSRMatrix *)vparcsr;
829   comm   = hypre_ParCSRMatrixComm(parcsr);
830   ierr   = PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
831   ierr   = PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
832   ierr   = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr);
833   ierr   = PetscStrcmp(mtype,MATHYPRE,&ishyp);CHKERRQ(ierr);
834   ierr   = PetscStrcmp(mtype,MATIS,&isis);CHKERRQ(ierr);
835   isaij  = (PetscBool)(isseqaij || ismpiaij || isaij);
836   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);
837   if (ishyp && copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");
838 
839   /* access ParCSRMatrix */
840   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
841   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
842   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
843   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
844   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
845   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);
846 
847   /* create PETSc matrix with MatHYPRE */
848   ierr = MatCreate(comm,&T);CHKERRQ(ierr);
849   ierr = MatSetSizes(T,rend-rstart+1,cend-cstart+1,M,N);CHKERRQ(ierr);
850   ierr = MatSetType(T,MATHYPRE);CHKERRQ(ierr);
851   ierr = MatSetUp(T);CHKERRQ(ierr);
852   hA   = (Mat_HYPRE*)(T->data);
853 
854   /* create HYPRE_IJMatrix */
855   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
856 
857   /* set ParCSR object */
858   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
859   hypre_IJMatrixObject(hA->ij) = parcsr;
860 
861   /* set assembled flag */
862   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
863   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
864   if (ishyp) {
865     /* prevent from freeing the pointer */
866     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
867     *A = T;
868     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
869     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
870   } else if (isaij) {
871     if (copymode != PETSC_OWN_POINTER) {
872       /* prevent from freeing the pointer */
873       hA->inner_free = PETSC_FALSE;
874       ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr);
875       ierr = MatDestroy(&T);CHKERRQ(ierr);
876     } else { /* AIJ return type with PETSC_OWN_POINTER */
877       ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr);
878       *A   = T;
879     }
880   } else if (isis) {
881     ierr = MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr);
882     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
883     ierr = MatDestroy(&T);CHKERRQ(ierr);
884   }
885   PetscFunctionReturn(0);
886 }
887 
888 #undef __FUNCT__
889 #define __FUNCT__ "MatCreate_HYPRE"
890 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
891 {
892   Mat_HYPRE      *hB;
893   PetscErrorCode ierr;
894 
895   PetscFunctionBegin;
896   ierr           = PetscNewLog(B,&hB);CHKERRQ(ierr);
897   hB->inner_free = PETSC_TRUE;
898 
899   B->data       = (void*)hB;
900   B->rmap->bs   = 1;
901   B->assembled  = PETSC_FALSE;
902 
903   B->ops->mult          = MatMult_HYPRE;
904   B->ops->multtranspose = MatMultTranspose_HYPRE;
905   B->ops->setup         = MatSetUp_HYPRE;
906   B->ops->destroy       = MatDestroy_HYPRE;
907   B->ops->assemblyend   = MatAssemblyEnd_HYPRE;
908   B->ops->ptap          = MatPtAP_HYPRE_HYPRE;
909 
910   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr);
911   ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr);
912   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr);
913   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr);
914   ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr);
915   ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr);
916   PetscFunctionReturn(0);
917 }
918 
919 #undef __FUNCT__
920 #define __FUNCT__ "hypre_array_destroy"
921 static PetscErrorCode hypre_array_destroy(void *ptr)
922 {
923    PetscFunctionBegin;
924    hypre_TFree(ptr);
925    PetscFunctionReturn(0);
926 }
927 
928 #if 0
929 /*
930     Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format
931 
932     This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first
933     which will corrupt the PETSc data structure if we did this. Need a work around to this problem.
934 */
935 #include <_hypre_IJ_mv.h>
936 #include <HYPRE_IJ_mv.h>
937 
938 #undef __FUNCT__
939 #define __FUNCT__ "MatHYPRE_IJMatrixLink"
940 PetscErrorCode MatHYPRE_IJMatrixLink(Mat A, HYPRE_IJMatrix *ij)
941 {
942   PetscErrorCode        ierr;
943   PetscInt              rstart,rend,cstart,cend;
944   PetscBool             flg;
945   hypre_AuxParCSRMatrix *aux_matrix;
946 
947   PetscFunctionBegin;
948   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
949   PetscValidType(A,1);
950   PetscValidPointer(ij,2);
951   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
952   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices");
953   ierr = MatSetUp(A);CHKERRQ(ierr);
954 
955   ierr   = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
956   rstart = A->rmap->rstart;
957   rend   = A->rmap->rend;
958   cstart = A->cmap->rstart;
959   cend   = A->cmap->rend;
960   PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij));
961   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR));
962 
963   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(*ij));
964   PetscStackCall("hypre_IJMatrixTranslator",aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij));
965 
966   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
967 
968   /* this is the Hack part where we monkey directly with the hypre datastructures */
969 
970   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(*ij));
971   ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
972   PetscFunctionReturn(0);
973 }
974 #endif
975