2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>

  7: static PetscErrorCode MatWrapCholmod_seqaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool  *aijalloc)
  8: {
  9:   Mat_SeqAIJ      *aij = (Mat_SeqAIJ*)A->data;
 10:   const PetscInt  *ai  = aij->i,*aj = aij->j,*adiag;
 11:   const MatScalar *aa  = aij->a;
 12:   PetscInt        m    = A->rmap->n,i,j,k,nz,*ci,*cj;
 13:   PetscScalar     *ca;
 14:   PetscErrorCode  ierr;

 17:   MatMarkDiagonal_SeqAIJ(A);
 18:   adiag = aij->diag;
 19:   for (i=0,nz=0; i<m; i++) nz += ai[i+1] - adiag[i];
 20:   PetscMalloc3(m+1,&ci,nz,&cj,values ? nz : 0,&ca);
 21:   for (i=0,k=0; i<m; i++) {
 22:     ci[i] = k;
 23:     for (j=adiag[i]; j<ai[i+1]; j++,k++) {
 24:       cj[k] = aj[j];
 25:       if (values) ca[k] = aa[j];
 26:     }
 27:   }
 28:   ci[i]     = k;
 29:   *aijalloc = PETSC_TRUE;

 31:   PetscMemzero(C,sizeof(*C));

 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:   return(0);
 46: }

 50: static PetscErrorCode MatFactorGetSolverPackage_seqaij_cholmod(Mat A,const MatSolverPackage *type)
 51: {
 53:   *type = MATSOLVERCHOLMOD;
 54:   return(0);
 55: }

 59: /* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */
 60: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat *F)
 61: {
 62:   Mat            B;
 63:   Mat_CHOLMOD    *chol;
 65:   PetscInt       m=A->rmap->n,n=A->cmap->n;

 68:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with AIJ, only %s",
 69:                                              MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]);
 70:   /* Create the factorization matrix F */
 71:   MatCreate(PetscObjectComm((PetscObject)A),&B);
 72:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
 73:   MatSetType(B,((PetscObject)A)->type_name);
 74:   MatSeqAIJSetPreallocation(B,0,NULL);
 75:   PetscNewLog(B,&chol);

 77:   chol->Wrap    = MatWrapCholmod_seqaij;
 78:   chol->Destroy = MatDestroy_SeqAIJ;
 79:   B->spptr      = chol;

 81:   B->ops->matsolve               = NULL;
 82:   B->ops->view                   = MatView_CHOLMOD;
 83:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
 84:   B->ops->destroy                = MatDestroy_CHOLMOD;

 86:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_cholmod);

 88:   B->factortype   = MAT_FACTOR_CHOLESKY;
 89:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
 90:   B->preallocated = PETSC_TRUE;

 92:   PetscFree(B->solvertype);
 93:   PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);

 95:   CholmodStart(B);
 96:   *F   = B;
 97:   return(0);
 98: }