1 #define PETSCMAT_DLL 2 3 #include "../src/mat/impls/aij/seq/aij.h" 4 #include "../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h" 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "MatWrapCholmod_seqaij" 8 static PetscErrorCode MatWrapCholmod_seqaij(Mat A,PetscTruth values,cholmod_sparse *C,PetscTruth *aijalloc) 9 { 10 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 11 const PetscInt *ai = aij->i,*aj = aij->j,*adiag; 12 const MatScalar *aa = aij->a; 13 PetscInt m = A->rmap->n,i,j,k,nz,*ci,*cj; 14 PetscScalar *ca; 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 19 adiag = aij->diag; 20 for (i=0,nz=0; i<m; i++) nz += ai[i+1] - adiag[i]; 21 ierr = PetscMalloc3(m+1,PetscInt,&ci,nz,PetscInt,&cj,values?nz:0,PetscScalar,&ca);CHKERRQ(ierr); 22 for (i=0,k=0; i<m; i++) { 23 ci[i] = k; 24 for (j=adiag[i]; j<ai[i+1]; j++,k++) { 25 cj[k] = aj[j]; 26 if (values) ca[k] = aa[j]; 27 } 28 } 29 ci[i] = k; 30 *aijalloc = PETSC_TRUE; 31 32 ierr = PetscMemzero(C,sizeof(*C));CHKERRQ(ierr); 33 C->nrow = (size_t)A->cmap->n; 34 C->ncol = (size_t)A->rmap->n; 35 C->nzmax = (size_t)nz; 36 C->p = ci; 37 C->i = cj; 38 C->x = values ? ca : 0; 39 C->stype = -1; 40 C->itype = CHOLMOD_INT_TYPE; 41 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 42 C->dtype = CHOLMOD_DOUBLE; 43 C->sorted = 1; 44 C->packed = 1; 45 PetscFunctionReturn(0); 46 } 47 48 EXTERN_C_BEGIN 49 #undef __FUNCT__ 50 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_cholmod" 51 PetscErrorCode MatFactorGetSolverPackage_seqaij_cholmod(Mat A,const MatSolverPackage *type) 52 { 53 PetscFunctionBegin; 54 *type = MATSOLVERCHOLMOD; 55 PetscFunctionReturn(0); 56 } 57 EXTERN_C_END 58 59 EXTERN_C_BEGIN 60 #undef __FUNCT__ 61 #define __FUNCT__ "MatGetFactor_seqaij_cholmod" 62 /* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */ 63 PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 64 { 65 Mat B; 66 Mat_CHOLMOD *chol; 67 PetscErrorCode ierr; 68 PetscInt m=A->rmap->n,n=A->cmap->n; 69 70 PetscFunctionBegin; 71 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with AIJ, only %s", 72 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 73 /* Create the factorization matrix F */ 74 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 75 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 76 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 77 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 78 ierr = PetscNewLog(B,Mat_CHOLMOD,&chol);CHKERRQ(ierr); 79 chol->Wrap = MatWrapCholmod_seqaij; 80 chol->Destroy = MatDestroy_SeqAIJ; 81 B->spptr = chol; 82 83 B->ops->view = MatView_CHOLMOD; 84 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 85 B->ops->destroy = MatDestroy_CHOLMOD; 86 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_cholmod",MatFactorGetSolverPackage_seqaij_cholmod);CHKERRQ(ierr); 87 B->factortype = MAT_FACTOR_CHOLESKY; 88 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 89 B->preallocated = PETSC_TRUE; 90 91 ierr = CholmodStart(B);CHKERRQ(ierr); 92 *F = B; 93 PetscFunctionReturn(0); 94 } 95 EXTERN_C_END 96