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) PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 35 36 PetscCall(PetscMemzero(C, sizeof(*C))); 37 38 C->nrow = (size_t)A->cmap->n; 39 C->ncol = (size_t)A->rmap->n; 40 C->nzmax = (size_t)nz; 41 C->p = ci; 42 C->i = cj; 43 C->x = values ? ca : 0; 44 C->stype = -1; 45 C->itype = CHOLMOD_INT_TYPE; 46 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 47 C->dtype = CHOLMOD_DOUBLE; 48 C->sorted = 1; 49 C->packed = 1; 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode MatFactorGetSolverType_seqaij_cholmod(Mat A, MatSolverType *type) 54 { 55 PetscFunctionBegin; 56 *type = MATSOLVERCHOLMOD; 57 PetscFunctionReturn(PETSC_SUCCESS); 58 } 59 60 /* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */ 61 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A, MatFactorType ftype, Mat *F) 62 { 63 Mat B; 64 Mat_CHOLMOD *chol; 65 PetscInt m = A->rmap->n, n = A->cmap->n; 66 67 PetscFunctionBegin; 68 #if defined(PETSC_USE_COMPLEX) 69 if (A->hermitian != PETSC_BOOL3_TRUE) { 70 PetscCall(PetscInfo(A, "Only for Hermitian matrices.\n")); 71 *F = NULL; 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 #endif 75 /* Create the factorization matrix F */ 76 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 77 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 78 PetscCall(PetscStrallocpy("cholmod", &((PetscObject)B)->type_name)); 79 PetscCall(MatSetUp(B)); 80 PetscCall(PetscNew(&chol)); 81 82 chol->Wrap = MatWrapCholmod_seqaij; 83 B->data = chol; 84 85 B->ops->getinfo = MatGetInfo_CHOLMOD; 86 B->ops->view = MatView_CHOLMOD; 87 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 88 B->ops->destroy = MatDestroy_CHOLMOD; 89 90 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_cholmod)); 91 92 B->factortype = MAT_FACTOR_CHOLESKY; 93 B->assembled = PETSC_TRUE; 94 B->preallocated = PETSC_TRUE; 95 96 PetscCall(PetscFree(B->solvertype)); 97 PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype)); 98 B->canuseordering = PETSC_TRUE; 99 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_CHOLESKY])); 100 PetscCall(CholmodStart(B)); 101 *F = B; 102 PetscFunctionReturn(PETSC_SUCCESS); 103 } 104