xref: /petsc/src/mat/impls/hypre/mhypre.c (revision 2df223495fa2a5750e9e2d07b4eba9c91a3b5014) 
1 
2 /*
3     Creates hypre ijmatrix from PETSc matrix
4 */
5 #include <petscsys.h>
6 #include <petsc/private/matimpl.h>
7 #include <../src/mat/impls/hypre/mhypre.h>
8 #include <../src/mat/impls/aij/mpi/mpiaij.h>
9 #include <HYPRE_struct_mv.h>
10 #include <HYPRE_struct_ls.h>
11 #include <_hypre_struct_mv.h>
12 
13 static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
14 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
15 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
16 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
17 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,Vec,Vec,PetscBool);
18 
19 #undef __FUNCT__
20 #define __FUNCT__ "MatHYPRE_IJMatrixPreallocate"
21 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
22 {
23   PetscErrorCode ierr;
24   PetscInt       i,n_d,n_o;
25   const PetscInt *ia_d,*ia_o;
26   PetscBool      done_d=PETSC_FALSE,done_o=PETSC_FALSE;
27   PetscInt       *nnz_d=NULL,*nnz_o=NULL;
28 
29   PetscFunctionBegin;
30   if (A_d) { /* determine number of nonzero entries in local diagonal part */
31     ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);CHKERRQ(ierr);
32     if (done_d) {
33       ierr = PetscMalloc1(n_d,&nnz_d);CHKERRQ(ierr);
34       for (i=0; i<n_d; i++) {
35         nnz_d[i] = ia_d[i+1] - ia_d[i];
36       }
37     }
38     ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);CHKERRQ(ierr);
39   }
40   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
41     ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
42     if (done_o) {
43       ierr = PetscMalloc1(n_o,&nnz_o);CHKERRQ(ierr);
44       for (i=0; i<n_o; i++) {
45         nnz_o[i] = ia_o[i+1] - ia_o[i];
46       }
47     }
48     ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr);
49   }
50   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
51     if (!done_o) { /* only diagonal part */
52       ierr = PetscMalloc1(n_d,&nnz_o);CHKERRQ(ierr);
53       for (i=0; i<n_d; i++) {
54         nnz_o[i] = 0;
55       }
56     }
57     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
58     ierr = PetscFree(nnz_d);CHKERRQ(ierr);
59     ierr = PetscFree(nnz_o);CHKERRQ(ierr);
60   }
61   PetscFunctionReturn(0);
62 }
63 
64 #undef __FUNCT__
65 #define __FUNCT__ "MatHYPRE_CreateFromMat"
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   = MatSetUp(A);CHKERRQ(ierr);
73   rstart = A->rmap->rstart;
74   rend   = A->rmap->rend;
75   cstart = A->cmap->rstart;
76   cend   = A->cmap->rend;
77   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
78   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
79   {
80     PetscBool      same;
81     Mat            A_d,A_o;
82     const PetscInt *colmap;
83     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr);
84     if (same) {
85       ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
86       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
87       PetscFunctionReturn(0);
88     }
89     ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr);
90     if (same) {
91       ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr);
92       ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr);
93       PetscFunctionReturn(0);
94     }
95     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr);
96     if (same) {
97       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
98       PetscFunctionReturn(0);
99     }
100     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr);
101     if (same) {
102       ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr);
103       PetscFunctionReturn(0);
104     }
105   }
106   PetscFunctionReturn(0);
107 }
108 
109 #undef __FUNCT__
110 #define __FUNCT__ "MatHYPRE_IJMatrixCopy"
111 static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
112 {
113   PetscErrorCode    ierr;
114   PetscInt          i,rstart,rend,ncols,nr,nc;
115   const PetscScalar *values;
116   const PetscInt    *cols;
117   PetscBool         flg;
118 
119   PetscFunctionBegin;
120   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
121   ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr);
122   if (flg && nr == nc) {
123     ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr);
124     PetscFunctionReturn(0);
125   }
126   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
127   if (flg) {
128     ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr);
129     PetscFunctionReturn(0);
130   }
131 
132   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
133   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
134   for (i=rstart; i<rend; i++) {
135     ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
136     PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values));
137     ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr);
138   }
139   PetscFunctionReturn(0);
140 }
141 
142 #undef __FUNCT__
143 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_SeqAIJ"
144 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
145 {
146   PetscErrorCode        ierr;
147   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ*)A->data;
148   int                   type;
149   hypre_ParCSRMatrix    *par_matrix;
150   hypre_AuxParCSRMatrix *aux_matrix;
151   hypre_CSRMatrix       *hdiag;
152 
153   PetscFunctionBegin;
154   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
155   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
156   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
157   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
158   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
159   /*
160        this is the Hack part where we monkey directly with the hypre datastructures
161   */
162   ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
163   ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
164   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
165 
166   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
167   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
168   PetscFunctionReturn(0);
169 }
170 
171 #undef __FUNCT__
172 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_MPIAIJ"
173 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
174 {
175   PetscErrorCode        ierr;
176   Mat_MPIAIJ            *pA = (Mat_MPIAIJ*)A->data;
177   Mat_SeqAIJ            *pdiag,*poffd;
178   PetscInt              i,*garray = pA->garray,*jj,cstart,*pjj;
179   int                   type;
180   hypre_ParCSRMatrix    *par_matrix;
181   hypre_AuxParCSRMatrix *aux_matrix;
182   hypre_CSRMatrix       *hdiag,*hoffd;
183 
184   PetscFunctionBegin;
185   pdiag = (Mat_SeqAIJ*) pA->A->data;
186   poffd = (Mat_SeqAIJ*) pA->B->data;
187   /* cstart is only valid for square MPIAIJ layed out in the usual way */
188   ierr = MatGetOwnershipRange(A,&cstart,NULL);CHKERRQ(ierr);
189 
190   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
191   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
192   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
193   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
194   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
195   hoffd = hypre_ParCSRMatrixOffd(par_matrix);
196 
197   /*
198        this is the Hack part where we monkey directly with the hypre datastructures
199   */
200   ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
201   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
202   jj  = (PetscInt*)hdiag->j;
203   pjj = (PetscInt*)pdiag->j;
204   for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
205   ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
206   ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
207   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
208      If we hacked a hypre a bit more we might be able to avoid this step */
209   jj  = (PetscInt*) hoffd->j;
210   pjj = (PetscInt*) poffd->j;
211   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
212   ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));
213 
214   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
215   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
216   PetscFunctionReturn(0);
217 }
218 
219 #undef __FUNCT__
220 #define __FUNCT__ "MatConvert_HYPRE_IS"
221 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
222 {
223   Mat_HYPRE*             mhA = (Mat_HYPRE*)(A->data);
224   Mat                    lA;
225   ISLocalToGlobalMapping rl2g,cl2g;
226   IS                     is;
227   hypre_ParCSRMatrix     *hA;
228   hypre_CSRMatrix        *hdiag,*hoffd;
229   MPI_Comm               comm;
230   PetscScalar            *hdd,*hod,*aa,*data;
231   PetscInt               *col_map_offd,*hdi,*hdj,*hoi,*hoj;
232   PetscInt               *ii,*jj,*iptr,*jptr;
233   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
234   int                    type;
235   PetscErrorCode         ierr;
236 
237   PetscFunctionBegin;
238   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
239   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
240   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
241   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
242   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
243   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
244   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
245   hdiag = hypre_ParCSRMatrixDiag(hA);
246   hoffd = hypre_ParCSRMatrixOffd(hA);
247   dr    = hypre_CSRMatrixNumRows(hdiag);
248   dc    = hypre_CSRMatrixNumCols(hdiag);
249   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
250   hdi   = hypre_CSRMatrixI(hdiag);
251   hdj   = hypre_CSRMatrixJ(hdiag);
252   hdd   = hypre_CSRMatrixData(hdiag);
253   oc    = hypre_CSRMatrixNumCols(hoffd);
254   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
255   hoi   = hypre_CSRMatrixI(hoffd);
256   hoj   = hypre_CSRMatrixJ(hoffd);
257   hod   = hypre_CSRMatrixData(hoffd);
258   if (reuse != MAT_REUSE_MATRIX) {
259     PetscInt *aux;
260 
261     /* generate l2g maps for rows and cols */
262     comm = PetscObjectComm((PetscObject)A);
263     ierr = ISCreateStride(comm,dr,str,1,&is);CHKERRQ(ierr);
264     ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr);
265     ierr = ISDestroy(&is);CHKERRQ(ierr);
266     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
267     ierr = PetscMalloc1(dc+oc,&aux);CHKERRQ(ierr);
268     for (i=0; i<dc; i++) aux[i] = i+stc;
269     for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
270     ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);CHKERRQ(ierr);
271     ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr);
272     ierr = ISDestroy(&is);CHKERRQ(ierr);
273     /* create MATIS object */
274     ierr = MatCreate(comm,B);CHKERRQ(ierr);
275     ierr = MatSetSizes(*B,dr,dc,M,N);CHKERRQ(ierr);
276     ierr = MatSetType(*B,MATIS);CHKERRQ(ierr);
277     ierr = MatSetLocalToGlobalMapping(*B,rl2g,cl2g);CHKERRQ(ierr);
278     ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr);
279     ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr);
280 
281     /* allocate CSR for local matrix */
282     ierr = PetscMalloc1(dr+1,&iptr);CHKERRQ(ierr);
283     ierr = PetscMalloc1(nnz,&jptr);CHKERRQ(ierr);
284     ierr = PetscMalloc1(nnz,&data);CHKERRQ(ierr);
285   } else {
286     PetscInt  nr;
287     PetscBool done;
288     ierr = MatISGetLocalMat(*B,&lA);CHKERRQ(ierr);
289     ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);CHKERRQ(ierr);
290     if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
291     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);
292     ierr = MatSeqAIJGetArray(lA,&data);CHKERRQ(ierr);
293   }
294   /* merge local matrices */
295   ii   = iptr;
296   jj   = jptr;
297   aa   = data;
298   *ii  = *(hdi++) + *(hoi++);
299   for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
300     PetscScalar *aold = aa;
301     PetscInt    *jold = jj,nc = jd+jo;
302     for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
303     for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
304     *(++ii) = *(hdi++) + *(hoi++);
305     ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);CHKERRQ(ierr);
306   }
307   for (; cum<dr; cum++) *(++ii) = nnz;
308   if (reuse != MAT_REUSE_MATRIX) {
309     ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr);
310     ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr);
311     ierr = MatDestroy(&lA);CHKERRQ(ierr);
312   }
313   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
314   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
315   if (reuse == MAT_INPLACE_MATRIX) {
316     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
317   }
318   PetscFunctionReturn(0);
319 }
320 
321 #undef __FUNCT__
322 #define __FUNCT__ "MatConvert_AIJ_HYPRE"
323 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
324 {
325   Mat_HYPRE      *hB;
326   MPI_Comm       comm = PetscObjectComm((PetscObject)A);
327   PetscErrorCode ierr;
328 
329   PetscFunctionBegin;
330   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
331   if (reuse == MAT_REUSE_MATRIX) {
332     /* always destroy the old matrix and create a new memory;
333        hope this does not churn the memory too much. The problem
334        is I do not know if it is possible to put the matrix back to
335        its initial state so that we can directly copy the values
336        the second time through. */
337     hB = (Mat_HYPRE*)((*B)->data);
338     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
339   } else {
340     ierr = MatCreate(comm,B);CHKERRQ(ierr);
341     ierr = MatSetType(*B,MATHYPRE);CHKERRQ(ierr);
342     ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
343     ierr = MatSetUp(*B);CHKERRQ(ierr);
344     hB   = (Mat_HYPRE*)((*B)->data);
345   }
346   ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr);
347   ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr);
348   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
349   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
350   PetscFunctionReturn(0);
351 }
352 
353 #undef __FUNCT__
354 #define __FUNCT__ "MatConvert_HYPRE_AIJ"
355 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
356 {
357   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
358   hypre_ParCSRMatrix *parcsr;
359   hypre_CSRMatrix    *hdiag,*hoffd;
360   MPI_Comm           comm;
361   PetscScalar        *da,*oa,*aptr;
362   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
363   PetscInt           i,dnnz,onnz,m,n;
364   int                type;
365   PetscMPIInt        size;
366   PetscErrorCode     ierr;
367 
368   PetscFunctionBegin;
369   comm = PetscObjectComm((PetscObject)A);
370   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
371   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
372   if (reuse == MAT_REUSE_MATRIX) {
373     PetscBool ismpiaij,isseqaij;
374     ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
375     ierr = PetscObjectTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
376     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
377   }
378   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
379 
380   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
381   hdiag = hypre_ParCSRMatrixDiag(parcsr);
382   hoffd = hypre_ParCSRMatrixOffd(parcsr);
383   m     = hypre_CSRMatrixNumRows(hdiag);
384   n     = hypre_CSRMatrixNumCols(hdiag);
385   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
386   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
387   if (reuse != MAT_REUSE_MATRIX) {
388     ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr);
389     ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr);
390     ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr);
391   } else {
392     PetscInt  nr;
393     PetscBool done;
394     if (size > 1) {
395       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
396 
397       ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
398       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows in diag part! %D != %D",nr,m);
399       if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in diag part! reuse %D hypre %D",dii[nr],dnnz);
400       ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr);
401     } else {
402       ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
403       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
404       if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros! reuse %D hypre %D",dii[nr],dnnz);
405       ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr);
406     }
407   }
408   ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
409   ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr);
410   ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr);
411   iptr = djj;
412   aptr = da;
413   for (i=0; i<m; i++) {
414     PetscInt nc = dii[i+1]-dii[i];
415     ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
416     iptr += nc;
417     aptr += nc;
418   }
419   if (size > 1) {
420     PetscInt *offdj,*coffd;
421 
422     if (reuse != MAT_REUSE_MATRIX) {
423       ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr);
424       ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr);
425       ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr);
426     } else {
427       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
428       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
429       PetscBool  done;
430 
431       ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr);
432       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);
433       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);
434       ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr);
435     }
436     ierr  = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
437     offdj = hypre_CSRMatrixJ(hoffd);
438     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
439     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
440     ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr);
441     iptr = ojj;
442     aptr = oa;
443     for (i=0; i<m; i++) {
444        PetscInt nc = oii[i+1]-oii[i];
445        ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
446        iptr += nc;
447        aptr += nc;
448     }
449     if (reuse != MAT_REUSE_MATRIX) {
450       Mat_MPIAIJ *b;
451       Mat_SeqAIJ *d,*o;
452       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,
453                                             djj,da,oii,ojj,oa,B);CHKERRQ(ierr);
454       /* hack MPIAIJ */
455       b          = (Mat_MPIAIJ*)((*B)->data);
456       d          = (Mat_SeqAIJ*)b->A->data;
457       o          = (Mat_SeqAIJ*)b->B->data;
458       d->free_a  = PETSC_TRUE;
459       d->free_ij = PETSC_TRUE;
460       o->free_a  = PETSC_TRUE;
461       o->free_ij = PETSC_TRUE;
462     }
463   } else if (reuse != MAT_REUSE_MATRIX) {
464     Mat_SeqAIJ* b;
465     ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr);
466     /* hack SeqAIJ */
467     b          = (Mat_SeqAIJ*)((*B)->data);
468     b->free_a  = PETSC_TRUE;
469     b->free_ij = PETSC_TRUE;
470   }
471   if (reuse == MAT_INPLACE_MATRIX) {
472     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
473   }
474   PetscFunctionReturn(0);
475 }
476 
477 #undef __FUNCT__
478 #define __FUNCT__ "MatMultTranspose_HYPRE"
479 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
480 {
481   PetscErrorCode ierr;
482 
483   PetscFunctionBegin;
484   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);CHKERRQ(ierr);
485   PetscFunctionReturn(0);
486 }
487 
488 #undef __FUNCT__
489 #define __FUNCT__ "MatMult_HYPRE"
490 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
491 {
492   PetscErrorCode ierr;
493 
494   PetscFunctionBegin;
495   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);CHKERRQ(ierr);
496   PetscFunctionReturn(0);
497 }
498 
499 #undef __FUNCT__
500 #define __FUNCT__ "MatHYPRE_MultKernel_Private"
501 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
502 {
503   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
504   hypre_ParCSRMatrix *parcsr;
505   hypre_ParVector    *hx,*hy;
506   PetscScalar        *ax,*ay,*sax,*say;
507   PetscErrorCode     ierr;
508 
509   PetscFunctionBegin;
510   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
511   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
512   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
513   ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
514   ierr = VecGetArray(y,&ay);CHKERRQ(ierr);
515   if (trans) {
516     HYPREReplacePointer(hA->x,ay,say);
517     HYPREReplacePointer(hA->b,ax,sax);
518     hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
519     HYPREReplacePointer(hA->x,say,ay);
520     HYPREReplacePointer(hA->b,sax,ax);
521   } else {
522     HYPREReplacePointer(hA->x,ax,sax);
523     HYPREReplacePointer(hA->b,ay,say);
524     hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
525     HYPREReplacePointer(hA->x,sax,ax);
526     HYPREReplacePointer(hA->b,say,ay);
527   }
528   ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
529   ierr = VecRestoreArray(y,&ay);CHKERRQ(ierr);
530   PetscFunctionReturn(0);
531 }
532 
533 #undef __FUNCT__
534 #define __FUNCT__ "MatDestroy_HYPRE"
535 static PetscErrorCode MatDestroy_HYPRE(Mat A)
536 {
537   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
538   PetscErrorCode ierr;
539 
540   PetscFunctionBegin;
541   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
542   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
543   if (hA->ij) {
544     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
545     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
546   }
547   if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); }
548   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr);
549   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr);
550   ierr = PetscFree(A->data);CHKERRQ(ierr);
551   PetscFunctionReturn(0);
552 }
553 
554 #undef __FUNCT__
555 #define __FUNCT__ "MatSetUp_HYPRE"
556 static PetscErrorCode MatSetUp_HYPRE(Mat A)
557 {
558   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
559   Vec                x,b;
560   PetscErrorCode     ierr;
561 
562   PetscFunctionBegin;
563   if (hA->x) PetscFunctionReturn(0);
564   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
565   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
566   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr);
567   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr);
568   ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr);
569   ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr);
570   ierr = VecDestroy(&x);CHKERRQ(ierr);
571   ierr = VecDestroy(&b);CHKERRQ(ierr);
572   PetscFunctionReturn(0);
573 }
574 
575 #undef __FUNCT__
576 #define __FUNCT__ "MatAssemblyEnd_HYPRE"
577 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
578 {
579   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
580   PetscErrorCode     ierr;
581 
582   PetscFunctionBegin;
583   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
584   PetscFunctionReturn(0);
585 }
586 
587 #undef __FUNCT__
588 #define __FUNCT__ "MatHYPRECreateFromParCSR"
589 PETSC_EXTERN PetscErrorCode MatHYPRECreateFromParCSR(void *vparcsr, PetscCopyMode copymode, Mat* A)
590 {
591   Mat_HYPRE             *hA;
592   hypre_ParCSRMatrix    *parcsr;
593   MPI_Comm              comm;
594   PetscInt              rstart,rend,cstart,cend,M,N;
595   PetscErrorCode        ierr;
596 
597   PetscFunctionBegin;
598   parcsr = (hypre_ParCSRMatrix *)vparcsr;
599   comm   = hypre_ParCSRMatrixComm(parcsr);
600   if (copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");
601 
602   /* access ParCSRMatrix */
603   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
604   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
605   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
606   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
607   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
608   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);
609 
610   /* create PETSc matrix with MatHYPRE */
611   ierr = MatCreate(comm,A);CHKERRQ(ierr);
612   ierr = MatSetSizes(*A,rend-rstart+1,cend-cstart+1,M,N);CHKERRQ(ierr);
613   ierr = MatSetType(*A,MATHYPRE);CHKERRQ(ierr);
614   ierr = MatSetUp(*A);CHKERRQ(ierr);
615   hA   = (Mat_HYPRE*)((*A)->data);
616 
617   /* create HYPRE_IJMatrix */
618   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
619 
620   /* set ParCSR object */
621   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
622   hypre_IJMatrixObject(hA->ij) = parcsr;
623 
624   /* set assembled flag */
625   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
626   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
627 
628   /* prevent from freeing the pointer */
629   if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
630 
631   ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
632   ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
633   PetscFunctionReturn(0);
634 }
635 
636 #undef __FUNCT__
637 #define __FUNCT__ "MatCreate_HYPRE"
638 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
639 {
640   Mat_HYPRE      *hB;
641   PetscErrorCode ierr;
642 
643   PetscFunctionBegin;
644   ierr           = PetscNewLog(B,&hB);CHKERRQ(ierr);
645   hB->inner_free = PETSC_TRUE;
646 
647   B->data       = (void*)hB;
648   B->rmap->bs   = 1;
649   B->assembled  = PETSC_FALSE;
650 
651   B->ops->mult          = MatMult_HYPRE;
652   B->ops->multtranspose = MatMultTranspose_HYPRE;
653   B->ops->setup         = MatSetUp_HYPRE;
654   B->ops->destroy       = MatDestroy_HYPRE;
655   B->ops->assemblyend   = MatAssemblyEnd_HYPRE;
656 
657   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr);
658   ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr);
659   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr);
660   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr);
661   PetscFunctionReturn(0);
662 }
663 
664 #if 0
665 /*
666     Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format
667 
668     This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first
669     which will corrupt the PETSc data structure if we did this. Need a work around to this problem.
670 */
671 #include <_hypre_IJ_mv.h>
672 #include <HYPRE_IJ_mv.h>
673 
674 #undef __FUNCT__
675 #define __FUNCT__ "MatHYPRE_IJMatrixLink"
676 PetscErrorCode MatHYPRE_IJMatrixLink(Mat A, HYPRE_IJMatrix *ij)
677 {
678   PetscErrorCode        ierr;
679   PetscInt              rstart,rend,cstart,cend;
680   PetscBool             flg;
681   hypre_AuxParCSRMatrix *aux_matrix;
682 
683   PetscFunctionBegin;
684   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
685   PetscValidType(A,1);
686   PetscValidPointer(ij,2);
687   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
688   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices");
689   ierr = MatSetUp(A);CHKERRQ(ierr);
690 
691   ierr   = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
692   rstart = A->rmap->rstart;
693   rend   = A->rmap->rend;
694   cstart = A->cmap->rstart;
695   cend   = A->cmap->rend;
696   PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij));
697   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR));
698 
699   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(*ij));
700   PetscStackCall("hypre_IJMatrixTranslator",aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij));
701 
702   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
703 
704   /* this is the Hack part where we monkey directly with the hypre datastructures */
705 
706   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(*ij));
707   ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
708   PetscFunctionReturn(0);
709 }
710 #endif
711