1 2 /* 3 Creates hypre ijmatrix from PETSc matrix 4 */ 5 6 /*MC 7 MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices 8 based on the hypre IJ interface. 9 10 Level: intermediate 11 12 .seealso: MatCreate() 13 M*/ 14 15 #include <petscmathypre.h> 16 #include <petsc/private/matimpl.h> 17 #include <../src/mat/impls/hypre/mhypre.h> 18 #include <../src/mat/impls/aij/mpi/mpiaij.h> 19 #include <../src/vec/vec/impls/hypre/vhyp.h> 20 #include <HYPRE.h> 21 #include <HYPRE_utilities.h> 22 #include <_hypre_parcsr_ls.h> 23 24 static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*); 25 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix); 26 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix); 27 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix); 28 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,Vec,Vec,PetscBool); 29 static PetscErrorCode hypre_array_destroy(void*); 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "MatHYPRE_IJMatrixPreallocate" 33 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij) 34 { 35 PetscErrorCode ierr; 36 PetscInt i,n_d,n_o; 37 const PetscInt *ia_d,*ia_o; 38 PetscBool done_d=PETSC_FALSE,done_o=PETSC_FALSE; 39 PetscInt *nnz_d=NULL,*nnz_o=NULL; 40 41 PetscFunctionBegin; 42 if (A_d) { /* determine number of nonzero entries in local diagonal part */ 43 ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);CHKERRQ(ierr); 44 if (done_d) { 45 ierr = PetscMalloc1(n_d,&nnz_d);CHKERRQ(ierr); 46 for (i=0; i<n_d; i++) { 47 nnz_d[i] = ia_d[i+1] - ia_d[i]; 48 } 49 } 50 ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);CHKERRQ(ierr); 51 } 52 if (A_o) { /* determine number of nonzero entries in local off-diagonal part */ 53 ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr); 54 if (done_o) { 55 ierr = PetscMalloc1(n_o,&nnz_o);CHKERRQ(ierr); 56 for (i=0; i<n_o; i++) { 57 nnz_o[i] = ia_o[i+1] - ia_o[i]; 58 } 59 } 60 ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);CHKERRQ(ierr); 61 } 62 if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */ 63 if (!done_o) { /* only diagonal part */ 64 ierr = PetscMalloc1(n_d,&nnz_o);CHKERRQ(ierr); 65 for (i=0; i<n_d; i++) { 66 nnz_o[i] = 0; 67 } 68 } 69 PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o)); 70 ierr = PetscFree(nnz_d);CHKERRQ(ierr); 71 ierr = PetscFree(nnz_o);CHKERRQ(ierr); 72 } 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "MatHYPRE_CreateFromMat" 78 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA) 79 { 80 PetscErrorCode ierr; 81 PetscInt rstart,rend,cstart,cend; 82 83 PetscFunctionBegin; 84 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 85 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 86 rstart = A->rmap->rstart; 87 rend = A->rmap->rend; 88 cstart = A->cmap->rstart; 89 cend = A->cmap->rend; 90 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij)); 91 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 92 { 93 PetscBool same; 94 Mat A_d,A_o; 95 const PetscInt *colmap; 96 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr); 97 if (same) { 98 ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 99 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr); 103 if (same) { 104 ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 105 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr); 109 if (same) { 110 ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr); 114 if (same) { 115 ierr = MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 } 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatHYPRE_IJMatrixCopy" 124 static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij) 125 { 126 PetscErrorCode ierr; 127 PetscInt i,rstart,rend,ncols,nr,nc; 128 const PetscScalar *values; 129 const PetscInt *cols; 130 PetscBool flg; 131 132 PetscFunctionBegin; 133 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 134 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 135 if (flg && nr == nc) { 136 ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 140 if (flg) { 141 ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 146 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 147 for (i=rstart; i<rend; i++) { 148 ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 149 PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values)); 150 ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 151 } 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_SeqAIJ" 157 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij) 158 { 159 PetscErrorCode ierr; 160 Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data; 161 HYPRE_Int type; 162 hypre_ParCSRMatrix *par_matrix; 163 hypre_AuxParCSRMatrix *aux_matrix; 164 hypre_CSRMatrix *hdiag; 165 166 PetscFunctionBegin; 167 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 168 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type)); 169 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 170 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix)); 171 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 172 /* 173 this is the Hack part where we monkey directly with the hypre datastructures 174 */ 175 ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt)); 176 ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt)); 177 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar)); 178 179 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 180 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_MPIAIJ" 186 static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij) 187 { 188 PetscErrorCode ierr; 189 Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data; 190 Mat_SeqAIJ *pdiag,*poffd; 191 PetscInt i,*garray = pA->garray,*jj,cstart,*pjj; 192 HYPRE_Int type; 193 hypre_ParCSRMatrix *par_matrix; 194 hypre_AuxParCSRMatrix *aux_matrix; 195 hypre_CSRMatrix *hdiag,*hoffd; 196 197 PetscFunctionBegin; 198 pdiag = (Mat_SeqAIJ*) pA->A->data; 199 poffd = (Mat_SeqAIJ*) pA->B->data; 200 /* cstart is only valid for square MPIAIJ layed out in the usual way */ 201 ierr = MatGetOwnershipRange(A,&cstart,NULL);CHKERRQ(ierr); 202 203 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij)); 204 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type)); 205 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 206 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix)); 207 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 208 hoffd = hypre_ParCSRMatrixOffd(par_matrix); 209 210 /* 211 this is the Hack part where we monkey directly with the hypre datastructures 212 */ 213 ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 214 /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */ 215 jj = (PetscInt*)hdiag->j; 216 pjj = (PetscInt*)pdiag->j; 217 for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i]; 218 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar)); 219 ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 220 /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this 221 If we hacked a hypre a bit more we might be able to avoid this step */ 222 jj = (PetscInt*) hoffd->j; 223 pjj = (PetscInt*) poffd->j; 224 for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]]; 225 ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar)); 226 227 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 228 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatConvert_HYPRE_IS" 234 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B) 235 { 236 Mat_HYPRE* mhA = (Mat_HYPRE*)(A->data); 237 Mat lA; 238 ISLocalToGlobalMapping rl2g,cl2g; 239 IS is; 240 hypre_ParCSRMatrix *hA; 241 hypre_CSRMatrix *hdiag,*hoffd; 242 MPI_Comm comm; 243 PetscScalar *hdd,*hod,*aa,*data; 244 HYPRE_Int *col_map_offd,*hdi,*hdj,*hoi,*hoj; 245 PetscInt *ii,*jj,*iptr,*jptr; 246 PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N; 247 HYPRE_Int type; 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 comm = PetscObjectComm((PetscObject)A); 252 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type)); 253 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 254 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA)); 255 M = hypre_ParCSRMatrixGlobalNumRows(hA); 256 N = hypre_ParCSRMatrixGlobalNumCols(hA); 257 str = hypre_ParCSRMatrixFirstRowIndex(hA); 258 stc = hypre_ParCSRMatrixFirstColDiag(hA); 259 hdiag = hypre_ParCSRMatrixDiag(hA); 260 hoffd = hypre_ParCSRMatrixOffd(hA); 261 dr = hypre_CSRMatrixNumRows(hdiag); 262 dc = hypre_CSRMatrixNumCols(hdiag); 263 nnz = hypre_CSRMatrixNumNonzeros(hdiag); 264 hdi = hypre_CSRMatrixI(hdiag); 265 hdj = hypre_CSRMatrixJ(hdiag); 266 hdd = hypre_CSRMatrixData(hdiag); 267 oc = hypre_CSRMatrixNumCols(hoffd); 268 nnz += hypre_CSRMatrixNumNonzeros(hoffd); 269 hoi = hypre_CSRMatrixI(hoffd); 270 hoj = hypre_CSRMatrixJ(hoffd); 271 hod = hypre_CSRMatrixData(hoffd); 272 if (reuse != MAT_REUSE_MATRIX) { 273 PetscInt *aux; 274 275 /* generate l2g maps for rows and cols */ 276 ierr = ISCreateStride(comm,dr,str,1,&is);CHKERRQ(ierr); 277 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 278 ierr = ISDestroy(&is);CHKERRQ(ierr); 279 col_map_offd = hypre_ParCSRMatrixColMapOffd(hA); 280 ierr = PetscMalloc1(dc+oc,&aux);CHKERRQ(ierr); 281 for (i=0; i<dc; i++) aux[i] = i+stc; 282 for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i]; 283 ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 284 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 285 ierr = ISDestroy(&is);CHKERRQ(ierr); 286 /* create MATIS object */ 287 ierr = MatCreate(comm,B);CHKERRQ(ierr); 288 ierr = MatSetSizes(*B,dr,dc,M,N);CHKERRQ(ierr); 289 ierr = MatSetType(*B,MATIS);CHKERRQ(ierr); 290 ierr = MatSetLocalToGlobalMapping(*B,rl2g,cl2g);CHKERRQ(ierr); 291 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 292 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 293 294 /* allocate CSR for local matrix */ 295 ierr = PetscMalloc1(dr+1,&iptr);CHKERRQ(ierr); 296 ierr = PetscMalloc1(nnz,&jptr);CHKERRQ(ierr); 297 ierr = PetscMalloc1(nnz,&data);CHKERRQ(ierr); 298 } else { 299 PetscInt nr; 300 PetscBool done; 301 ierr = MatISGetLocalMat(*B,&lA);CHKERRQ(ierr); 302 ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);CHKERRQ(ierr); 303 if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr); 304 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); 305 ierr = MatSeqAIJGetArray(lA,&data);CHKERRQ(ierr); 306 } 307 /* merge local matrices */ 308 ii = iptr; 309 jj = jptr; 310 aa = data; 311 *ii = *(hdi++) + *(hoi++); 312 for (jd=0,jo=0,cum=0; *ii<nnz; cum++) { 313 PetscScalar *aold = aa; 314 PetscInt *jold = jj,nc = jd+jo; 315 for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; } 316 for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; } 317 *(++ii) = *(hdi++) + *(hoi++); 318 ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);CHKERRQ(ierr); 319 } 320 for (; cum<dr; cum++) *(++ii) = nnz; 321 if (reuse != MAT_REUSE_MATRIX) { 322 Mat_SeqAIJ* a; 323 324 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);CHKERRQ(ierr); 325 ierr = MatISSetLocalMat(*B,lA);CHKERRQ(ierr); 326 /* hack SeqAIJ */ 327 a = (Mat_SeqAIJ*)(lA->data); 328 a->free_a = PETSC_TRUE; 329 a->free_ij = PETSC_TRUE; 330 ierr = MatDestroy(&lA);CHKERRQ(ierr); 331 } 332 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 333 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 334 if (reuse == MAT_INPLACE_MATRIX) { 335 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 336 } 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "MatConvert_AIJ_HYPRE" 342 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B) 343 { 344 Mat_HYPRE *hB; 345 MPI_Comm comm = PetscObjectComm((PetscObject)A); 346 PetscErrorCode ierr; 347 348 PetscFunctionBegin; 349 if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX"); 350 if (reuse == MAT_REUSE_MATRIX) { 351 /* always destroy the old matrix and create a new memory; 352 hope this does not churn the memory too much. The problem 353 is I do not know if it is possible to put the matrix back to 354 its initial state so that we can directly copy the values 355 the second time through. */ 356 hB = (Mat_HYPRE*)((*B)->data); 357 PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij)); 358 } else { 359 ierr = MatCreate(comm,B);CHKERRQ(ierr); 360 ierr = MatSetType(*B,MATHYPRE);CHKERRQ(ierr); 361 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 362 hB = (Mat_HYPRE*)((*B)->data); 363 } 364 ierr = MatHYPRE_CreateFromMat(A,hB);CHKERRQ(ierr); 365 ierr = MatHYPRE_IJMatrixCopy(A,hB->ij);CHKERRQ(ierr); 366 (*B)->preallocated = PETSC_TRUE; 367 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 368 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 369 PetscFunctionReturn(0); 370 } 371 372 #undef __FUNCT__ 373 #define __FUNCT__ "MatConvert_HYPRE_AIJ" 374 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B) 375 { 376 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 377 hypre_ParCSRMatrix *parcsr; 378 hypre_CSRMatrix *hdiag,*hoffd; 379 MPI_Comm comm; 380 PetscScalar *da,*oa,*aptr; 381 PetscInt *dii,*djj,*oii,*ojj,*iptr; 382 PetscInt i,dnnz,onnz,m,n; 383 HYPRE_Int type; 384 PetscMPIInt size; 385 PetscErrorCode ierr; 386 387 PetscFunctionBegin; 388 comm = PetscObjectComm((PetscObject)A); 389 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 390 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 391 if (reuse == MAT_REUSE_MATRIX) { 392 PetscBool ismpiaij,isseqaij; 393 ierr = PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 394 ierr = PetscObjectTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 395 if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported"); 396 } 397 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 398 399 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr)); 400 hdiag = hypre_ParCSRMatrixDiag(parcsr); 401 hoffd = hypre_ParCSRMatrixOffd(parcsr); 402 m = hypre_CSRMatrixNumRows(hdiag); 403 n = hypre_CSRMatrixNumCols(hdiag); 404 dnnz = hypre_CSRMatrixNumNonzeros(hdiag); 405 onnz = hypre_CSRMatrixNumNonzeros(hoffd); 406 if (reuse == MAT_INITIAL_MATRIX) { 407 ierr = PetscMalloc1(m+1,&dii);CHKERRQ(ierr); 408 ierr = PetscMalloc1(dnnz,&djj);CHKERRQ(ierr); 409 ierr = PetscMalloc1(dnnz,&da);CHKERRQ(ierr); 410 } else if (reuse == MAT_REUSE_MATRIX) { 411 PetscInt nr; 412 PetscBool done; 413 if (size > 1) { 414 Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data); 415 416 ierr = MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr); 417 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); 418 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); 419 ierr = MatSeqAIJGetArray(b->A,&da);CHKERRQ(ierr); 420 } else { 421 ierr = MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);CHKERRQ(ierr); 422 if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m); 423 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); 424 ierr = MatSeqAIJGetArray(*B,&da);CHKERRQ(ierr); 425 } 426 } else { /* MAT_INPLACE_MATRIX */ 427 dii = (PetscInt*)hypre_CSRMatrixI(hdiag); 428 djj = (PetscInt*)hypre_CSRMatrixJ(hdiag); 429 da = hypre_CSRMatrixData(hdiag); 430 } 431 ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 432 ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));CHKERRQ(ierr); 433 ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));CHKERRQ(ierr); 434 iptr = djj; 435 aptr = da; 436 for (i=0; i<m; i++) { 437 PetscInt nc = dii[i+1]-dii[i]; 438 ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr); 439 iptr += nc; 440 aptr += nc; 441 } 442 if (size > 1) { 443 HYPRE_Int *offdj,*coffd; 444 445 if (reuse == MAT_INITIAL_MATRIX) { 446 ierr = PetscMalloc1(m+1,&oii);CHKERRQ(ierr); 447 ierr = PetscMalloc1(onnz,&ojj);CHKERRQ(ierr); 448 ierr = PetscMalloc1(onnz,&oa);CHKERRQ(ierr); 449 } else if (reuse == MAT_REUSE_MATRIX) { 450 Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data); 451 PetscInt nr,hr = hypre_CSRMatrixNumRows(hoffd); 452 PetscBool done; 453 454 ierr = MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);CHKERRQ(ierr); 455 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); 456 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); 457 ierr = MatSeqAIJGetArray(b->B,&oa);CHKERRQ(ierr); 458 } else { /* MAT_INPLACE_MATRIX */ 459 oii = (PetscInt*)hypre_CSRMatrixI(hoffd); 460 ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd); 461 oa = hypre_CSRMatrixData(hoffd); 462 } 463 ierr = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 464 offdj = hypre_CSRMatrixJ(hoffd); 465 coffd = hypre_ParCSRMatrixColMapOffd(parcsr); 466 for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]]; 467 ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));CHKERRQ(ierr); 468 iptr = ojj; 469 aptr = oa; 470 for (i=0; i<m; i++) { 471 PetscInt nc = oii[i+1]-oii[i]; 472 ierr = PetscSortIntWithScalarArray(nc,iptr,aptr);CHKERRQ(ierr); 473 iptr += nc; 474 aptr += nc; 475 } 476 if (reuse == MAT_INITIAL_MATRIX) { 477 Mat_MPIAIJ *b; 478 Mat_SeqAIJ *d,*o; 479 480 ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);CHKERRQ(ierr); 481 /* hack MPIAIJ */ 482 b = (Mat_MPIAIJ*)((*B)->data); 483 d = (Mat_SeqAIJ*)b->A->data; 484 o = (Mat_SeqAIJ*)b->B->data; 485 d->free_a = PETSC_TRUE; 486 d->free_ij = PETSC_TRUE; 487 o->free_a = PETSC_TRUE; 488 o->free_ij = PETSC_TRUE; 489 } else if (reuse == MAT_INPLACE_MATRIX) { 490 Mat T; 491 ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);CHKERRQ(ierr); 492 hypre_CSRMatrixI(hdiag) = NULL; 493 hypre_CSRMatrixJ(hdiag) = NULL; 494 hypre_CSRMatrixData(hdiag) = NULL; 495 hypre_CSRMatrixI(hoffd) = NULL; 496 hypre_CSRMatrixJ(hoffd) = NULL; 497 hypre_CSRMatrixData(hoffd) = NULL; 498 ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr); 499 } 500 } else { 501 oii = NULL; 502 ojj = NULL; 503 oa = NULL; 504 if (reuse == MAT_INITIAL_MATRIX) { 505 Mat_SeqAIJ* b; 506 ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);CHKERRQ(ierr); 507 /* hack SeqAIJ */ 508 b = (Mat_SeqAIJ*)((*B)->data); 509 b->free_a = PETSC_TRUE; 510 b->free_ij = PETSC_TRUE; 511 } else if (reuse == MAT_INPLACE_MATRIX) { 512 Mat T; 513 ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);CHKERRQ(ierr); 514 hypre_CSRMatrixI(hdiag) = NULL; 515 hypre_CSRMatrixJ(hdiag) = NULL; 516 hypre_CSRMatrixData(hdiag) = NULL; 517 ierr = MatHeaderReplace(A,&T);CHKERRQ(ierr); 518 } 519 } 520 521 /* we have to use hypre_Tfree to free the arrays */ 522 if (reuse == MAT_INPLACE_MATRIX) { 523 void *ptrs[6] = {dii,djj,da,oii,ojj,oa}; 524 const char *names[6] = {"_hypre_csr_dii", 525 "_hypre_csr_djj", 526 "_hypre_csr_da", 527 "_hypre_csr_oii", 528 "_hypre_csr_ojj", 529 "_hypre_csr_oa"}; 530 for (i=0; i<6; i++) { 531 PetscContainer c; 532 533 ierr = PetscContainerCreate(comm,&c);CHKERRQ(ierr); 534 ierr = PetscContainerSetPointer(c,ptrs[i]);CHKERRQ(ierr); 535 ierr = PetscContainerSetUserDestroy(c,hypre_array_destroy);CHKERRQ(ierr); 536 ierr = PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);CHKERRQ(ierr); 537 ierr = PetscContainerDestroy(&c);CHKERRQ(ierr); 538 } 539 } 540 PetscFunctionReturn(0); 541 } 542 543 #undef __FUNCT__ 544 #define __FUNCT__ "MatPtAP_AIJ_HYPRE" 545 static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C) 546 { 547 Mat_HYPRE *hP = (Mat_HYPRE*)P->data; 548 hypre_ParCSRMatrix *tA,*Pparcsr,*ptapparcsr; 549 hypre_CSRMatrix *hdiag,*hoffd; 550 Mat_SeqAIJ *diag,*offd; 551 PetscInt *garray,noffd,dnnz,onnz,*row_starts,*col_starts; 552 HYPRE_Int type,P_owns_col_starts; 553 PetscBool ismpiaij,isseqaij; 554 MPI_Comm comm = PetscObjectComm((PetscObject)A); 555 char mtype[256]; 556 PetscErrorCode ierr; 557 558 PetscFunctionBegin; 559 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type)); 560 if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 561 if (scall == MAT_REUSE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_REUSE_MATRIX"); 562 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 563 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 564 if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type); 565 566 /* It looks like we don't need to have the diagonal entries 567 ordered first in the rows of the diagonal part 568 for boomerAMGBuildCoarseOperator to work */ 569 if (ismpiaij) { 570 Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data); 571 572 diag = (Mat_SeqAIJ*)a->A->data; 573 offd = (Mat_SeqAIJ*)a->B->data; 574 garray = a->garray; 575 noffd = a->B->cmap->N; 576 dnnz = diag->nz; 577 onnz = offd->nz; 578 } else { 579 diag = (Mat_SeqAIJ*)A->data; 580 offd = NULL; 581 garray = NULL; 582 noffd = 0; 583 dnnz = diag->nz; 584 onnz = 0; 585 } 586 587 /* create a temporary ParCSR */ 588 if (HYPRE_AssumedPartitionCheck()) { 589 PetscMPIInt myid; 590 591 ierr = MPI_Comm_rank(comm,&myid);CHKERRQ(ierr); 592 row_starts = A->rmap->range + myid; 593 col_starts = A->cmap->range + myid; 594 } else { 595 row_starts = A->rmap->range; 596 col_starts = A->cmap->range; 597 } 598 tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_Int*)row_starts,(HYPRE_Int*)col_starts,noffd,dnnz,onnz); 599 hypre_ParCSRMatrixSetRowStartsOwner(tA,0); 600 hypre_ParCSRMatrixSetColStartsOwner(tA,0); 601 602 /* set diagonal part */ 603 hdiag = hypre_ParCSRMatrixDiag(tA); 604 hypre_CSRMatrixI(hdiag) = (HYPRE_Int*)diag->i; 605 hypre_CSRMatrixJ(hdiag) = (HYPRE_Int*)diag->j; 606 hypre_CSRMatrixData(hdiag) = diag->a; 607 hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz; 608 hypre_CSRMatrixSetRownnz(hdiag); 609 hypre_CSRMatrixSetDataOwner(hdiag,0); 610 611 /* set offdiagonal part */ 612 hoffd = hypre_ParCSRMatrixOffd(tA); 613 if (offd) { 614 hypre_CSRMatrixI(hoffd) = (HYPRE_Int*)offd->i; 615 hypre_CSRMatrixJ(hoffd) = (HYPRE_Int*)offd->j; 616 hypre_CSRMatrixData(hoffd) = offd->a; 617 hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz; 618 hypre_CSRMatrixSetRownnz(hoffd); 619 hypre_CSRMatrixSetDataOwner(hoffd,0); 620 hypre_ParCSRMatrixSetNumNonzeros(tA); 621 hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_Int*)garray; 622 } 623 624 /* call RAP from BoomerAMG */ 625 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts 626 from Pparcsr (even if it does not own them)! */ 627 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr)); 628 P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(Pparcsr); 629 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 630 PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(Pparcsr,tA,Pparcsr,&ptapparcsr)); 631 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 632 hypre_ParCSRMatrixSetRowStartsOwner(ptapparcsr,0); 633 hypre_ParCSRMatrixSetColStartsOwner(ptapparcsr,0); 634 if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(Pparcsr, 1); 635 636 /* set pointers to NULL before destroying tA */ 637 hypre_CSRMatrixI(hdiag) = NULL; 638 hypre_CSRMatrixJ(hdiag) = NULL; 639 hypre_CSRMatrixData(hdiag) = NULL; 640 hypre_CSRMatrixI(hoffd) = NULL; 641 hypre_CSRMatrixJ(hoffd) = NULL; 642 hypre_CSRMatrixData(hoffd) = NULL; 643 hypre_ParCSRMatrixColMapOffd(tA) = NULL; 644 hypre_ParCSRMatrixDestroy(tA); 645 646 /* create C depending on mtype */ 647 sprintf(mtype,MATAIJ); 648 ierr = PetscOptionsGetString(((PetscObject)A)->options,((PetscObject)A)->prefix,"-matptap_hypre_outtype",mtype,256,NULL);CHKERRQ(ierr); 649 ierr = MatCreateFromParCSR(ptapparcsr,mtype,PETSC_OWN_POINTER,C);CHKERRQ(ierr); 650 PetscFunctionReturn(0); 651 } 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "MatPtAP_HYPRE_HYPRE" 655 static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C) 656 { 657 hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr; 658 Mat_HYPRE *hA = (Mat_HYPRE*)A->data, *hP = (Mat_HYPRE*)P->data; 659 HYPRE_Int type,P_owns_col_starts; 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 if (scall == MAT_REUSE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Unsupported MAT_REUSE_MATRIX"); 664 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 665 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 666 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type)); 667 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 668 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr)); 669 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr)); 670 671 /* call RAP from BoomerAMG */ 672 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts 673 from Pparcsr (even if it does not own them)! */ 674 P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(Pparcsr); 675 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 676 PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr)); 677 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 678 hypre_ParCSRMatrixSetRowStartsOwner(ptapparcsr,0); 679 hypre_ParCSRMatrixSetColStartsOwner(ptapparcsr,0); 680 if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(Pparcsr, 1); 681 682 /* create MatHYPRE */ 683 ierr = MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,C);CHKERRQ(ierr); 684 PetscFunctionReturn(0); 685 } 686 687 #undef __FUNCT__ 688 #define __FUNCT__ "MatMultTranspose_HYPRE" 689 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y) 690 { 691 PetscErrorCode ierr; 692 693 PetscFunctionBegin; 694 ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 #undef __FUNCT__ 699 #define __FUNCT__ "MatMult_HYPRE" 700 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y) 701 { 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 ierr = MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNCT__ 710 #define __FUNCT__ "MatHYPRE_MultKernel_Private" 711 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans) 712 { 713 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 714 hypre_ParCSRMatrix *parcsr; 715 hypre_ParVector *hx,*hy; 716 PetscScalar *ax,*ay,*sax,*say; 717 PetscErrorCode ierr; 718 719 PetscFunctionBegin; 720 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr)); 721 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx)); 722 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy)); 723 ierr = VecGetArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr); 724 ierr = VecGetArray(y,&ay);CHKERRQ(ierr); 725 if (trans) { 726 VecHYPRE_ParVectorReplacePointer(hA->x,ay,say); 727 VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax); 728 hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx); 729 VecHYPRE_ParVectorReplacePointer(hA->x,say,ay); 730 VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax); 731 } else { 732 VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax); 733 VecHYPRE_ParVectorReplacePointer(hA->b,ay,say); 734 hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy); 735 VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax); 736 VecHYPRE_ParVectorReplacePointer(hA->b,say,ay); 737 } 738 ierr = VecRestoreArrayRead(x,(const PetscScalar**)&ax);CHKERRQ(ierr); 739 ierr = VecRestoreArray(y,&ay);CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743 #undef __FUNCT__ 744 #define __FUNCT__ "MatDestroy_HYPRE" 745 static PetscErrorCode MatDestroy_HYPRE(Mat A) 746 { 747 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 748 PetscErrorCode ierr; 749 750 PetscFunctionBegin; 751 if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x)); 752 if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b)); 753 if (hA->ij) { 754 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL; 755 PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij)); 756 } 757 if (hA->comm) { ierr = MPI_Comm_free(&hA->comm);CHKERRQ(ierr); } 758 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);CHKERRQ(ierr); 759 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);CHKERRQ(ierr); 760 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);CHKERRQ(ierr); 761 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);CHKERRQ(ierr); 762 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);CHKERRQ(ierr); 763 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);CHKERRQ(ierr); 764 ierr = PetscFree(A->data);CHKERRQ(ierr); 765 PetscFunctionReturn(0); 766 } 767 768 #undef __FUNCT__ 769 #define __FUNCT__ "MatSetUp_HYPRE" 770 static PetscErrorCode MatSetUp_HYPRE(Mat A) 771 { 772 PetscErrorCode ierr; 773 774 PetscFunctionBegin; 775 ierr = MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "MatAssemblyEnd_HYPRE" 781 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode) 782 { 783 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 784 Vec x,b; 785 PetscErrorCode ierr; 786 787 PetscFunctionBegin; 788 if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE"); 789 PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij)); 790 if (hA->x) PetscFunctionReturn(0); 791 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 792 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 793 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);CHKERRQ(ierr); 794 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);CHKERRQ(ierr); 795 ierr = VecHYPRE_IJVectorCreate(x,&hA->x);CHKERRQ(ierr); 796 ierr = VecHYPRE_IJVectorCreate(b,&hA->b);CHKERRQ(ierr); 797 ierr = VecDestroy(&x);CHKERRQ(ierr); 798 ierr = VecDestroy(&b);CHKERRQ(ierr); 799 PetscFunctionReturn(0); 800 } 801 802 #define MATHYPRE_SCRATCH 2048 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "MatSetValues_HYPRE" 806 PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins) 807 { 808 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 809 PetscScalar *vals = (PetscScalar *)v; 810 PetscScalar sscr[MATHYPRE_SCRATCH]; 811 HYPRE_Int cscr[2][MATHYPRE_SCRATCH]; 812 HYPRE_Int i,nzc; 813 PetscErrorCode ierr; 814 815 PetscFunctionBegin; 816 for (i=0,nzc=0;i<nc;i++) { 817 if (cols[i] >= 0) { 818 cscr[0][nzc ] = cols[i]; 819 cscr[1][nzc++] = i; 820 } 821 } 822 if (!nzc) PetscFunctionReturn(0); 823 824 if (ins == ADD_VALUES) { 825 for (i=0;i<nr;i++) { 826 if (rows[i] >= 0) { 827 PetscInt j; 828 for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]]; 829 PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&nzc,(HYPRE_Int*)(rows+i),cscr[0],sscr)); 830 } 831 vals += nc; 832 } 833 } else { /* INSERT_VALUES */ 834 #if defined(PETSC_USE_DEBUG) 835 /* Insert values cannot be used to insert offproc entries */ 836 PetscInt rst,ren; 837 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr); 838 for (i=0;i<nr;i++) 839 if (rows[i] >= 0 && (rows[i] < rst || rows[i] >= ren)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use INSERT_VALUES for off-proc entries with MatHYPRE. Use ADD_VALUES instead"); 840 #endif 841 for (i=0;i<nr;i++) { 842 if (rows[i] >= 0) { 843 PetscInt j; 844 for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]]; 845 PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&nzc,(HYPRE_Int*)(rows+i),cscr[0],sscr)); 846 } 847 vals += nc; 848 } 849 } 850 PetscFunctionReturn(0); 851 } 852 853 #undef __FUNCT__ 854 #define __FUNCT__ "MatHYPRESetPreallocation_HYPRE" 855 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 856 { 857 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 858 HYPRE_Int *hdnnz,*honnz; 859 PetscInt i,rs,re,cs,ce,bs; 860 PetscMPIInt size; 861 PetscErrorCode ierr; 862 863 PetscFunctionBegin; 864 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 865 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 866 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 867 rs = A->rmap->rstart; 868 re = A->rmap->rend; 869 cs = A->cmap->rstart; 870 ce = A->cmap->rend; 871 if (!hA->ij) { 872 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij)); 873 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 874 } else { 875 HYPRE_Int hrs,hre,hcs,hce; 876 PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce)); 877 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); 878 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); 879 } 880 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij)); 881 882 if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs; 883 if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs; 884 885 if (!dnnz) { 886 ierr = PetscMalloc1(A->rmap->n,&hdnnz);CHKERRQ(ierr); 887 for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz; 888 } else { 889 hdnnz = (HYPRE_Int*)dnnz; 890 } 891 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 892 if (size > 1) { 893 if (!onnz) { 894 ierr = PetscMalloc1(A->rmap->n,&honnz);CHKERRQ(ierr); 895 for (i=0;i<A->rmap->n;i++) honnz[i] = onz; 896 } else { 897 honnz = (HYPRE_Int*)onnz; 898 } 899 PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz)); 900 } else { 901 honnz = NULL; 902 PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz)); 903 } 904 if (!dnnz) { 905 ierr = PetscFree(hdnnz);CHKERRQ(ierr); 906 } 907 if (!onnz && honnz) { 908 ierr = PetscFree(honnz);CHKERRQ(ierr); 909 } 910 A->preallocated = PETSC_TRUE; 911 912 /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0 */ 913 { 914 hypre_AuxParCSRMatrix *aux_matrix; 915 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij); 916 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; 917 } 918 PetscFunctionReturn(0); 919 } 920 921 /*@C 922 MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format 923 924 Collective on Mat 925 926 Input Parameters: 927 + A - the matrix 928 . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix 929 (same value is used for all local rows) 930 . dnnz - array containing the number of nonzeros in the various rows of the 931 DIAGONAL portion of the local submatrix (possibly different for each row) 932 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 933 The size of this array is equal to the number of local rows, i.e 'm'. 934 For matrices that will be factored, you must leave room for (and set) 935 the diagonal entry even if it is zero. 936 . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local 937 submatrix (same value is used for all local rows). 938 - onnz - array containing the number of nonzeros in the various rows of the 939 OFF-DIAGONAL portion of the local submatrix (possibly different for 940 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 941 structure. The size of this array is equal to the number 942 of local rows, i.e 'm'. 943 944 Notes: If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored. 945 946 Level: intermediate 947 948 .keywords: matrix, aij, compressed row, sparse, parallel 949 950 .seealso: MatCreate(), MatMPIAIJSetPreallocation, MATHYPRE 951 @*/ 952 #undef __FUNCT__ 953 #define __FUNCT__ "MatHYPRESetPreallocation" 954 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 955 { 956 PetscErrorCode ierr; 957 958 PetscFunctionBegin; 959 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 960 PetscValidType(A,1); 961 ierr = PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 /* 966 MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix 967 968 Collective 969 970 Input Parameters: 971 + vparcsr - the pointer to the hypre_ParCSRMatrix 972 . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported. 973 - copymode - PETSc copying options 974 975 Output Parameter: 976 . A - the matrix 977 978 Level: intermediate 979 980 .seealso: MatHYPRE, PetscCopyMode 981 */ 982 #undef __FUNCT__ 983 #define __FUNCT__ "MatCreateFromParCSR" 984 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *vparcsr, MatType mtype, PetscCopyMode copymode, Mat* A) 985 { 986 Mat T; 987 Mat_HYPRE *hA; 988 hypre_ParCSRMatrix *parcsr; 989 MPI_Comm comm; 990 PetscInt rstart,rend,cstart,cend,M,N; 991 PetscBool isseqaij,ismpiaij,isaij,ishyp,isis; 992 PetscErrorCode ierr; 993 994 PetscFunctionBegin; 995 parcsr = (hypre_ParCSRMatrix *)vparcsr; 996 comm = hypre_ParCSRMatrixComm(parcsr); 997 ierr = PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 998 ierr = PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 999 ierr = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr); 1000 ierr = PetscStrcmp(mtype,MATHYPRE,&ishyp);CHKERRQ(ierr); 1001 ierr = PetscStrcmp(mtype,MATIS,&isis);CHKERRQ(ierr); 1002 isaij = (PetscBool)(isseqaij || ismpiaij || isaij); 1003 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); 1004 if (ishyp && copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES"); 1005 1006 /* access ParCSRMatrix */ 1007 rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr); 1008 rend = hypre_ParCSRMatrixLastRowIndex(parcsr); 1009 cstart = hypre_ParCSRMatrixFirstColDiag(parcsr); 1010 cend = hypre_ParCSRMatrixLastColDiag(parcsr); 1011 M = hypre_ParCSRMatrixGlobalNumRows(parcsr); 1012 N = hypre_ParCSRMatrixGlobalNumCols(parcsr); 1013 1014 /* create PETSc matrix with MatHYPRE */ 1015 ierr = MatCreate(comm,&T);CHKERRQ(ierr); 1016 ierr = MatSetSizes(T,rend-rstart+1,cend-cstart+1,M,N);CHKERRQ(ierr); 1017 ierr = MatSetType(T,MATHYPRE);CHKERRQ(ierr); 1018 hA = (Mat_HYPRE*)(T->data); 1019 1020 /* create HYPRE_IJMatrix */ 1021 PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij)); 1022 1023 /* set ParCSR object */ 1024 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR)); 1025 hypre_IJMatrixObject(hA->ij) = parcsr; 1026 T->preallocated = PETSC_TRUE; 1027 1028 /* set assembled flag */ 1029 hypre_IJMatrixAssembleFlag(hA->ij) = 1; 1030 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij)); 1031 if (ishyp) { 1032 /* prevent from freeing the pointer */ 1033 if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE; 1034 *A = T; 1035 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1036 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1037 } else if (isaij) { 1038 if (copymode != PETSC_OWN_POINTER) { 1039 /* prevent from freeing the pointer */ 1040 hA->inner_free = PETSC_FALSE; 1041 ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr); 1042 ierr = MatDestroy(&T);CHKERRQ(ierr); 1043 } else { /* AIJ return type with PETSC_OWN_POINTER */ 1044 ierr = MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr); 1045 *A = T; 1046 } 1047 } else if (isis) { 1048 ierr = MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);CHKERRQ(ierr); 1049 if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE; 1050 ierr = MatDestroy(&T);CHKERRQ(ierr); 1051 } 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "MatHYPREGetParCSR_HYPRE" 1057 PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr) 1058 { 1059 Mat_HYPRE* hA = (Mat_HYPRE*)A->data; 1060 HYPRE_Int type; 1061 PetscErrorCode ierr; 1062 1063 PetscFunctionBegin; 1064 if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present"); 1065 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type)); 1066 if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR"); 1067 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr)); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 /* 1072 MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix 1073 1074 Not collective 1075 1076 Input Parameters: 1077 + A - the MATHYPRE object 1078 1079 Output Parameter: 1080 . parcsr - the pointer to the hypre_ParCSRMatrix 1081 1082 Level: intermediate 1083 1084 .seealso: MatHYPRE, PetscCopyMode 1085 */ 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "MatHYPREGetParCSR" 1088 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr) 1089 { 1090 PetscErrorCode ierr; 1091 1092 PetscFunctionBegin; 1093 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1094 PetscValidType(A,1); 1095 ierr = PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));CHKERRQ(ierr); 1096 PetscFunctionReturn(0); 1097 } 1098 1099 #undef __FUNCT__ 1100 #define __FUNCT__ "MatCreate_HYPRE" 1101 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B) 1102 { 1103 Mat_HYPRE *hB; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 ierr = PetscNewLog(B,&hB);CHKERRQ(ierr); 1108 hB->inner_free = PETSC_TRUE; 1109 1110 B->data = (void*)hB; 1111 B->rmap->bs = 1; 1112 B->assembled = PETSC_FALSE; 1113 1114 B->ops->mult = MatMult_HYPRE; 1115 B->ops->multtranspose = MatMultTranspose_HYPRE; 1116 B->ops->setup = MatSetUp_HYPRE; 1117 B->ops->destroy = MatDestroy_HYPRE; 1118 B->ops->assemblyend = MatAssemblyEnd_HYPRE; 1119 B->ops->ptap = MatPtAP_HYPRE_HYPRE; 1120 B->ops->setvalues = MatSetValues_HYPRE; 1121 1122 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);CHKERRQ(ierr); 1123 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);CHKERRQ(ierr); 1124 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);CHKERRQ(ierr); 1125 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);CHKERRQ(ierr); 1126 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr); 1127 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);CHKERRQ(ierr); 1128 ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);CHKERRQ(ierr); 1129 ierr = PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);CHKERRQ(ierr); 1130 PetscFunctionReturn(0); 1131 } 1132 1133 #undef __FUNCT__ 1134 #define __FUNCT__ "hypre_array_destroy" 1135 static PetscErrorCode hypre_array_destroy(void *ptr) 1136 { 1137 PetscFunctionBegin; 1138 hypre_TFree(ptr); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #if 0 1143 /* 1144 Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format 1145 1146 This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first 1147 which will corrupt the PETSc data structure if we did this. Need a work around to this problem. 1148 */ 1149 #include <_hypre_IJ_mv.h> 1150 #include <HYPRE_IJ_mv.h> 1151 1152 #undef __FUNCT__ 1153 #define __FUNCT__ "MatHYPRE_IJMatrixLink" 1154 PetscErrorCode MatHYPRE_IJMatrixLink(Mat A, HYPRE_IJMatrix *ij) 1155 { 1156 PetscErrorCode ierr; 1157 PetscInt rstart,rend,cstart,cend; 1158 PetscBool flg; 1159 hypre_AuxParCSRMatrix *aux_matrix; 1160 1161 PetscFunctionBegin; 1162 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1163 PetscValidType(A,1); 1164 PetscValidPointer(ij,2); 1165 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 1166 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices"); 1167 ierr = MatSetUp(A);CHKERRQ(ierr); 1168 1169 ierr = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 1170 rstart = A->rmap->rstart; 1171 rend = A->rmap->rend; 1172 cstart = A->cmap->rstart; 1173 cend = A->cmap->rend; 1174 PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij)); 1175 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR)); 1176 1177 PetscStackCallStandard(HYPRE_IJMatrixInitialize,(*ij)); 1178 PetscStackCall("hypre_IJMatrixTranslator",aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij)); 1179 1180 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 1181 1182 /* this is the Hack part where we monkey directly with the hypre datastructures */ 1183 1184 PetscStackCallStandard(HYPRE_IJMatrixAssemble,(*ij)); 1185 ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 1186 PetscFunctionReturn(0); 1187 } 1188 #endif 1189