xref: /petsc/src/mat/impls/hypre/mhypre.c (revision a033916d1bf723f6259b6e1f5089433bc451d4ab) 
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   comm = PetscObjectComm((PetscObject)A);
239   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
240   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
241   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
242   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
243   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
244   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
245   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
246   hdiag = hypre_ParCSRMatrixDiag(hA);
247   hoffd = hypre_ParCSRMatrixOffd(hA);
248   dr    = hypre_CSRMatrixNumRows(hdiag);
249   dc    = hypre_CSRMatrixNumCols(hdiag);
250   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
251   hdi   = hypre_CSRMatrixI(hdiag);
252   hdj   = hypre_CSRMatrixJ(hdiag);
253   hdd   = hypre_CSRMatrixData(hdiag);
254   oc    = hypre_CSRMatrixNumCols(hoffd);
255   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
256   hoi   = hypre_CSRMatrixI(hoffd);
257   hoj   = hypre_CSRMatrixJ(hoffd);
258   hod   = hypre_CSRMatrixData(hoffd);
259   if (reuse != MAT_REUSE_MATRIX) {
260     PetscInt *aux;
261 
262     /* generate l2g maps for rows and cols */
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     Mat_SeqAIJ* a;
310 
311     ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr);
312     ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr);
313     /* hack SeqAIJ */
314     a          = (Mat_SeqAIJ*)(lA->data);
315     a->free_a  = PETSC_TRUE;
316     a->free_ij = PETSC_TRUE;
317     ierr = MatDestroy(&lA);CHKERRQ(ierr);
318   }
319   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
320   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
321   if (reuse == MAT_INPLACE_MATRIX) {
322     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
323   }
324   PetscFunctionReturn(0);
325 }
326 
327 #undef __FUNCT__
328 #define __FUNCT__ "MatConvert_AIJ_HYPRE"
329 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
330 {
331   Mat_HYPRE      *hB;
332   MPI_Comm       comm = PetscObjectComm((PetscObject)A);
333   PetscErrorCode ierr;
334 
335   PetscFunctionBegin;
336   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
337   if (reuse == MAT_REUSE_MATRIX) {
338     /* always destroy the old matrix and create a new memory;
339        hope this does not churn the memory too much. The problem
340        is I do not know if it is possible to put the matrix back to
341        its initial state so that we can directly copy the values
342        the second time through. */
343     hB = (Mat_HYPRE*)((*B)->data);
344     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
345   } else {
346     ierr = MatCreate(comm,B);CHKERRQ(ierr);
347     ierr = MatSetType(*B,MATHYPRE);CHKERRQ(ierr);
348     ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
349     ierr = MatSetUp(*B);CHKERRQ(ierr);
350     hB   = (Mat_HYPRE*)((*B)->data);
351   }
352   ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr);
353   ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr);
354   ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
355   ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
356   PetscFunctionReturn(0);
357 }
358 
359 #undef __FUNCT__
360 #define __FUNCT__ "MatConvert_HYPRE_AIJ"
361 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
362 {
363   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
364   hypre_ParCSRMatrix *parcsr;
365   hypre_CSRMatrix    *hdiag,*hoffd;
366   MPI_Comm           comm;
367   PetscScalar        *da,*oa,*aptr;
368   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
369   PetscInt           i,dnnz,onnz,m,n;
370   int                type;
371   PetscMPIInt        size;
372   PetscErrorCode     ierr;
373 
374   PetscFunctionBegin;
375   comm = PetscObjectComm((PetscObject)A);
376   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
377   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
378   if (reuse == MAT_REUSE_MATRIX) {
379     PetscBool ismpiaij,isseqaij;
380     ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr);
381     ierr = PetscObjectTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
382     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
383   }
384   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
385 
386   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
387   hdiag = hypre_ParCSRMatrixDiag(parcsr);
388   hoffd = hypre_ParCSRMatrixOffd(parcsr);
389   m     = hypre_CSRMatrixNumRows(hdiag);
390   n     = hypre_CSRMatrixNumCols(hdiag);
391   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
392   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
393   if (reuse != MAT_REUSE_MATRIX) {
394     ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr);
395     ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr);
396     ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr);
397   } else {
398     PetscInt  nr;
399     PetscBool done;
400     if (size > 1) {
401       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
402 
403       ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
404       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);
405       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);
406       ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr);
407     } else {
408       ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr);
409       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
410       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);
411       ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr);
412     }
413   }
414   ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
415   ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr);
416   ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr);
417   iptr = djj;
418   aptr = da;
419   for (i=0; i<m; i++) {
420     PetscInt nc = dii[i+1]-dii[i];
421     ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
422     iptr += nc;
423     aptr += nc;
424   }
425   if (size > 1) {
426     PetscInt *offdj,*coffd;
427 
428     if (reuse != MAT_REUSE_MATRIX) {
429       ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr);
430       ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr);
431       ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr);
432     } else {
433       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
434       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
435       PetscBool  done;
436 
437       ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr);
438       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);
439       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);
440       ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr);
441     }
442     ierr  = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
443     offdj = hypre_CSRMatrixJ(hoffd);
444     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
445     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
446     ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr);
447     iptr = ojj;
448     aptr = oa;
449     for (i=0; i<m; i++) {
450        PetscInt nc = oii[i+1]-oii[i];
451        ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr);
452        iptr += nc;
453        aptr += nc;
454     }
455     if (reuse != MAT_REUSE_MATRIX) {
456       Mat_MPIAIJ *b;
457       Mat_SeqAIJ *d,*o;
458       ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,
459                                             djj,da,oii,ojj,oa,B);CHKERRQ(ierr);
460       /* hack MPIAIJ */
461       b          = (Mat_MPIAIJ*)((*B)->data);
462       d          = (Mat_SeqAIJ*)b->A->data;
463       o          = (Mat_SeqAIJ*)b->B->data;
464       d->free_a  = PETSC_TRUE;
465       d->free_ij = PETSC_TRUE;
466       o->free_a  = PETSC_TRUE;
467       o->free_ij = PETSC_TRUE;
468     }
469   } else if (reuse != MAT_REUSE_MATRIX) {
470     Mat_SeqAIJ* b;
471     ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr);
472     /* hack SeqAIJ */
473     b          = (Mat_SeqAIJ*)((*B)->data);
474     b->free_a  = PETSC_TRUE;
475     b->free_ij = PETSC_TRUE;
476   }
477   if (reuse == MAT_INPLACE_MATRIX) {
478     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
479   }
480   PetscFunctionReturn(0);
481 }
482 
483 #undef __FUNCT__
484 #define __FUNCT__ "MatMultTranspose_HYPRE"
485 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
486 {
487   PetscErrorCode ierr;
488 
489   PetscFunctionBegin;
490   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);CHKERRQ(ierr);
491   PetscFunctionReturn(0);
492 }
493 
494 #undef __FUNCT__
495 #define __FUNCT__ "MatMult_HYPRE"
496 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
497 {
498   PetscErrorCode ierr;
499 
500   PetscFunctionBegin;
501   ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);CHKERRQ(ierr);
502   PetscFunctionReturn(0);
503 }
504 
505 #undef __FUNCT__
506 #define __FUNCT__ "MatHYPRE_MultKernel_Private"
507 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
508 {
509   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
510   hypre_ParCSRMatrix *parcsr;
511   hypre_ParVector    *hx,*hy;
512   PetscScalar        *ax,*ay,*sax,*say;
513   PetscErrorCode     ierr;
514 
515   PetscFunctionBegin;
516   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
517   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
518   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
519   ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
520   ierr = VecGetArray(y,&ay);CHKERRQ(ierr);
521   if (trans) {
522     HYPREReplacePointer(hA->x,ay,say);
523     HYPREReplacePointer(hA->b,ax,sax);
524     hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
525     HYPREReplacePointer(hA->x,say,ay);
526     HYPREReplacePointer(hA->b,sax,ax);
527   } else {
528     HYPREReplacePointer(hA->x,ax,sax);
529     HYPREReplacePointer(hA->b,ay,say);
530     hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
531     HYPREReplacePointer(hA->x,sax,ax);
532     HYPREReplacePointer(hA->b,say,ay);
533   }
534   ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr);
535   ierr = VecRestoreArray(y,&ay);CHKERRQ(ierr);
536   PetscFunctionReturn(0);
537 }
538 
539 #undef __FUNCT__
540 #define __FUNCT__ "MatDestroy_HYPRE"
541 static PetscErrorCode MatDestroy_HYPRE(Mat A)
542 {
543   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
544   PetscErrorCode ierr;
545 
546   PetscFunctionBegin;
547   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
548   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
549   if (hA->ij) {
550     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
551     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
552   }
553   if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); }
554   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr);
555   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr);
556   ierr = PetscFree(A->data);CHKERRQ(ierr);
557   PetscFunctionReturn(0);
558 }
559 
560 #undef __FUNCT__
561 #define __FUNCT__ "MatSetUp_HYPRE"
562 static PetscErrorCode MatSetUp_HYPRE(Mat A)
563 {
564   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
565   Vec                x,b;
566   PetscErrorCode     ierr;
567 
568   PetscFunctionBegin;
569   if (hA->x) PetscFunctionReturn(0);
570   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
571   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
572   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr);
573   ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr);
574   ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr);
575   ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr);
576   ierr = VecDestroy(&x);CHKERRQ(ierr);
577   ierr = VecDestroy(&b);CHKERRQ(ierr);
578   PetscFunctionReturn(0);
579 }
580 
581 #undef __FUNCT__
582 #define __FUNCT__ "MatAssemblyEnd_HYPRE"
583 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
584 {
585   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
586   PetscErrorCode     ierr;
587 
588   PetscFunctionBegin;
589   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
590   PetscFunctionReturn(0);
591 }
592 
593 #undef __FUNCT__
594 #define __FUNCT__ "MatHYPRECreateFromParCSR"
595 PETSC_EXTERN PetscErrorCode MatHYPRECreateFromParCSR(void *vparcsr, PetscCopyMode copymode, Mat* A)
596 {
597   Mat_HYPRE             *hA;
598   hypre_ParCSRMatrix    *parcsr;
599   MPI_Comm              comm;
600   PetscInt              rstart,rend,cstart,cend,M,N;
601   PetscErrorCode        ierr;
602 
603   PetscFunctionBegin;
604   parcsr = (hypre_ParCSRMatrix *)vparcsr;
605   comm   = hypre_ParCSRMatrixComm(parcsr);
606   if (copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");
607 
608   /* access ParCSRMatrix */
609   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
610   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
611   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
612   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
613   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
614   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);
615 
616   /* create PETSc matrix with MatHYPRE */
617   ierr = MatCreate(comm,A);CHKERRQ(ierr);
618   ierr = MatSetSizes(*A,rend-rstart+1,cend-cstart+1,M,N);CHKERRQ(ierr);
619   ierr = MatSetType(*A,MATHYPRE);CHKERRQ(ierr);
620   ierr = MatSetUp(*A);CHKERRQ(ierr);
621   hA   = (Mat_HYPRE*)((*A)->data);
622 
623   /* create HYPRE_IJMatrix */
624   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
625 
626   /* set ParCSR object */
627   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
628   hypre_IJMatrixObject(hA->ij) = parcsr;
629 
630   /* set assembled flag */
631   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
632   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
633 
634   /* prevent from freeing the pointer */
635   if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
636 
637   ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
638   ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
639   PetscFunctionReturn(0);
640 }
641 
642 #undef __FUNCT__
643 #define __FUNCT__ "MatCreate_HYPRE"
644 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
645 {
646   Mat_HYPRE      *hB;
647   PetscErrorCode ierr;
648 
649   PetscFunctionBegin;
650   ierr           = PetscNewLog(B,&hB);CHKERRQ(ierr);
651   hB->inner_free = PETSC_TRUE;
652 
653   B->data       = (void*)hB;
654   B->rmap->bs   = 1;
655   B->assembled  = PETSC_FALSE;
656 
657   B->ops->mult          = MatMult_HYPRE;
658   B->ops->multtranspose = MatMultTranspose_HYPRE;
659   B->ops->setup         = MatSetUp_HYPRE;
660   B->ops->destroy       = MatDestroy_HYPRE;
661   B->ops->assemblyend   = MatAssemblyEnd_HYPRE;
662 
663   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr);
664   ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr);
665   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr);
666   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr);
667   PetscFunctionReturn(0);
668 }
669 
670 #if 0
671 /*
672     Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format
673 
674     This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first
675     which will corrupt the PETSc data structure if we did this. Need a work around to this problem.
676 */
677 #include <_hypre_IJ_mv.h>
678 #include <HYPRE_IJ_mv.h>
679 
680 #undef __FUNCT__
681 #define __FUNCT__ "MatHYPRE_IJMatrixLink"
682 PetscErrorCode MatHYPRE_IJMatrixLink(Mat A, HYPRE_IJMatrix *ij)
683 {
684   PetscErrorCode        ierr;
685   PetscInt              rstart,rend,cstart,cend;
686   PetscBool             flg;
687   hypre_AuxParCSRMatrix *aux_matrix;
688 
689   PetscFunctionBegin;
690   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
691   PetscValidType(A,1);
692   PetscValidPointer(ij,2);
693   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
694   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices");
695   ierr = MatSetUp(A);CHKERRQ(ierr);
696 
697   ierr   = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
698   rstart = A->rmap->rstart;
699   rend   = A->rmap->rend;
700   cstart = A->cmap->rstart;
701   cend   = A->cmap->rend;
702   PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij));
703   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR));
704 
705   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(*ij));
706   PetscStackCall("hypre_IJMatrixTranslator",aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij));
707 
708   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
709 
710   /* this is the Hack part where we monkey directly with the hypre datastructures */
711 
712   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(*ij));
713   ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr);
714   PetscFunctionReturn(0);
715 }
716 #endif
717