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