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