1
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>
4
MatWrapCholmod_seqaij(Mat A,PetscBool values,cholmod_sparse * C,PetscBool * aijalloc,PetscBool * valloc)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
MatFactorGetSolverType_seqaij_cholmod(Mat A,MatSolverType * type)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 */
MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat * F)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->useordering = PETSC_TRUE;
102 ierr = CholmodStart(B);CHKERRQ(ierr);
103 *F = B;
104 PetscFunctionReturn(0);
105 }
106