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