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