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