1 2 /* 3 Creates hypre ijmatrix from PETSc matrix 4 */ 5 6 #include <petscpkg_version.h> 7 #include <petsc/private/petschypre.h> 8 #include <petscmathypre.h> 9 #include <petsc/private/matimpl.h> 10 #include <../src/mat/impls/hypre/mhypre.h> 11 #include <../src/mat/impls/aij/mpi/mpiaij.h> 12 #include <../src/vec/vec/impls/hypre/vhyp.h> 13 #include <HYPRE.h> 14 #include <HYPRE_utilities.h> 15 #include <_hypre_parcsr_ls.h> 16 #include <_hypre_sstruct_ls.h> 17 18 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 19 20 static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*); 21 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix); 22 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix); 23 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix); 24 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,HYPRE_Complex,Vec,HYPRE_Complex,Vec,PetscBool); 25 static PetscErrorCode hypre_array_destroy(void*); 26 PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode ins); 27 28 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij) 29 { 30 PetscErrorCode ierr; 31 PetscInt i,n_d,n_o; 32 const PetscInt *ia_d,*ia_o; 33 PetscBool done_d=PETSC_FALSE,done_o=PETSC_FALSE; 34 HYPRE_Int *nnz_d=NULL,*nnz_o=NULL; 35 36 PetscFunctionBegin; 37 if (A_d) { /* determine number of nonzero entries in local diagonal part */ 38 ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);CHKERRQ(ierr); 39 if (done_d) { 40 ierr = PetscMalloc1(n_d,&nnz_d);CHKERRQ(ierr); 41 for (i=0; i<n_d; i++) { 42 nnz_d[i] = ia_d[i+1] - ia_d[i]; 43 } 44 } 45 ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);CHKERRQ(ierr); 46 } 47 if (A_o) { /* determine number of nonzero entries in local off-diagonal part */ 48 ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr); 49 if (done_o) { 50 ierr = PetscMalloc1(n_o,&nnz_o);CHKERRQ(ierr); 51 for (i=0; i<n_o; i++) { 52 nnz_o[i] = ia_o[i+1] - ia_o[i]; 53 } 54 } 55 ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr); 56 } 57 if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */ 58 if (!done_o) { /* only diagonal part */ 59 ierr = PetscMalloc1(n_d,&nnz_o);CHKERRQ(ierr); 60 for (i=0; i<n_d; i++) { 61 nnz_o[i] = 0; 62 } 63 } 64 #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0) 65 { /* If we don't do this, the columns of the matrix will be all zeros! */ 66 hypre_AuxParCSRMatrix *aux_matrix; 67 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 68 hypre_AuxParCSRMatrixDestroy(aux_matrix); 69 hypre_IJMatrixTranslator(ij) = NULL; 70 PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,nnz_d,nnz_o)); 71 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 72 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; 73 } 74 #else 75 PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,nnz_d,nnz_o)); 76 #endif 77 ierr = PetscFree(nnz_d);CHKERRQ(ierr); 78 ierr = PetscFree(nnz_o);CHKERRQ(ierr); 79 } 80 PetscFunctionReturn(0); 81 } 82 83 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA) 84 { 85 PetscErrorCode ierr; 86 PetscInt rstart,rend,cstart,cend; 87 88 PetscFunctionBegin; 89 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 90 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 91 rstart = A->rmap->rstart; 92 rend = A->rmap->rend; 93 cstart = A->cmap->rstart; 94 cend = A->cmap->rend; 95 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij)); 96 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 97 { 98 PetscBool same; 99 Mat A_d,A_o; 100 const PetscInt *colmap; 101 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr); 102 if (same) { 103 ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 104 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr); 108 if (same) { 109 ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 110 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr); 114 if (same) { 115 ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr); 119 if (same) { 120 ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 } 124 PetscFunctionReturn(0); 125 } 126 127 static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij) 128 { 129 PetscErrorCode ierr; 130 PetscInt i,rstart,rend,ncols,nr,nc; 131 const PetscScalar *values; 132 const PetscInt *cols; 133 PetscBool flg; 134 135 PetscFunctionBegin; 136 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 137 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 138 if (flg && nr == nc) { 139 ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 143 if (flg) { 144 ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 149 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 150 for (i=rstart; i<rend; i++) { 151 ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 152 if (ncols) { 153 HYPRE_Int nc = (HYPRE_Int)ncols; 154 155 if ((PetscInt)nc != ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",ncols,i); 156 PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,&nc,(HYPRE_BigInt *)&i,(HYPRE_BigInt *)cols,(HYPRE_Complex *)values)); 157 } 158 ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 159 } 160 PetscFunctionReturn(0); 161 } 162 163 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij) 164 { 165 PetscErrorCode ierr; 166 Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data; 167 HYPRE_Int type; 168 hypre_ParCSRMatrix *par_matrix; 169 hypre_AuxParCSRMatrix *aux_matrix; 170 hypre_CSRMatrix *hdiag; 171 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 172 173 PetscFunctionBegin; 174 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 175 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type)); 176 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 177 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix)); 178 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 179 /* 180 this is the Hack part where we monkey directly with the hypre datastructures 181 */ 182 if (sameint) { 183 ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));CHKERRQ(ierr); 184 ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));CHKERRQ(ierr); 185 } else { 186 PetscInt i; 187 188 for (i=0;i<A->rmap->n + 1;i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i]; 189 for (i=0;i<pdiag->nz;i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i]; 190 } 191 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));CHKERRQ(ierr); 192 193 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 194 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij) 199 { 200 PetscErrorCode ierr; 201 Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data; 202 Mat_SeqAIJ *pdiag,*poffd; 203 PetscInt i,*garray = pA->garray,*jj,cstart,*pjj; 204 HYPRE_Int *hjj,type; 205 hypre_ParCSRMatrix *par_matrix; 206 hypre_AuxParCSRMatrix *aux_matrix; 207 hypre_CSRMatrix *hdiag,*hoffd; 208 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 209 210 PetscFunctionBegin; 211 pdiag = (Mat_SeqAIJ*) pA->A->data; 212 poffd = (Mat_SeqAIJ*) pA->B->data; 213 /* cstart is only valid for square MPIAIJ layed out in the usual way */ 214 ierr = MatGetOwnershipRange(A,&cstart,NULL);CHKERRQ(ierr); 215 216 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 217 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type)); 218 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 219 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix)); 220 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 221 hoffd = hypre_ParCSRMatrixOffd(par_matrix); 222 223 /* 224 this is the Hack part where we monkey directly with the hypre datastructures 225 */ 226 if (sameint) { 227 ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 228 } else { 229 for (i=0; i<pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)(pdiag->i[i]); 230 } 231 /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */ 232 hjj = hdiag->j; 233 pjj = pdiag->j; 234 #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0) 235 for (i=0; i<pdiag->nz; i++) hjj[i] = pjj[i]; 236 #else 237 for (i=0; i<pdiag->nz; i++) hjj[i] = cstart + pjj[i]; 238 #endif 239 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar)); 240 if (sameint) { 241 ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 242 } else { 243 for (i=0; i<pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)(poffd->i[i]); 244 } 245 246 /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this 247 If we hacked a hypre a bit more we might be able to avoid this step */ 248 #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0) 249 PetscStackCallStandard(hypre_CSRMatrixBigInitialize,(hoffd)); 250 jj = (PetscInt*) hoffd->big_j; 251 #else 252 jj = (PetscInt*) hoffd->j; 253 #endif 254 pjj = poffd->j; 255 for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]]; 256 257 ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar)); 258 259 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 260 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 261 PetscFunctionReturn(0); 262 } 263 264 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B) 265 { 266 Mat_HYPRE* mhA = (Mat_HYPRE*)(A->data); 267 Mat lA; 268 ISLocalToGlobalMapping rl2g,cl2g; 269 IS is; 270 hypre_ParCSRMatrix *hA; 271 hypre_CSRMatrix *hdiag,*hoffd; 272 MPI_Comm comm; 273 HYPRE_Complex *hdd,*hod,*aa; 274 PetscScalar *data; 275 HYPRE_BigInt *col_map_offd; 276 HYPRE_Int *hdi,*hdj,*hoi,*hoj; 277 PetscInt *ii,*jj,*iptr,*jptr; 278 PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N; 279 HYPRE_Int type; 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 comm = PetscObjectComm((PetscObject)A); 284 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type)); 285 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 286 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA)); 287 M = hypre_ParCSRMatrixGlobalNumRows(hA); 288 N = hypre_ParCSRMatrixGlobalNumCols(hA); 289 str = hypre_ParCSRMatrixFirstRowIndex(hA); 290 stc = hypre_ParCSRMatrixFirstColDiag(hA); 291 hdiag = hypre_ParCSRMatrixDiag(hA); 292 hoffd = hypre_ParCSRMatrixOffd(hA); 293 dr = hypre_CSRMatrixNumRows(hdiag); 294 dc = hypre_CSRMatrixNumCols(hdiag); 295 nnz = hypre_CSRMatrixNumNonzeros(hdiag); 296 hdi = hypre_CSRMatrixI(hdiag); 297 hdj = hypre_CSRMatrixJ(hdiag); 298 hdd = hypre_CSRMatrixData(hdiag); 299 oc = hypre_CSRMatrixNumCols(hoffd); 300 nnz += hypre_CSRMatrixNumNonzeros(hoffd); 301 hoi = hypre_CSRMatrixI(hoffd); 302 hoj = hypre_CSRMatrixJ(hoffd); 303 hod = hypre_CSRMatrixData(hoffd); 304 if (reuse != MAT_REUSE_MATRIX) { 305 PetscInt *aux; 306 307 /* generate l2g maps for rows and cols */ 308 ierr = ISCreateStride(comm,dr,str,1,&is);CHKERRQ(ierr); 309 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 310 ierr = ISDestroy(&is);CHKERRQ(ierr); 311 col_map_offd = hypre_ParCSRMatrixColMapOffd(hA); 312 ierr = PetscMalloc1(dc+oc,&aux);CHKERRQ(ierr); 313 for (i=0; i<dc; i++) aux[i] = i+stc; 314 for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i]; 315 ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 316 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 317 ierr = ISDestroy(&is);CHKERRQ(ierr); 318 /* create MATIS object */ 319 ierr = MatCreate(comm,B);CHKERRQ(ierr); 320 ierr = MatSetSizes(*B,dr,dc,M,N);CHKERRQ(ierr); 321 ierr = MatSetType(*B,MATIS);CHKERRQ(ierr); 322 ierr = MatSetLocalToGlobalMapping(*B,rl2g,cl2g);CHKERRQ(ierr); 323 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 324 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 325 326 /* allocate CSR for local matrix */ 327 ierr = PetscMalloc1(dr+1,&iptr);CHKERRQ(ierr); 328 ierr = PetscMalloc1(nnz,&jptr);CHKERRQ(ierr); 329 ierr = PetscMalloc1(nnz,&data);CHKERRQ(ierr); 330 } else { 331 PetscInt nr; 332 PetscBool done; 333 ierr = MatISGetLocalMat(*B,&lA);CHKERRQ(ierr); 334 ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);CHKERRQ(ierr); 335 if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr); 336 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); 337 ierr = MatSeqAIJGetArray(lA,&data);CHKERRQ(ierr); 338 } 339 /* merge local matrices */ 340 ii = iptr; 341 jj = jptr; 342 aa = (HYPRE_Complex*)data; /* this cast fixes the clang error when doing the assignments below: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */ 343 *ii = *(hdi++) + *(hoi++); 344 for (jd=0,jo=0,cum=0; *ii<nnz; cum++) { 345 PetscScalar *aold = (PetscScalar*)aa; 346 PetscInt *jold = jj,nc = jd+jo; 347 for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; } 348 for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; } 349 *(++ii) = *(hdi++) + *(hoi++); 350 ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);CHKERRQ(ierr); 351 } 352 for (; cum<dr; cum++) *(++ii) = nnz; 353 if (reuse != MAT_REUSE_MATRIX) { 354 Mat_SeqAIJ* a; 355 356 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr); 357 ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr); 358 /* hack SeqAIJ */ 359 a = (Mat_SeqAIJ*)(lA->data); 360 a->free_a = PETSC_TRUE; 361 a->free_ij = PETSC_TRUE; 362 ierr = MatDestroy(&lA);CHKERRQ(ierr); 363 } 364 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 365 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 366 if (reuse == MAT_INPLACE_MATRIX) { 367 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 368 } 369 PetscFunctionReturn(0); 370 } 371 372 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B) 373 { 374 Mat M = NULL; 375 Mat_HYPRE *hB; 376 MPI_Comm comm = PetscObjectComm((PetscObject)A); 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 if (reuse == MAT_REUSE_MATRIX) { 381 /* always destroy the old matrix and create a new memory; 382 hope this does not churn the memory too much. The problem 383 is I do not know if it is possible to put the matrix back to 384 its initial state so that we can directly copy the values 385 the second time through. */ 386 hB = (Mat_HYPRE*)((*B)->data); 387 PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij)); 388 } else { 389 ierr = MatCreate(comm,&M);CHKERRQ(ierr); 390 ierr = MatSetType(M,MATHYPRE);CHKERRQ(ierr); 391 ierr = MatSetSizes(M,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 392 hB = (Mat_HYPRE*)(M->data); 393 if (reuse == MAT_INITIAL_MATRIX) *B = M; 394 } 395 ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr); 396 ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr); 397 if (reuse == MAT_INPLACE_MATRIX) { 398 ierr = MatHeaderReplace(A,&M);CHKERRQ(ierr); 399 } 400 (*B)->preallocated = PETSC_TRUE; 401 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 402 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 403 PetscFunctionReturn(0); 404 } 405 406 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B) 407 { 408 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 409 hypre_ParCSRMatrix *parcsr; 410 hypre_CSRMatrix *hdiag,*hoffd; 411 MPI_Comm comm; 412 PetscScalar *da,*oa,*aptr; 413 PetscInt *dii,*djj,*oii,*ojj,*iptr; 414 PetscInt i,dnnz,onnz,m,n; 415 HYPRE_Int type; 416 PetscMPIInt size; 417 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 comm = PetscObjectComm((PetscObject)A); 422 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 423 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 424 if (reuse == MAT_REUSE_MATRIX) { 425 PetscBool ismpiaij,isseqaij; 426 ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 427 ierr = PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 428 if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported"); 429 } 430 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 431 432 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr)); 433 hdiag = hypre_ParCSRMatrixDiag(parcsr); 434 hoffd = hypre_ParCSRMatrixOffd(parcsr); 435 m = hypre_CSRMatrixNumRows(hdiag); 436 n = hypre_CSRMatrixNumCols(hdiag); 437 dnnz = hypre_CSRMatrixNumNonzeros(hdiag); 438 onnz = hypre_CSRMatrixNumNonzeros(hoffd); 439 if (reuse == MAT_INITIAL_MATRIX) { 440 ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr); 441 ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr); 442 ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr); 443 } else if (reuse == MAT_REUSE_MATRIX) { 444 PetscInt nr; 445 PetscBool done; 446 if (size > 1) { 447 Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data); 448 449 ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr); 450 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); 451 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); 452 ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr); 453 } else { 454 ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr); 455 if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m); 456 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); 457 ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr); 458 } 459 } else { /* MAT_INPLACE_MATRIX */ 460 if (!sameint) { 461 ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr); 462 ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr); 463 } else { 464 dii = (PetscInt*)hypre_CSRMatrixI(hdiag); 465 djj = (PetscInt*)hypre_CSRMatrixJ(hdiag); 466 } 467 da = (PetscScalar*)hypre_CSRMatrixData(hdiag); 468 } 469 470 if (!sameint) { 471 for (i=0;i<m+1;i++) dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]); 472 for (i=0;i<dnnz;i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]); 473 } else { 474 ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 475 ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr); 476 } 477 ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr); 478 iptr = djj; 479 aptr = da; 480 for (i=0; i<m; i++) { 481 PetscInt nc = dii[i+1]-dii[i]; 482 ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr); 483 iptr += nc; 484 aptr += nc; 485 } 486 if (size > 1) { 487 HYPRE_BigInt *coffd; 488 HYPRE_Int *offdj; 489 490 if (reuse == MAT_INITIAL_MATRIX) { 491 ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr); 492 ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr); 493 ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr); 494 } else if (reuse == MAT_REUSE_MATRIX) { 495 Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data); 496 PetscInt nr,hr = hypre_CSRMatrixNumRows(hoffd); 497 PetscBool done; 498 499 ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr); 500 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); 501 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); 502 ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr); 503 } else { /* MAT_INPLACE_MATRIX */ 504 if (!sameint) { 505 ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr); 506 ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr); 507 } else { 508 oii = (PetscInt*)hypre_CSRMatrixI(hoffd); 509 ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd); 510 } 511 oa = (PetscScalar*)hypre_CSRMatrixData(hoffd); 512 } 513 if (!sameint) { 514 for (i=0;i<m+1;i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]); 515 } else { 516 ierr = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 517 } 518 offdj = hypre_CSRMatrixJ(hoffd); 519 coffd = hypre_ParCSRMatrixColMapOffd(parcsr); 520 for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]]; 521 ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr); 522 iptr = ojj; 523 aptr = oa; 524 for (i=0; i<m; i++) { 525 PetscInt nc = oii[i+1]-oii[i]; 526 ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr); 527 iptr += nc; 528 aptr += nc; 529 } 530 if (reuse == MAT_INITIAL_MATRIX) { 531 Mat_MPIAIJ *b; 532 Mat_SeqAIJ *d,*o; 533 534 ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);CHKERRQ(ierr); 535 /* hack MPIAIJ */ 536 b = (Mat_MPIAIJ*)((*B)->data); 537 d = (Mat_SeqAIJ*)b->A->data; 538 o = (Mat_SeqAIJ*)b->B->data; 539 d->free_a = PETSC_TRUE; 540 d->free_ij = PETSC_TRUE; 541 o->free_a = PETSC_TRUE; 542 o->free_ij = PETSC_TRUE; 543 } else if (reuse == MAT_INPLACE_MATRIX) { 544 Mat T; 545 546 ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);CHKERRQ(ierr); 547 if (sameint) { /* ownership of CSR pointers is transferred to PETSc */ 548 hypre_CSRMatrixI(hdiag) = NULL; 549 hypre_CSRMatrixJ(hdiag) = NULL; 550 hypre_CSRMatrixI(hoffd) = NULL; 551 hypre_CSRMatrixJ(hoffd) = NULL; 552 } else { /* Hack MPIAIJ -> free ij but not a */ 553 Mat_MPIAIJ *b = (Mat_MPIAIJ*)(T->data); 554 Mat_SeqAIJ *d = (Mat_SeqAIJ*)(b->A->data); 555 Mat_SeqAIJ *o = (Mat_SeqAIJ*)(b->B->data); 556 557 d->free_ij = PETSC_TRUE; 558 o->free_ij = PETSC_TRUE; 559 } 560 hypre_CSRMatrixData(hdiag) = NULL; 561 hypre_CSRMatrixData(hoffd) = NULL; 562 ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr); 563 } 564 } else { 565 oii = NULL; 566 ojj = NULL; 567 oa = NULL; 568 if (reuse == MAT_INITIAL_MATRIX) { 569 Mat_SeqAIJ* b; 570 571 ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr); 572 /* hack SeqAIJ */ 573 b = (Mat_SeqAIJ*)((*B)->data); 574 b->free_a = PETSC_TRUE; 575 b->free_ij = PETSC_TRUE; 576 } else if (reuse == MAT_INPLACE_MATRIX) { 577 Mat T; 578 579 ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);CHKERRQ(ierr); 580 if (sameint) { /* ownership of CSR pointers is transferred to PETSc */ 581 hypre_CSRMatrixI(hdiag) = NULL; 582 hypre_CSRMatrixJ(hdiag) = NULL; 583 } else { /* free ij but not a */ 584 Mat_SeqAIJ* b = (Mat_SeqAIJ*)(T->data); 585 586 b->free_ij = PETSC_TRUE; 587 } 588 hypre_CSRMatrixData(hdiag) = NULL; 589 ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr); 590 } 591 } 592 593 /* we have to use hypre_Tfree to free the HYPRE arrays 594 that PETSc now onws */ 595 if (reuse == MAT_INPLACE_MATRIX) { 596 PetscInt nh; 597 void *ptrs[6] = {da,oa,dii,djj,oii,ojj}; 598 const char *names[6] = {"_hypre_csr_da", 599 "_hypre_csr_oa", 600 "_hypre_csr_dii", 601 "_hypre_csr_djj", 602 "_hypre_csr_oii", 603 "_hypre_csr_ojj"}; 604 nh = sameint ? 6 : 2; 605 for (i=0; i<nh; i++) { 606 PetscContainer c; 607 608 ierr = PetscContainerCreate(comm,&c);CHKERRQ(ierr); 609 ierr = PetscContainerSetPointer(c,ptrs[i]);CHKERRQ(ierr); 610 ierr = PetscContainerSetUserDestroy(c,hypre_array_destroy);CHKERRQ(ierr); 611 ierr = PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);CHKERRQ(ierr); 612 ierr = PetscContainerDestroy(&c);CHKERRQ(ierr); 613 } 614 } 615 PetscFunctionReturn(0); 616 } 617 618 static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA) 619 { 620 hypre_ParCSRMatrix *tA; 621 hypre_CSRMatrix *hdiag,*hoffd; 622 Mat_SeqAIJ *diag,*offd; 623 PetscInt *garray,i,noffd,dnnz,onnz,*row_starts,*col_starts; 624 MPI_Comm comm = PetscObjectComm((PetscObject)A); 625 PetscBool ismpiaij,isseqaij; 626 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 627 PetscErrorCode ierr; 628 629 PetscFunctionBegin; 630 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 631 ierr = PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 632 if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type); 633 if (ismpiaij) { 634 Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data); 635 636 diag = (Mat_SeqAIJ*)a->A->data; 637 offd = (Mat_SeqAIJ*)a->B->data; 638 garray = a->garray; 639 noffd = a->B->cmap->N; 640 dnnz = diag->nz; 641 onnz = offd->nz; 642 } else { 643 diag = (Mat_SeqAIJ*)A->data; 644 offd = NULL; 645 garray = NULL; 646 noffd = 0; 647 dnnz = diag->nz; 648 onnz = 0; 649 } 650 651 /* create a temporary ParCSR */ 652 if (HYPRE_AssumedPartitionCheck()) { 653 PetscMPIInt myid; 654 655 ierr = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr); 656 row_starts = A->rmap->range + myid; 657 col_starts = A->cmap->range + myid; 658 } else { 659 row_starts = A->rmap->range; 660 col_starts = A->cmap->range; 661 } 662 tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_BigInt*)row_starts,(HYPRE_BigInt*)col_starts,noffd,dnnz,onnz); 663 hypre_ParCSRMatrixSetRowStartsOwner(tA,0); 664 hypre_ParCSRMatrixSetColStartsOwner(tA,0); 665 666 /* set diagonal part */ 667 hdiag = hypre_ParCSRMatrixDiag(tA); 668 if (sameint) { /* reuse CSR pointers */ 669 hypre_CSRMatrixI(hdiag) = (HYPRE_Int*)diag->i; 670 hypre_CSRMatrixJ(hdiag) = (HYPRE_Int*)diag->j; 671 } else { /* malloc CSR pointers */ 672 HYPRE_Int *hi,*hj; 673 674 ierr = PetscMalloc2(A->rmap->n+1,&hi,dnnz,&hj);CHKERRQ(ierr); 675 for (i = 0; i < A->rmap->n+1; i++) hi[i] = (HYPRE_Int)(diag->i[i]); 676 for (i = 0; i < dnnz; i++) hj[i] = (HYPRE_Int)(diag->j[i]); 677 hypre_CSRMatrixI(hdiag) = hi; 678 hypre_CSRMatrixJ(hdiag) = hj; 679 } 680 hypre_CSRMatrixData(hdiag) = (HYPRE_Complex*)diag->a; 681 hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz; 682 hypre_CSRMatrixSetRownnz(hdiag); 683 hypre_CSRMatrixSetDataOwner(hdiag,0); 684 685 /* set offdiagonal part */ 686 hoffd = hypre_ParCSRMatrixOffd(tA); 687 if (offd) { 688 if (sameint) { /* reuse CSR pointers */ 689 hypre_CSRMatrixI(hoffd) = (HYPRE_Int*)offd->i; 690 hypre_CSRMatrixJ(hoffd) = (HYPRE_Int*)offd->j; 691 } else { /* malloc CSR pointers */ 692 HYPRE_Int *hi,*hj; 693 694 ierr = PetscMalloc2(A->rmap->n+1,&hi,onnz,&hj);CHKERRQ(ierr); 695 for (i = 0; i < A->rmap->n+1; i++) hi[i] = (HYPRE_Int)(offd->i[i]); 696 for (i = 0; i < onnz; i++) hj[i] = (HYPRE_Int)(offd->j[i]); 697 hypre_CSRMatrixI(hoffd) = hi; 698 hypre_CSRMatrixJ(hoffd) = hj; 699 } 700 hypre_CSRMatrixData(hoffd) = (HYPRE_Complex*)offd->a; 701 hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz; 702 hypre_CSRMatrixSetRownnz(hoffd); 703 hypre_CSRMatrixSetDataOwner(hoffd,0); 704 hypre_ParCSRMatrixSetNumNonzeros(tA); 705 hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt*)garray; 706 } 707 *hA = tA; 708 PetscFunctionReturn(0); 709 } 710 711 static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA) 712 { 713 hypre_CSRMatrix *hdiag,*hoffd; 714 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 715 PetscErrorCode ierr; 716 717 PetscFunctionBegin; 718 hdiag = hypre_ParCSRMatrixDiag(*hA); 719 hoffd = hypre_ParCSRMatrixOffd(*hA); 720 /* free temporary memory allocated by PETSc */ 721 if (!sameint) { 722 HYPRE_Int *hi,*hj; 723 724 hi = hypre_CSRMatrixI(hdiag); 725 hj = hypre_CSRMatrixJ(hdiag); 726 ierr = PetscFree2(hi,hj);CHKERRQ(ierr); 727 if (hoffd) { 728 hi = hypre_CSRMatrixI(hoffd); 729 hj = hypre_CSRMatrixJ(hoffd); 730 ierr = PetscFree2(hi,hj);CHKERRQ(ierr); 731 } 732 } 733 /* set pointers to NULL before destroying tA */ 734 hypre_CSRMatrixI(hdiag) = NULL; 735 hypre_CSRMatrixJ(hdiag) = NULL; 736 hypre_CSRMatrixData(hdiag) = NULL; 737 hypre_CSRMatrixI(hoffd) = NULL; 738 hypre_CSRMatrixJ(hoffd) = NULL; 739 hypre_CSRMatrixData(hoffd) = NULL; 740 hypre_ParCSRMatrixColMapOffd(*hA) = NULL; 741 hypre_ParCSRMatrixDestroy(*hA); 742 *hA = NULL; 743 PetscFunctionReturn(0); 744 } 745 746 /* calls RAP from BoomerAMG: 747 the resulting ParCSR will not own the column and row starts 748 It looks like we don't need to have the diagonal entries 749 ordered first in the rows of the diagonal part 750 for boomerAMGBuildCoarseOperator to work */ 751 static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP) 752 { 753 HYPRE_Int P_owns_col_starts,R_owns_row_starts; 754 755 PetscFunctionBegin; 756 P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP); 757 R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR); 758 PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP)); 759 PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP)); 760 /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */ 761 hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0); 762 hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0); 763 if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1); 764 if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1); 765 PetscFunctionReturn(0); 766 } 767 768 static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C) 769 { 770 Mat B; 771 hypre_ParCSRMatrix *hA,*hP,*hPtAP; 772 PetscErrorCode ierr; 773 774 PetscFunctionBegin; 775 ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr); 776 ierr = MatAIJGetParCSR_Private(P,&hP);CHKERRQ(ierr); 777 ierr = MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);CHKERRQ(ierr); 778 ierr = MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);CHKERRQ(ierr); 779 ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr); 780 ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr); 781 ierr = MatAIJRestoreParCSR_Private(P,&hP);CHKERRQ(ierr); 782 PetscFunctionReturn(0); 783 } 784 785 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat *C) 786 { 787 PetscErrorCode ierr; 788 789 PetscFunctionBegin; 790 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 791 ierr = MatSetType(*C,MATAIJ);CHKERRQ(ierr); 792 (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE; 793 PetscFunctionReturn(0); 794 } 795 796 static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C) 797 { 798 Mat B; 799 Mat_HYPRE *hP; 800 hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr; 801 HYPRE_Int type; 802 MPI_Comm comm = PetscObjectComm((PetscObject)A); 803 PetscBool ishypre; 804 PetscErrorCode ierr; 805 806 PetscFunctionBegin; 807 ierr = PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);CHKERRQ(ierr); 808 if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE); 809 hP = (Mat_HYPRE*)P->data; 810 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type)); 811 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 812 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr)); 813 814 ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr); 815 ierr = MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);CHKERRQ(ierr); 816 ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr); 817 818 /* create temporary matrix and merge to C */ 819 ierr = MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);CHKERRQ(ierr); 820 ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C) 825 { 826 PetscErrorCode ierr; 827 828 PetscFunctionBegin; 829 if (scall == MAT_INITIAL_MATRIX) { 830 const char *deft = MATAIJ; 831 char type[256]; 832 PetscBool flg; 833 834 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 835 ierr = PetscOptionsFList("-matptap_hypre_outtype","Matrix type",NULL,MatList,deft,type,256,&flg);CHKERRQ(ierr); 836 ierr = PetscOptionsEnd();CHKERRQ(ierr); 837 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 838 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 839 if (flg) { 840 ierr = MatSetType(*C,type);CHKERRQ(ierr); 841 } else { 842 ierr = MatSetType(*C,deft);CHKERRQ(ierr); 843 } 844 (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE; 845 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 846 } 847 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 848 ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr); 849 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 850 PetscFunctionReturn(0); 851 } 852 853 static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C) 854 { 855 Mat B; 856 hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr; 857 Mat_HYPRE *hA,*hP; 858 PetscBool ishypre; 859 HYPRE_Int type; 860 PetscErrorCode ierr; 861 862 PetscFunctionBegin; 863 ierr = PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);CHKERRQ(ierr); 864 if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE); 865 ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr); 866 if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE); 867 hA = (Mat_HYPRE*)A->data; 868 hP = (Mat_HYPRE*)P->data; 869 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 870 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 871 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type)); 872 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 873 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr)); 874 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr)); 875 ierr = MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);CHKERRQ(ierr); 876 ierr = MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);CHKERRQ(ierr); 877 ierr = MatHeaderMerge(C,&B);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881 static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C) 882 { 883 PetscErrorCode ierr; 884 885 PetscFunctionBegin; 886 if (scall == MAT_INITIAL_MATRIX) { 887 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 888 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 889 ierr = MatSetType(*C,MATHYPRE);CHKERRQ(ierr); 890 (*C)->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE; 891 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 892 } 893 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 894 ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr); 895 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 896 PetscFunctionReturn(0); 897 } 898 899 /* calls hypre_ParMatmul 900 hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA 901 hypre_ParMatrixCreate does not duplicate the communicator 902 It looks like we don't need to have the diagonal entries 903 ordered first in the rows of the diagonal part 904 for boomerAMGBuildCoarseOperator to work */ 905 static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB) 906 { 907 PetscFunctionBegin; 908 PetscStackPush("hypre_ParMatmul"); 909 *hAB = hypre_ParMatmul(hA,hB); 910 PetscStackPop; 911 PetscFunctionReturn(0); 912 } 913 914 static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C) 915 { 916 Mat D; 917 hypre_ParCSRMatrix *hA,*hB,*hAB = NULL; 918 PetscErrorCode ierr; 919 920 PetscFunctionBegin; 921 ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr); 922 ierr = MatAIJGetParCSR_Private(B,&hB);CHKERRQ(ierr); 923 ierr = MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);CHKERRQ(ierr); 924 ierr = MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);CHKERRQ(ierr); 925 ierr = MatHeaderMerge(C,&D);CHKERRQ(ierr); 926 ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr); 927 ierr = MatAIJRestoreParCSR_Private(B,&hB);CHKERRQ(ierr); 928 PetscFunctionReturn(0); 929 } 930 931 PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat *C) 932 { 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 937 ierr = MatSetType(*C,MATAIJ);CHKERRQ(ierr); 938 (*C)->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE; 939 PetscFunctionReturn(0); 940 } 941 942 static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C) 943 { 944 Mat D; 945 hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL; 946 Mat_HYPRE *hA,*hB; 947 PetscBool ishypre; 948 HYPRE_Int type; 949 PetscErrorCode ierr; 950 951 PetscFunctionBegin; 952 ierr = PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);CHKERRQ(ierr); 953 if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE); 954 ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr); 955 if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE); 956 hA = (Mat_HYPRE*)A->data; 957 hB = (Mat_HYPRE*)B->data; 958 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 959 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 960 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type)); 961 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 962 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr)); 963 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr)); 964 ierr = MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);CHKERRQ(ierr); 965 ierr = MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);CHKERRQ(ierr); 966 /* need to use HeaderReplace because HeaderMerge messes up with the communicator */ 967 ierr = MatHeaderReplace(C,&D);CHKERRQ(ierr); 968 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE; 969 PetscFunctionReturn(0); 970 } 971 972 static PetscErrorCode MatMatMult_HYPRE_HYPRE(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 973 { 974 PetscErrorCode ierr; 975 976 PetscFunctionBegin; 977 if (scall == MAT_INITIAL_MATRIX) { 978 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 979 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 980 ierr = MatSetType(*C,MATHYPRE);CHKERRQ(ierr); 981 (*C)->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE; 982 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 983 } 984 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 985 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 986 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 987 PetscFunctionReturn(0); 988 } 989 990 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D) 991 { 992 Mat E; 993 hypre_ParCSRMatrix *hA,*hB,*hC,*hABC; 994 PetscErrorCode ierr; 995 996 PetscFunctionBegin; 997 ierr = MatAIJGetParCSR_Private(A,&hA);CHKERRQ(ierr); 998 ierr = MatAIJGetParCSR_Private(B,&hB);CHKERRQ(ierr); 999 ierr = MatAIJGetParCSR_Private(C,&hC);CHKERRQ(ierr); 1000 ierr = MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);CHKERRQ(ierr); 1001 ierr = MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);CHKERRQ(ierr); 1002 ierr = MatHeaderMerge(D,&E);CHKERRQ(ierr); 1003 ierr = MatAIJRestoreParCSR_Private(A,&hA);CHKERRQ(ierr); 1004 ierr = MatAIJRestoreParCSR_Private(B,&hB);CHKERRQ(ierr); 1005 ierr = MatAIJRestoreParCSR_Private(C,&hC);CHKERRQ(ierr); 1006 PetscFunctionReturn(0); 1007 } 1008 1009 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat *D) 1010 { 1011 PetscErrorCode ierr; 1012 1013 PetscFunctionBegin; 1014 ierr = MatCreate(PetscObjectComm((PetscObject)A),D);CHKERRQ(ierr); 1015 ierr = MatSetType(*D,MATAIJ);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y) 1020 { 1021 PetscErrorCode ierr; 1022 1023 PetscFunctionBegin; 1024 ierr = MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 1028 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y) 1029 { 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 ierr = MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);CHKERRQ(ierr); 1034 PetscFunctionReturn(0); 1035 } 1036 1037 static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1038 { 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 if (y != z) { 1043 ierr = VecCopy(y,z);CHKERRQ(ierr); 1044 } 1045 ierr = MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);CHKERRQ(ierr); 1046 PetscFunctionReturn(0); 1047 } 1048 1049 static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1050 { 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 if (y != z) { 1055 ierr = VecCopy(y,z);CHKERRQ(ierr); 1056 } 1057 ierr = MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);CHKERRQ(ierr); 1058 PetscFunctionReturn(0); 1059 } 1060 1061 /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */ 1062 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans) 1063 { 1064 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1065 hypre_ParCSRMatrix *parcsr; 1066 hypre_ParVector *hx,*hy; 1067 HYPRE_Complex *ax,*ay,*sax,*say; 1068 PetscErrorCode ierr; 1069 1070 PetscFunctionBegin; 1071 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr)); 1072 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx)); 1073 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy)); 1074 ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr); 1075 ierr = VecGetArray(y,(PetscScalar**)&ay);CHKERRQ(ierr); 1076 if (trans) { 1077 VecHYPRE_ParVectorReplacePointer(hA->x,ay,say); 1078 VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax); 1079 hypre_ParCSRMatrixMatvecT(a,parcsr,hy,b,hx); 1080 VecHYPRE_ParVectorReplacePointer(hA->x,say,ay); 1081 VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax); 1082 } else { 1083 VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax); 1084 VecHYPRE_ParVectorReplacePointer(hA->b,ay,say); 1085 hypre_ParCSRMatrixMatvec(a,parcsr,hx,b,hy); 1086 VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax); 1087 VecHYPRE_ParVectorReplacePointer(hA->b,say,ay); 1088 } 1089 ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr); 1090 ierr = VecRestoreArray(y,(PetscScalar**)&ay);CHKERRQ(ierr); 1091 PetscFunctionReturn(0); 1092 } 1093 1094 static PetscErrorCode MatDestroy_HYPRE(Mat A) 1095 { 1096 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1097 PetscErrorCode ierr; 1098 1099 PetscFunctionBegin; 1100 if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x)); 1101 if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b)); 1102 if (hA->ij) { 1103 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL; 1104 PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij)); 1105 } 1106 if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); } 1107 1108 ierr = MatStashDestroy_Private(&A->stash);CHKERRQ(ierr); 1109 1110 ierr = PetscFree(hA->array);CHKERRQ(ierr); 1111 1112 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr); 1113 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr); 1114 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);CHKERRQ(ierr); 1115 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);CHKERRQ(ierr); 1116 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);CHKERRQ(ierr); 1117 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);CHKERRQ(ierr); 1118 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_hypre_C",NULL);CHKERRQ(ierr); 1119 ierr = PetscFree(A->data);CHKERRQ(ierr); 1120 PetscFunctionReturn(0); 1121 } 1122 1123 static PetscErrorCode MatSetUp_HYPRE(Mat A) 1124 { 1125 PetscErrorCode ierr; 1126 1127 PetscFunctionBegin; 1128 ierr = MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode) 1133 { 1134 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1135 Vec x,b; 1136 PetscMPIInt n; 1137 PetscInt i,j,rstart,ncols,flg; 1138 PetscInt *row,*col; 1139 PetscScalar *val; 1140 PetscErrorCode ierr; 1141 1142 PetscFunctionBegin; 1143 if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE"); 1144 1145 if (!A->nooffprocentries) { 1146 while (1) { 1147 ierr = MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);CHKERRQ(ierr); 1148 if (!flg) break; 1149 1150 for (i=0; i<n; ) { 1151 /* Now identify the consecutive vals belonging to the same row */ 1152 for (j=i,rstart=row[j]; j<n; j++) { 1153 if (row[j] != rstart) break; 1154 } 1155 if (j < n) ncols = j-i; 1156 else ncols = n-i; 1157 /* Now assemble all these values with a single function call */ 1158 ierr = MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);CHKERRQ(ierr); 1159 1160 i = j; 1161 } 1162 } 1163 ierr = MatStashScatterEnd_Private(&A->stash);CHKERRQ(ierr); 1164 } 1165 1166 PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij)); 1167 /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitalize */ 1168 { 1169 hypre_AuxParCSRMatrix *aux_matrix; 1170 1171 /* call destroy just to make sure we do not leak anything */ 1172 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij); 1173 PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,(aux_matrix)); 1174 hypre_IJMatrixTranslator(hA->ij) = NULL; 1175 1176 /* Initialize with assembled flag -> it only recreates the aux_par_matrix */ 1177 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij)); 1178 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij); 1179 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */ 1180 PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,(aux_matrix)); 1181 } 1182 if (hA->x) PetscFunctionReturn(0); 1183 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1184 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1185 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr); 1186 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr); 1187 ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr); 1188 ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr); 1189 ierr = VecDestroy(&x);CHKERRQ(ierr); 1190 ierr = VecDestroy(&b);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 1194 static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array) 1195 { 1196 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1197 PetscErrorCode ierr; 1198 1199 PetscFunctionBegin; 1200 if (!hA->available) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Temporary space is in use"); 1201 1202 if (hA->size >= size) { 1203 *array = hA->array; 1204 } else { 1205 ierr = PetscFree(hA->array);CHKERRQ(ierr); 1206 hA->size = size; 1207 ierr = PetscMalloc(hA->size,&hA->array);CHKERRQ(ierr); 1208 *array = hA->array; 1209 } 1210 1211 hA->available = PETSC_FALSE; 1212 PetscFunctionReturn(0); 1213 } 1214 1215 static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array) 1216 { 1217 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1218 1219 PetscFunctionBegin; 1220 *array = NULL; 1221 hA->available = PETSC_TRUE; 1222 PetscFunctionReturn(0); 1223 } 1224 1225 PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins) 1226 { 1227 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1228 PetscScalar *vals = (PetscScalar *)v; 1229 HYPRE_Complex *sscr; 1230 PetscInt *cscr[2]; 1231 PetscInt i,nzc; 1232 void *array = NULL; 1233 PetscErrorCode ierr; 1234 1235 PetscFunctionBegin; 1236 ierr = MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(HYPRE_Complex)*nc*nr,&array);CHKERRQ(ierr); 1237 cscr[0] = (PetscInt*)array; 1238 cscr[1] = ((PetscInt*)array)+nc; 1239 sscr = (HYPRE_Complex*)(((PetscInt*)array)+nc*2); 1240 for (i=0,nzc=0;i<nc;i++) { 1241 if (cols[i] >= 0) { 1242 cscr[0][nzc ] = cols[i]; 1243 cscr[1][nzc++] = i; 1244 } 1245 } 1246 if (!nzc) { 1247 ierr = MatRestoreArray_HYPRE(A,&array);CHKERRQ(ierr); 1248 PetscFunctionReturn(0); 1249 } 1250 1251 if (ins == ADD_VALUES) { 1252 for (i=0;i<nr;i++) { 1253 if (rows[i] >= 0 && nzc) { 1254 PetscInt j; 1255 HYPRE_Int hnc = (HYPRE_Int)nzc; 1256 1257 if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]); 1258 for (j=0;j<nzc;j++) { ierr = PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); } 1259 PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr)); 1260 } 1261 vals += nc; 1262 } 1263 } else { /* INSERT_VALUES */ 1264 PetscInt rst,ren; 1265 1266 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr); 1267 for (i=0;i<nr;i++) { 1268 if (rows[i] >= 0 && nzc) { 1269 PetscInt j; 1270 HYPRE_Int hnc = (HYPRE_Int)nzc; 1271 1272 if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]); 1273 for (j=0;j<nzc;j++) { ierr = PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); } 1274 /* nonlocal values */ 1275 if (rows[i] < rst || rows[i] >= ren) { ierr = MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],(PetscScalar*)sscr,PETSC_FALSE);CHKERRQ(ierr); } 1276 /* local values */ 1277 else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr)); 1278 } 1279 vals += nc; 1280 } 1281 } 1282 1283 ierr = MatRestoreArray_HYPRE(A,&array);CHKERRQ(ierr); 1284 PetscFunctionReturn(0); 1285 } 1286 1287 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1288 { 1289 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1290 HYPRE_Int *hdnnz,*honnz; 1291 PetscInt i,rs,re,cs,ce,bs; 1292 PetscMPIInt size; 1293 PetscErrorCode ierr; 1294 1295 PetscFunctionBegin; 1296 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1297 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1298 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1299 rs = A->rmap->rstart; 1300 re = A->rmap->rend; 1301 cs = A->cmap->rstart; 1302 ce = A->cmap->rend; 1303 if (!hA->ij) { 1304 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij)); 1305 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 1306 } else { 1307 HYPRE_BigInt hrs,hre,hcs,hce; 1308 PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce)); 1309 if (hre-hrs+1 != re -rs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local rows: IJMatrix [%D,%D), PETSc [%D,%D)",hrs,hre+1,rs,re); 1310 if (hce-hcs+1 != ce -cs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local cols: IJMatrix [%D,%D), PETSc [%D,%D)",hcs,hce+1,cs,ce); 1311 } 1312 if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs; 1313 if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs; 1314 1315 if (!dnnz) { 1316 ierr = PetscMalloc1(A->rmap->n,&hdnnz);CHKERRQ(ierr); 1317 for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz; 1318 } else { 1319 hdnnz = (HYPRE_Int*)dnnz; 1320 } 1321 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1322 if (size > 1) { 1323 hypre_AuxParCSRMatrix *aux_matrix; 1324 if (!onnz) { 1325 ierr = PetscMalloc1(A->rmap->n,&honnz);CHKERRQ(ierr); 1326 for (i=0;i<A->rmap->n;i++) honnz[i] = onz; 1327 } else { 1328 honnz = (HYPRE_Int*)onnz; 1329 } 1330 /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems 1331 they assume the user will input the entire row values, properly sorted 1332 In PETSc, we don't make such an assumption, and we instead set this flag to 1 1333 Also, to avoid possible memory leaks, we destroy and recreate the translator 1334 This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize 1335 the IJ matrix for us */ 1336 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij); 1337 hypre_AuxParCSRMatrixDestroy(aux_matrix); 1338 hypre_IJMatrixTranslator(hA->ij) = NULL; 1339 PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz)); 1340 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij); 1341 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; 1342 } else { 1343 honnz = NULL; 1344 PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz)); 1345 } 1346 1347 /* reset assembled flag and call the initialize method */ 1348 hypre_IJMatrixAssembleFlag(hA->ij) = 0; 1349 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij)); 1350 if (!dnnz) { 1351 ierr = PetscFree(hdnnz);CHKERRQ(ierr); 1352 } 1353 if (!onnz && honnz) { 1354 ierr = PetscFree(honnz);CHKERRQ(ierr); 1355 } 1356 1357 /* Match AIJ logic */ 1358 A->preallocated = PETSC_TRUE; 1359 A->assembled = PETSC_FALSE; 1360 PetscFunctionReturn(0); 1361 } 1362 1363 /*@C 1364 MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format 1365 1366 Collective on Mat 1367 1368 Input Parameters: 1369 + A - the matrix 1370 . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix 1371 (same value is used for all local rows) 1372 . dnnz - array containing the number of nonzeros in the various rows of the 1373 DIAGONAL portion of the local submatrix (possibly different for each row) 1374 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1375 The size of this array is equal to the number of local rows, i.e 'm'. 1376 For matrices that will be factored, you must leave room for (and set) 1377 the diagonal entry even if it is zero. 1378 . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1379 submatrix (same value is used for all local rows). 1380 - onnz - array containing the number of nonzeros in the various rows of the 1381 OFF-DIAGONAL portion of the local submatrix (possibly different for 1382 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1383 structure. The size of this array is equal to the number 1384 of local rows, i.e 'm'. 1385 1386 Notes: 1387 If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored. 1388 1389 Level: intermediate 1390 1391 .keywords: matrix, aij, compressed row, sparse, parallel 1392 1393 .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE 1394 @*/ 1395 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1396 { 1397 PetscErrorCode ierr; 1398 1399 PetscFunctionBegin; 1400 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1401 PetscValidType(A,1); 1402 ierr = PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));CHKERRQ(ierr); 1403 PetscFunctionReturn(0); 1404 } 1405 1406 /* 1407 MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix 1408 1409 Collective 1410 1411 Input Parameters: 1412 + parcsr - the pointer to the hypre_ParCSRMatrix 1413 . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported. 1414 - copymode - PETSc copying options 1415 1416 Output Parameter: 1417 . A - the matrix 1418 1419 Level: intermediate 1420 1421 .seealso: MatHYPRE, PetscCopyMode 1422 */ 1423 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A) 1424 { 1425 Mat T; 1426 Mat_HYPRE *hA; 1427 MPI_Comm comm; 1428 PetscInt rstart,rend,cstart,cend,M,N; 1429 PetscBool isseqaij,ismpiaij,isaij,ishyp,isis; 1430 PetscErrorCode ierr; 1431 1432 PetscFunctionBegin; 1433 comm = hypre_ParCSRMatrixComm(parcsr); 1434 ierr = PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1435 ierr = PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 1436 ierr = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr); 1437 ierr = PetscStrcmp(mtype,MATHYPRE,&ishyp);CHKERRQ(ierr); 1438 ierr = PetscStrcmp(mtype,MATIS,&isis);CHKERRQ(ierr); 1439 isaij = (PetscBool)(isseqaij || ismpiaij || isaij); 1440 if (!isaij && !ishyp && !isis) SETERRQ6(comm,PETSC_ERR_SUP,"Unsupported MatType %s! Supported types are %s, %s, %s, %s, and %s",mtype,MATAIJ,MATSEQAIJ,MATMPIAIJ,MATIS,MATHYPRE); 1441 /* access ParCSRMatrix */ 1442 rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr); 1443 rend = hypre_ParCSRMatrixLastRowIndex(parcsr); 1444 cstart = hypre_ParCSRMatrixFirstColDiag(parcsr); 1445 cend = hypre_ParCSRMatrixLastColDiag(parcsr); 1446 M = hypre_ParCSRMatrixGlobalNumRows(parcsr); 1447 N = hypre_ParCSRMatrixGlobalNumCols(parcsr); 1448 1449 /* fix for empty local rows/columns */ 1450 if (rend < rstart) rend = rstart; 1451 if (cend < cstart) cend = cstart; 1452 1453 /* PETSc convention */ 1454 rend++; 1455 cend++; 1456 rend = PetscMin(rend,M); 1457 cend = PetscMin(cend,N); 1458 1459 /* create PETSc matrix with MatHYPRE */ 1460 ierr = MatCreate(comm,&T);CHKERRQ(ierr); 1461 ierr = MatSetSizes(T,rend-rstart,cend-cstart,M,N);CHKERRQ(ierr); 1462 ierr = MatSetType(T,MATHYPRE);CHKERRQ(ierr); 1463 hA = (Mat_HYPRE*)(T->data); 1464 1465 /* create HYPRE_IJMatrix */ 1466 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij)); 1467 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 1468 1469 /* create new ParCSR object if needed */ 1470 if (ishyp && copymode == PETSC_COPY_VALUES) { 1471 hypre_ParCSRMatrix *new_parcsr; 1472 hypre_CSRMatrix *hdiag,*hoffd,*ndiag,*noffd; 1473 1474 new_parcsr = hypre_ParCSRMatrixCompleteClone(parcsr); 1475 hdiag = hypre_ParCSRMatrixDiag(parcsr); 1476 hoffd = hypre_ParCSRMatrixOffd(parcsr); 1477 ndiag = hypre_ParCSRMatrixDiag(new_parcsr); 1478 noffd = hypre_ParCSRMatrixOffd(new_parcsr); 1479 ierr = PetscMemcpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag)*sizeof(HYPRE_Complex));CHKERRQ(ierr); 1480 ierr = PetscMemcpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd)*sizeof(HYPRE_Complex));CHKERRQ(ierr); 1481 parcsr = new_parcsr; 1482 copymode = PETSC_OWN_POINTER; 1483 } 1484 1485 /* set ParCSR object */ 1486 hypre_IJMatrixObject(hA->ij) = parcsr; 1487 T->preallocated = PETSC_TRUE; 1488 1489 /* set assembled flag */ 1490 hypre_IJMatrixAssembleFlag(hA->ij) = 1; 1491 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij)); 1492 if (ishyp) { 1493 PetscMPIInt myid = 0; 1494 1495 /* make sure we always have row_starts and col_starts available */ 1496 if (HYPRE_AssumedPartitionCheck()) { 1497 ierr = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr); 1498 } 1499 if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) { 1500 PetscLayout map; 1501 1502 ierr = MatGetLayouts(T,NULL,&map);CHKERRQ(ierr); 1503 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 1504 hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid); 1505 } 1506 if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) { 1507 PetscLayout map; 1508 1509 ierr = MatGetLayouts(T,&map,NULL);CHKERRQ(ierr); 1510 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 1511 hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid); 1512 } 1513 /* prevent from freeing the pointer */ 1514 if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE; 1515 *A = T; 1516 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1517 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1518 } else if (isaij) { 1519 if (copymode != PETSC_OWN_POINTER) { 1520 /* prevent from freeing the pointer */ 1521 hA->inner_free = PETSC_FALSE; 1522 ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr); 1523 ierr = MatDestroy(&T);CHKERRQ(ierr); 1524 } else { /* AIJ return type with PETSC_OWN_POINTER */ 1525 ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr); 1526 *A = T; 1527 } 1528 } else if (isis) { 1529 ierr = MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr); 1530 if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE; 1531 ierr = MatDestroy(&T);CHKERRQ(ierr); 1532 } 1533 PetscFunctionReturn(0); 1534 } 1535 1536 static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr) 1537 { 1538 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1539 HYPRE_Int type; 1540 1541 PetscFunctionBegin; 1542 if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present"); 1543 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 1544 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR"); 1545 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr)); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 /* 1550 MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix 1551 1552 Not collective 1553 1554 Input Parameters: 1555 + A - the MATHYPRE object 1556 1557 Output Parameter: 1558 . parcsr - the pointer to the hypre_ParCSRMatrix 1559 1560 Level: intermediate 1561 1562 .seealso: MatHYPRE, PetscCopyMode 1563 */ 1564 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr) 1565 { 1566 PetscErrorCode ierr; 1567 1568 PetscFunctionBegin; 1569 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1570 PetscValidType(A,1); 1571 ierr = PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));CHKERRQ(ierr); 1572 PetscFunctionReturn(0); 1573 } 1574 1575 static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd) 1576 { 1577 hypre_ParCSRMatrix *parcsr; 1578 hypre_CSRMatrix *ha; 1579 PetscInt rst; 1580 PetscErrorCode ierr; 1581 1582 PetscFunctionBegin; 1583 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks"); 1584 ierr = MatGetOwnershipRange(A,&rst,NULL);CHKERRQ(ierr); 1585 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1586 if (missing) *missing = PETSC_FALSE; 1587 if (dd) *dd = -1; 1588 ha = hypre_ParCSRMatrixDiag(parcsr); 1589 if (ha) { 1590 PetscInt size,i; 1591 HYPRE_Int *ii,*jj; 1592 1593 size = hypre_CSRMatrixNumRows(ha); 1594 ii = hypre_CSRMatrixI(ha); 1595 jj = hypre_CSRMatrixJ(ha); 1596 for (i = 0; i < size; i++) { 1597 PetscInt j; 1598 PetscBool found = PETSC_FALSE; 1599 1600 for (j = ii[i]; j < ii[i+1] && !found; j++) 1601 found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE; 1602 1603 if (!found) { 1604 PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i); 1605 if (missing) *missing = PETSC_TRUE; 1606 if (dd) *dd = i+rst; 1607 PetscFunctionReturn(0); 1608 } 1609 } 1610 if (!size) { 1611 PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n"); 1612 if (missing) *missing = PETSC_TRUE; 1613 if (dd) *dd = rst; 1614 } 1615 } else { 1616 PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n"); 1617 if (missing) *missing = PETSC_TRUE; 1618 if (dd) *dd = rst; 1619 } 1620 PetscFunctionReturn(0); 1621 } 1622 1623 static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s) 1624 { 1625 hypre_ParCSRMatrix *parcsr; 1626 hypre_CSRMatrix *ha; 1627 PetscErrorCode ierr; 1628 HYPRE_Complex hs; 1629 1630 PetscFunctionBegin; 1631 ierr = PetscHYPREScalarCast(s,&hs);CHKERRQ(ierr); 1632 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1633 /* diagonal part */ 1634 ha = hypre_ParCSRMatrixDiag(parcsr); 1635 if (ha) { 1636 PetscInt size,i; 1637 HYPRE_Int *ii; 1638 HYPRE_Complex *a; 1639 1640 size = hypre_CSRMatrixNumRows(ha); 1641 a = hypre_CSRMatrixData(ha); 1642 ii = hypre_CSRMatrixI(ha); 1643 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1644 } 1645 /* offdiagonal part */ 1646 ha = hypre_ParCSRMatrixOffd(parcsr); 1647 if (ha) { 1648 PetscInt size,i; 1649 HYPRE_Int *ii; 1650 HYPRE_Complex *a; 1651 1652 size = hypre_CSRMatrixNumRows(ha); 1653 a = hypre_CSRMatrixData(ha); 1654 ii = hypre_CSRMatrixI(ha); 1655 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1656 } 1657 PetscFunctionReturn(0); 1658 } 1659 1660 static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1661 { 1662 hypre_ParCSRMatrix *parcsr; 1663 HYPRE_Int *lrows; 1664 PetscInt rst,ren,i; 1665 PetscErrorCode ierr; 1666 1667 PetscFunctionBegin; 1668 if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented"); 1669 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1670 ierr = PetscMalloc1(numRows,&lrows);CHKERRQ(ierr); 1671 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr); 1672 for (i=0;i<numRows;i++) { 1673 if (rows[i] < rst || rows[i] >= ren) 1674 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported"); 1675 lrows[i] = rows[i] - rst; 1676 } 1677 PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows)); 1678 ierr = PetscFree(lrows);CHKERRQ(ierr); 1679 PetscFunctionReturn(0); 1680 } 1681 1682 static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha) 1683 { 1684 PetscErrorCode ierr; 1685 1686 PetscFunctionBegin; 1687 if (ha) { 1688 HYPRE_Int *ii, size; 1689 HYPRE_Complex *a; 1690 1691 size = hypre_CSRMatrixNumRows(ha); 1692 a = hypre_CSRMatrixData(ha); 1693 ii = hypre_CSRMatrixI(ha); 1694 1695 if (a) { ierr = PetscMemzero(a,(ii[size])*sizeof(HYPRE_Complex));CHKERRQ(ierr); } 1696 } 1697 PetscFunctionReturn(0); 1698 } 1699 1700 PetscErrorCode MatZeroEntries_HYPRE(Mat A) 1701 { 1702 hypre_ParCSRMatrix *parcsr; 1703 PetscErrorCode ierr; 1704 1705 PetscFunctionBegin; 1706 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1707 /* diagonal part */ 1708 ierr = MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));CHKERRQ(ierr); 1709 /* off-diagonal part */ 1710 ierr = MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));CHKERRQ(ierr); 1711 PetscFunctionReturn(0); 1712 } 1713 1714 static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag) 1715 { 1716 PetscInt ii; 1717 HYPRE_Int *i, *j; 1718 HYPRE_Complex *a; 1719 1720 PetscFunctionBegin; 1721 if (!hA) PetscFunctionReturn(0); 1722 1723 i = hypre_CSRMatrixI(hA); 1724 j = hypre_CSRMatrixJ(hA); 1725 a = hypre_CSRMatrixData(hA); 1726 1727 for (ii = 0; ii < N; ii++) { 1728 HYPRE_Int jj, ibeg, iend, irow; 1729 1730 irow = rows[ii]; 1731 ibeg = i[irow]; 1732 iend = i[irow+1]; 1733 for (jj = ibeg; jj < iend; jj++) 1734 if (j[jj] == irow) a[jj] = diag; 1735 else a[jj] = 0.0; 1736 } 1737 PetscFunctionReturn(0); 1738 } 1739 1740 static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1741 { 1742 hypre_ParCSRMatrix *parcsr; 1743 PetscInt *lrows,len; 1744 HYPRE_Complex hdiag; 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size"); 1749 ierr = PetscHYPREScalarCast(diag,&hdiag);CHKERRQ(ierr); 1750 /* retrieve the internal matrix */ 1751 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1752 /* get locally owned rows */ 1753 ierr = MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);CHKERRQ(ierr); 1754 /* zero diagonal part */ 1755 ierr = MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag);CHKERRQ(ierr); 1756 /* zero off-diagonal part */ 1757 ierr = MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);CHKERRQ(ierr); 1758 1759 ierr = PetscFree(lrows);CHKERRQ(ierr); 1760 PetscFunctionReturn(0); 1761 } 1762 1763 static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode) 1764 { 1765 PetscErrorCode ierr; 1766 1767 PetscFunctionBegin; 1768 if (mat->nooffprocentries) PetscFunctionReturn(0); 1769 1770 ierr = MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);CHKERRQ(ierr); 1771 PetscFunctionReturn(0); 1772 } 1773 1774 static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1775 { 1776 hypre_ParCSRMatrix *parcsr; 1777 HYPRE_Int hnz; 1778 PetscErrorCode ierr; 1779 1780 PetscFunctionBegin; 1781 /* retrieve the internal matrix */ 1782 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1783 /* call HYPRE API */ 1784 PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v)); 1785 if (nz) *nz = (PetscInt)hnz; 1786 PetscFunctionReturn(0); 1787 } 1788 1789 static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1790 { 1791 hypre_ParCSRMatrix *parcsr; 1792 HYPRE_Int hnz; 1793 PetscErrorCode ierr; 1794 1795 PetscFunctionBegin; 1796 /* retrieve the internal matrix */ 1797 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1798 /* call HYPRE API */ 1799 hnz = nz ? (HYPRE_Int)(*nz) : 0; 1800 PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v)); 1801 PetscFunctionReturn(0); 1802 } 1803 1804 static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 1805 { 1806 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1807 PetscInt i; 1808 1809 PetscFunctionBegin; 1810 if (!m || !n) PetscFunctionReturn(0); 1811 /* Ignore negative row indices 1812 * And negative column indices should be automatically ignored in hypre 1813 * */ 1814 for (i=0; i<m; i++) { 1815 if (idxm[i] >= 0) { 1816 HYPRE_Int hn = (HYPRE_Int)n; 1817 PetscStackCallStandard(HYPRE_IJMatrixGetValues,(hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n))); 1818 } 1819 } 1820 PetscFunctionReturn(0); 1821 } 1822 1823 static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg) 1824 { 1825 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1826 1827 PetscFunctionBegin; 1828 switch (op) { 1829 case MAT_NO_OFF_PROC_ENTRIES: 1830 if (flg) { 1831 PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,(hA->ij,0)); 1832 } 1833 break; 1834 default: 1835 break; 1836 } 1837 PetscFunctionReturn(0); 1838 } 1839 1840 static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view) 1841 { 1842 hypre_ParCSRMatrix *parcsr; 1843 PetscErrorCode ierr; 1844 Mat B; 1845 PetscViewerFormat format; 1846 PetscErrorCode (*mview)(Mat,PetscViewer) = NULL; 1847 1848 PetscFunctionBegin; 1849 ierr = PetscViewerGetFormat(view,&format);CHKERRQ(ierr); 1850 if (format != PETSC_VIEWER_NATIVE) { 1851 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1852 ierr = MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);CHKERRQ(ierr); 1853 ierr = MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);CHKERRQ(ierr); 1854 if (!mview) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation");CHKERRQ(ierr); 1855 ierr = (*mview)(B,view);CHKERRQ(ierr); 1856 ierr = MatDestroy(&B);CHKERRQ(ierr); 1857 } else { 1858 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1859 PetscMPIInt size; 1860 PetscBool isascii; 1861 const char *filename; 1862 1863 /* HYPRE uses only text files */ 1864 ierr = PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 1865 if (!isascii) SETERRQ1(PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name); 1866 ierr = PetscViewerFileGetName(view,&filename);CHKERRQ(ierr); 1867 PetscStackCallStandard(HYPRE_IJMatrixPrint,(hA->ij,filename)); 1868 ierr = MPI_Comm_size(hA->comm,&size);CHKERRQ(ierr); 1869 if (size > 1) { 1870 ierr = PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);CHKERRQ(ierr); 1871 } else { 1872 ierr = PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);CHKERRQ(ierr); 1873 } 1874 } 1875 PetscFunctionReturn(0); 1876 } 1877 1878 static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B) 1879 { 1880 hypre_ParCSRMatrix *parcsr; 1881 PetscErrorCode ierr; 1882 PetscCopyMode cpmode; 1883 1884 PetscFunctionBegin; 1885 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1886 if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) { 1887 parcsr = hypre_ParCSRMatrixCompleteClone(parcsr); 1888 cpmode = PETSC_OWN_POINTER; 1889 } else { 1890 cpmode = PETSC_COPY_VALUES; 1891 } 1892 ierr = MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 1896 static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str) 1897 { 1898 hypre_ParCSRMatrix *acsr,*bcsr; 1899 PetscErrorCode ierr; 1900 1901 PetscFunctionBegin; 1902 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1903 ierr = MatHYPREGetParCSR_HYPRE(A,&acsr);CHKERRQ(ierr); 1904 ierr = MatHYPREGetParCSR_HYPRE(B,&bcsr);CHKERRQ(ierr); 1905 PetscStackCallStandard(hypre_ParCSRMatrixCopy,(acsr,bcsr,1)); 1906 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1907 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1908 } else { 1909 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1910 } 1911 PetscFunctionReturn(0); 1912 } 1913 1914 static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d) 1915 { 1916 hypre_ParCSRMatrix *parcsr; 1917 hypre_CSRMatrix *dmat; 1918 HYPRE_Complex *a; 1919 HYPRE_Complex *data = NULL; 1920 HYPRE_Int *diag = NULL; 1921 PetscInt i; 1922 PetscBool cong; 1923 PetscErrorCode ierr; 1924 1925 PetscFunctionBegin; 1926 ierr = MatHasCongruentLayouts(A,&cong);CHKERRQ(ierr); 1927 if (!cong) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns"); 1928 #if defined(PETSC_USE_DEBUG) 1929 { 1930 PetscBool miss; 1931 ierr = MatMissingDiagonal(A,&miss,NULL);CHKERRQ(ierr); 1932 if (miss && A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing"); 1933 } 1934 #endif 1935 ierr = MatHYPREGetParCSR_HYPRE(A,&parcsr);CHKERRQ(ierr); 1936 dmat = hypre_ParCSRMatrixDiag(parcsr); 1937 if (dmat) { 1938 /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */ 1939 ierr = VecGetArray(d,(PetscScalar**)&a);CHKERRQ(ierr); 1940 diag = hypre_CSRMatrixI(dmat); 1941 data = hypre_CSRMatrixData(dmat); 1942 for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]]; 1943 ierr = VecRestoreArray(d,(PetscScalar**)&a);CHKERRQ(ierr); 1944 } 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #include <petscblaslapack.h> 1949 1950 static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str) 1951 { 1952 PetscErrorCode ierr; 1953 1954 PetscFunctionBegin; 1955 if (str == SAME_NONZERO_PATTERN) { 1956 hypre_ParCSRMatrix *x,*y; 1957 hypre_CSRMatrix *xloc,*yloc; 1958 PetscInt xnnz,ynnz; 1959 HYPRE_Complex *xarr,*yarr; 1960 PetscBLASInt one=1,bnz; 1961 1962 ierr = MatHYPREGetParCSR(Y,&y);CHKERRQ(ierr); 1963 ierr = MatHYPREGetParCSR(X,&x);CHKERRQ(ierr); 1964 1965 /* diagonal block */ 1966 xloc = hypre_ParCSRMatrixDiag(x); 1967 yloc = hypre_ParCSRMatrixDiag(y); 1968 xnnz = 0; 1969 ynnz = 0; 1970 xarr = NULL; 1971 yarr = NULL; 1972 if (xloc) { 1973 xarr = hypre_CSRMatrixData(xloc); 1974 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 1975 } 1976 if (yloc) { 1977 yarr = hypre_CSRMatrixData(yloc); 1978 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 1979 } 1980 if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %D != %D",xnnz,ynnz); 1981 ierr = PetscBLASIntCast(xnnz,&bnz);CHKERRQ(ierr); 1982 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one)); 1983 1984 /* off-diagonal block */ 1985 xloc = hypre_ParCSRMatrixOffd(x); 1986 yloc = hypre_ParCSRMatrixOffd(y); 1987 xnnz = 0; 1988 ynnz = 0; 1989 xarr = NULL; 1990 yarr = NULL; 1991 if (xloc) { 1992 xarr = hypre_CSRMatrixData(xloc); 1993 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 1994 } 1995 if (yloc) { 1996 yarr = hypre_CSRMatrixData(yloc); 1997 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 1998 } 1999 if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %D != %D",xnnz,ynnz); 2000 ierr = PetscBLASIntCast(xnnz,&bnz);CHKERRQ(ierr); 2001 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one)); 2002 } else if (str == SUBSET_NONZERO_PATTERN) { 2003 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2004 } else { 2005 Mat B; 2006 2007 ierr = MatAXPY_Basic_Preallocate(Y,X,&B);CHKERRQ(ierr); 2008 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2009 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2010 } 2011 PetscFunctionReturn(0); 2012 } 2013 2014 /*MC 2015 MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices 2016 based on the hypre IJ interface. 2017 2018 Level: intermediate 2019 2020 .seealso: MatCreate() 2021 M*/ 2022 2023 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B) 2024 { 2025 Mat_HYPRE *hB; 2026 PetscErrorCode ierr; 2027 2028 PetscFunctionBegin; 2029 ierr = PetscNewLog(B,&hB);CHKERRQ(ierr); 2030 hB->inner_free = PETSC_TRUE; 2031 hB->available = PETSC_TRUE; 2032 hB->size = 0; 2033 hB->array = NULL; 2034 2035 B->data = (void*)hB; 2036 B->rmap->bs = 1; 2037 B->assembled = PETSC_FALSE; 2038 2039 ierr = PetscMemzero(B->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2040 B->ops->mult = MatMult_HYPRE; 2041 B->ops->multtranspose = MatMultTranspose_HYPRE; 2042 B->ops->multadd = MatMultAdd_HYPRE; 2043 B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE; 2044 B->ops->setup = MatSetUp_HYPRE; 2045 B->ops->destroy = MatDestroy_HYPRE; 2046 B->ops->assemblyend = MatAssemblyEnd_HYPRE; 2047 B->ops->assemblybegin = MatAssemblyBegin_HYPRE; 2048 B->ops->ptap = MatPtAP_HYPRE_HYPRE; 2049 B->ops->matmult = MatMatMult_HYPRE_HYPRE; 2050 B->ops->setvalues = MatSetValues_HYPRE; 2051 B->ops->missingdiagonal = MatMissingDiagonal_HYPRE; 2052 B->ops->scale = MatScale_HYPRE; 2053 B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE; 2054 B->ops->zeroentries = MatZeroEntries_HYPRE; 2055 B->ops->zerorows = MatZeroRows_HYPRE; 2056 B->ops->getrow = MatGetRow_HYPRE; 2057 B->ops->restorerow = MatRestoreRow_HYPRE; 2058 B->ops->getvalues = MatGetValues_HYPRE; 2059 B->ops->setoption = MatSetOption_HYPRE; 2060 B->ops->duplicate = MatDuplicate_HYPRE; 2061 B->ops->copy = MatCopy_HYPRE; 2062 B->ops->view = MatView_HYPRE; 2063 B->ops->getdiagonal = MatGetDiagonal_HYPRE; 2064 B->ops->axpy = MatAXPY_HYPRE; 2065 2066 /* build cache for off array entries formed */ 2067 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);CHKERRQ(ierr); 2068 2069 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr); 2070 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr); 2071 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr); 2072 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr); 2073 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr); 2074 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr); 2075 ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);CHKERRQ(ierr); 2076 ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);CHKERRQ(ierr); 2077 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_hypre_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 2078 PetscFunctionReturn(0); 2079 } 2080 2081 static PetscErrorCode hypre_array_destroy(void *ptr) 2082 { 2083 PetscFunctionBegin; 2084 hypre_TFree(ptr,HYPRE_MEMORY_HOST); 2085 PetscFunctionReturn(0); 2086 } 2087