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