38 template <
class chol_type,
class chol_elt,
class p_type>
41 (
const cholmod_sparse* S)
53 Sx =
static_cast<chol_elt *
>(S->x);
57 for (k = 0; k < ncol; k++)
65 if (CHOLMOD_IS_NONZERO (sik))
80 template <
class chol_type,
class chol_elt,
class p_type>
83 (
const chol_type& a,
bool natural,
bool force)
93 (*current_liboctave_error_handler)
94 (
"SparseCHOL requires square matrix");
98 cholmod_common *cm = &Common;
101 CHOLMOD_NAME(start) (cm);
102 cm->prefer_zomplex =
false;
108 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
112 cm->print =
static_cast<int> (spu) + 2;
113 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
117 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
118 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
120 cm->final_asis =
false;
121 cm->final_super =
false;
123 cm->final_pack =
true;
124 cm->final_monotonic =
true;
125 cm->final_resymbol =
false;
128 cholmod_sparse *ac = &
A;
135 ac->nzmax = a.nnz ();
139 #ifdef USE_64_BIT_IDX_T
140 ac->itype = CHOLMOD_LONG;
142 ac->itype = CHOLMOD_INT;
144 ac->dtype = CHOLMOD_DOUBLE;
146 #ifdef OCTAVE_CHOLMOD_TYPE
149 ac->xtype = CHOLMOD_REAL;
161 cm->method[0].ordering = CHOLMOD_NATURAL ;
162 cm->postorder = false ;
165 cholmod_factor *Lfactor;
167 Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
168 CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
171 is_pd = cm->status == CHOLMOD_OK;
172 info = (is_pd ? 0 : cm->status);
177 cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
180 minor_p = Lfactor->minor;
183 Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm);
186 if (minor_p > 0 && minor_p < a_nr)
188 size_t n1 = a_nr + 1;
189 Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
191 Lsparse->p, &n1, cm);
193 CHOLMOD_NAME(reallocate_sparse)
196 Lsparse->ncol = minor_p;
199 drop_zeros (Lsparse);
208 static char tmp[] =
" ";
211 CHOLMOD_NAME(free_factor) (&Lfactor, cm);
212 CHOLMOD_NAME(finish) (cm);
213 CHOLMOD_NAME(print_common) (tmp, cm);
217 (*current_liboctave_error_handler)
218 (
"Missing CHOLMOD. Sparse cholesky factorization disabled");
223 template <
class chol_type,
class chol_elt,
class p_type>
228 cholmod_sparse *m = rep->L ();
231 chol_type ret (m->nrow, nc, nnz);
237 ret.xdata (i) =
static_cast<chol_elt *
>(m->x)[i];
245 template <
class chol_type,
class chol_elt,
class p_type>
268 template <
class chol_type,
class chol_elt,
class p_type>
274 cholmod_sparse *m = rep->L ();
281 chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
286 retval = Qc * linv * linv.hermitian () * Qc.transpose ();
289 retval = linv * linv.hermitian ();
#define OCTAVE_CHOLMOD_TYPE
void drop_zeros(const cholmod_sparse *S)
void SparseCholError(int status, char *file, int line, char *message)
F77_RET_T const octave_idx_type Complex * A
static double get_key(const std::string &key)
chol_type inverse(void) const
octave_idx_type init(const chol_type &a, bool natural, bool force)
octave_idx_type length(void) const
Number of elements in the array.
int SparseCholPrint(const char *fmt,...)
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double * Q
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE