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,PetscBool *valloc) 6 { 7 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 8 const PetscScalar *aa; 9 PetscScalar *ca; 10 const PetscInt *ai = aij->i,*aj = aij->j,*adiag; 11 PetscInt m = A->rmap->n,i,j,k,nz,*ci,*cj; 12 PetscBool vain = PETSC_FALSE; 13 14 PetscFunctionBegin; 15 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 16 adiag = aij->diag; 17 for (i=0,nz=0; i<m; i++) nz += ai[i+1] - adiag[i]; 18 PetscCall(PetscMalloc2(m+1,&ci,nz,&cj)); 19 if (values) { 20 vain = PETSC_TRUE; 21 PetscCall(PetscMalloc1(nz,&ca)); 22 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 23 } 24 for (i=0,k=0; i<m; i++) { 25 ci[i] = k; 26 for (j=adiag[i]; j<ai[i+1]; j++,k++) { 27 cj[k] = aj[j]; 28 if (values) ca[k] = PetscConj(aa[j]); 29 } 30 } 31 ci[i] = k; 32 *aijalloc = PETSC_TRUE; 33 *valloc = vain; 34 if (values) { 35 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 36 } 37 38 PetscCall(PetscMemzero(C,sizeof(*C))); 39 40 C->nrow = (size_t)A->cmap->n; 41 C->ncol = (size_t)A->rmap->n; 42 C->nzmax = (size_t)nz; 43 C->p = ci; 44 C->i = cj; 45 C->x = values ? ca : 0; 46 C->stype = -1; 47 C->itype = CHOLMOD_INT_TYPE; 48 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 49 C->dtype = CHOLMOD_DOUBLE; 50 C->sorted = 1; 51 C->packed = 1; 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode MatFactorGetSolverType_seqaij_cholmod(Mat A,MatSolverType *type) 56 { 57 PetscFunctionBegin; 58 *type = MATSOLVERCHOLMOD; 59 PetscFunctionReturn(0); 60 } 61 62 /* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */ 63 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 64 { 65 Mat B; 66 Mat_CHOLMOD *chol; 67 PetscInt m=A->rmap->n,n=A->cmap->n; 68 const char *prefix; 69 70 PetscFunctionBegin; 71 #if defined(PETSC_USE_COMPLEX) 72 PetscCheck(A->hermitian,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for hermitian matrices"); 73 #endif 74 /* Create the factorization matrix F */ 75 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 76 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 77 PetscCall(PetscStrallocpy("cholmod",&((PetscObject)B)->type_name)); 78 PetscCall(MatGetOptionsPrefix(A,&prefix)); 79 PetscCall(MatSetOptionsPrefix(B,prefix)); 80 PetscCall(MatSetUp(B)); 81 PetscCall(PetscNewLog(B,&chol)); 82 83 chol->Wrap = MatWrapCholmod_seqaij; 84 B->data = chol; 85 86 B->ops->getinfo = MatGetInfo_CHOLMOD; 87 B->ops->view = MatView_CHOLMOD; 88 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 89 B->ops->destroy = MatDestroy_CHOLMOD; 90 91 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_cholmod)); 92 93 B->factortype = MAT_FACTOR_CHOLESKY; 94 B->assembled = PETSC_TRUE; 95 B->preallocated = PETSC_TRUE; 96 97 PetscCall(PetscFree(B->solvertype)); 98 PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype)); 99 B->canuseordering = PETSC_TRUE; 100 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_CHOLESKY])); 101 PetscCall(CholmodStart(B)); 102 *F = B; 103 PetscFunctionReturn(0); 104 } 105