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