1586621ddSJed Brown 2c6db04a5SJed Brown #include <../src/mat/impls/aij/seq/aij.h> 3c6db04a5SJed Brown #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 4586621ddSJed Brown 5586621ddSJed Brown #undef __FUNCT__ 6586621ddSJed Brown #define __FUNCT__ "MatWrapCholmod_seqaij" 7ace3abfcSBarry Smith static PetscErrorCode MatWrapCholmod_seqaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc) 8586621ddSJed Brown { 9586621ddSJed Brown Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 10586621ddSJed Brown const PetscInt *ai = aij->i,*aj = aij->j,*adiag; 11586621ddSJed Brown const MatScalar *aa = aij->a; 12586621ddSJed Brown PetscInt m = A->rmap->n,i,j,k,nz,*ci,*cj; 13586621ddSJed Brown PetscScalar *ca; 14586621ddSJed Brown PetscErrorCode ierr; 15586621ddSJed Brown 16586621ddSJed Brown PetscFunctionBegin; 17586621ddSJed Brown ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 18586621ddSJed Brown adiag = aij->diag; 19586621ddSJed Brown for (i=0,nz=0; i<m; i++) nz += ai[i+1] - adiag[i]; 20*dcca6d9dSJed Brown ierr = PetscMalloc3(m+1,&ci,nz,&cj,values ? nz : 0,&ca);CHKERRQ(ierr); 21586621ddSJed Brown for (i=0,k=0; i<m; i++) { 22586621ddSJed Brown ci[i] = k; 23586621ddSJed Brown for (j=adiag[i]; j<ai[i+1]; j++,k++) { 24586621ddSJed Brown cj[k] = aj[j]; 25586621ddSJed Brown if (values) ca[k] = aa[j]; 26586621ddSJed Brown } 27586621ddSJed Brown } 28586621ddSJed Brown ci[i] = k; 29586621ddSJed Brown *aijalloc = PETSC_TRUE; 30586621ddSJed Brown 31586621ddSJed Brown ierr = PetscMemzero(C,sizeof(*C));CHKERRQ(ierr); 322205254eSKarl Rupp 33586621ddSJed Brown C->nrow = (size_t)A->cmap->n; 34586621ddSJed Brown C->ncol = (size_t)A->rmap->n; 35586621ddSJed Brown C->nzmax = (size_t)nz; 36586621ddSJed Brown C->p = ci; 37586621ddSJed Brown C->i = cj; 38586621ddSJed Brown C->x = values ? ca : 0; 39586621ddSJed Brown C->stype = -1; 40586621ddSJed Brown C->itype = CHOLMOD_INT_TYPE; 41586621ddSJed Brown C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 42586621ddSJed Brown C->dtype = CHOLMOD_DOUBLE; 43586621ddSJed Brown C->sorted = 1; 44586621ddSJed Brown C->packed = 1; 45586621ddSJed Brown PetscFunctionReturn(0); 46586621ddSJed Brown } 47586621ddSJed Brown 48586621ddSJed Brown #undef __FUNCT__ 49586621ddSJed Brown #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_cholmod" 50586621ddSJed Brown PetscErrorCode MatFactorGetSolverPackage_seqaij_cholmod(Mat A,const MatSolverPackage *type) 51586621ddSJed Brown { 52586621ddSJed Brown PetscFunctionBegin; 53586621ddSJed Brown *type = MATSOLVERCHOLMOD; 54586621ddSJed Brown PetscFunctionReturn(0); 55586621ddSJed Brown } 56586621ddSJed Brown 57586621ddSJed Brown #undef __FUNCT__ 58586621ddSJed Brown #define __FUNCT__ "MatGetFactor_seqaij_cholmod" 59586621ddSJed Brown /* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */ 608cc058d9SJed Brown PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 61586621ddSJed Brown { 62586621ddSJed Brown Mat B; 63586621ddSJed Brown Mat_CHOLMOD *chol; 64586621ddSJed Brown PetscErrorCode ierr; 65586621ddSJed Brown PetscInt m=A->rmap->n,n=A->cmap->n; 66586621ddSJed Brown 67586621ddSJed Brown PetscFunctionBegin; 68586621ddSJed Brown if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with AIJ, only %s", 69586621ddSJed Brown MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 70586621ddSJed Brown /* Create the factorization matrix F */ 71ce94432eSBarry Smith ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 72586621ddSJed Brown ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 73586621ddSJed Brown ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 740298fd71SBarry Smith ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 75586621ddSJed Brown ierr = PetscNewLog(B,Mat_CHOLMOD,&chol);CHKERRQ(ierr); 762205254eSKarl Rupp 77586621ddSJed Brown chol->Wrap = MatWrapCholmod_seqaij; 78586621ddSJed Brown chol->Destroy = MatDestroy_SeqAIJ; 79586621ddSJed Brown B->spptr = chol; 80586621ddSJed Brown 81586621ddSJed Brown B->ops->view = MatView_CHOLMOD; 82586621ddSJed Brown B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 83586621ddSJed Brown B->ops->destroy = MatDestroy_CHOLMOD; 842205254eSKarl Rupp 85bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_cholmod);CHKERRQ(ierr); 862205254eSKarl Rupp 87586621ddSJed Brown B->factortype = MAT_FACTOR_CHOLESKY; 88586621ddSJed Brown B->assembled = PETSC_TRUE; /* required by -ksp_view */ 89586621ddSJed Brown B->preallocated = PETSC_TRUE; 90586621ddSJed Brown 91586621ddSJed Brown ierr = CholmodStart(B);CHKERRQ(ierr); 92586621ddSJed Brown *F = B; 93586621ddSJed Brown PetscFunctionReturn(0); 94586621ddSJed Brown } 95