25 #if defined (HAVE_CONFIG_H)
41 template <
typename chol_type>
47 : count (1), is_pd (
false), minor_p (0), perms (), cond (0)
49 , Lsparse (0), Common ()
54 : count (1), is_pd (
false), minor_p (0), perms (), cond (0)
56 , Lsparse (0), Common ()
59 init (a, natural, force);
63 bool natural,
bool force)
64 : count (1), is_pd (
false), minor_p (0), perms (), cond (0)
66 , Lsparse (0), Common ()
69 info = init (a, natural, force);
74 #if defined (HAVE_CHOLMOD)
76 CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
78 CHOLMOD_NAME(finish) (&Common);
82 #if defined (HAVE_CHOLMOD)
83 cholmod_sparse *
L (
void)
const
91 #if defined (HAVE_CHOLMOD)
92 return (minor_p == static_cast<octave_idx_type>(Lsparse->ncol) ?
105 double rcond (
void)
const {
return cond; }
119 #if defined (HAVE_CHOLMOD)
124 void drop_zeros (
const cholmod_sparse *S);
136 #if defined (HAVE_CHOLMOD)
141 template <
typename chol_type>
161 for (; p < pend; p++)
165 if (CHOLMOD_IS_NONZERO (sik))
182 template <
typename T>
197 return CHOLMOD_COMPLEX;
202 template <
typename chol_type>
205 bool natural,
bool force)
209 #if defined (HAVE_CHOLMOD)
215 (*current_liboctave_error_handler) (
"sparse_chol requires square matrix");
217 cholmod_common *cm = &Common;
221 CHOLMOD_NAME(
start) (cm);
222 cm->prefer_zomplex =
false;
229 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
233 cm->print =
static_cast<int> (spu) + 2;
234 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
239 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
240 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
242 cm->final_asis =
false;
243 cm->final_super =
false;
245 cm->final_pack =
true;
246 cm->final_monotonic =
true;
247 cm->final_resymbol =
false;
250 cholmod_sparse *ac = &
A;
258 ac->nzmax = a.nnz ();
262 #if defined (OCTAVE_ENABLE_64)
263 ac->itype = CHOLMOD_LONG;
265 ac->itype = CHOLMOD_INT;
267 ac->dtype = CHOLMOD_DOUBLE;
269 ac->xtype = get_xtype<chol_elt> ();
280 cm->method[0].ordering = CHOLMOD_NATURAL;
281 cm->postorder =
false;
284 cholmod_factor *Lfactor;
286 Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
287 CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
290 is_pd = cm->status == CHOLMOD_OK;
291 info = (is_pd ? 0 : cm->status);
296 cond = CHOLMOD_NAME(
rcond) (Lfactor, cm);
299 minor_p = Lfactor->minor;
302 Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm);
305 if (minor_p > 0 && minor_p < a_nr)
307 size_t n1 = a_nr + 1;
308 Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
310 Lsparse->p, &n1, cm);
312 CHOLMOD_NAME(reallocate_sparse)
316 Lsparse->ncol = minor_p;
319 drop_zeros (Lsparse);
330 static char blank_name[] =
" ";
333 CHOLMOD_NAME(print_common) (blank_name, cm);
334 CHOLMOD_NAME(free_factor) (&Lfactor, cm);
341 octave_unused_parameter (a);
342 octave_unused_parameter (natural);
343 octave_unused_parameter (force);
345 (*current_liboctave_error_handler)
346 (
"support for CHOLMOD was unavailable or disabled when liboctave was built");
353 template <
typename chol_type>
357 #if defined (HAVE_CHOLMOD)
380 template <
typename chol_type>
385 template <
typename chol_type>
392 template <
typename chol_type>
394 bool natural,
bool force)
399 template <
typename chol_type>
406 template <
typename chol_type>
412 template <
typename chol_type>
419 template <
typename chol_type>
422 if (--rep->count == 0)
426 template <
typename chol_type>
432 if (--rep->count == 0)
442 template <
typename chol_type>
446 #if defined (HAVE_CHOLMOD)
448 cholmod_sparse *
m = rep->L ();
453 chol_type ret (m->nrow, nc, nnz);
461 ret.xdata (
i) =
static_cast<chol_elt *
>(m->x)[
i];
473 template <
typename chol_type>
480 template <
typename chol_type>
487 template <
typename chol_type>
494 template <
typename chol_type>
498 return rep->is_positive_definite ();
501 template <
typename chol_type>
505 return rep->rcond ();
508 template <
typename chol_type>
514 #if defined (HAVE_CHOLMOD)
516 cholmod_sparse *
m = rep->L ();
522 chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
524 if (perms.
numel () == n)
538 template <
typename chol_type>
547 (*current_liboctave_error_handler) (
"U must be a square matrix");
550 int typ = mattype.type (
false);
553 chol_type rtra, multip;
557 rtra = r.transpose ();
562 rtra = r.transpose ();
569 retval = multip.inverse (mattypenew, info, rcond,
true,
false);
octave_idx_type * xridx(void)
octave_idx_type P(void) const
sparse_chol_rep(const chol_type &a, bool natural, bool force)
Octave interface to the compression and uncompression libraries.
SparseMatrix Q(void) const
SparseMatrix transpose(void) const
octave_idx_type numel(void) const
Number of elements in the array.
chol_type::element_type chol_elt
bool is_positive_definite(void) const
void SparseCholError(int status, char *file, int line, char *message)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
RowVector perm(void) const
octave_idx_type * xcidx(void)
chol_type inverse(void) const
SparseMatrix hermitian(void) const
sparse_chol_rep(const chol_type &a, octave_idx_type &info, bool natural, bool force)
static double get_key(const std::string &key)
int get_xtype< Complex >(void)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
octave_idx_type init(const chol_type &a, bool natural, bool force)
if(nargin< 2) print_usage()
nd deftypefn *octave_map m
sparse_chol & operator=(const sparse_chol &a)
bool is_positive_definite(void) const
template SparseMatrix chol2inv< SparseMatrix >(const SparseMatrix &r)
cholmod_sparse * L(void) const
SparseMatrix Q(void) const
void drop_zeros(const cholmod_sparse *S)
int get_xtype< double >(void)
defaults to zero A value of zero computes the digamma a value the trigamma and so on The digamma function is defined
int SparseCholPrint(const char *fmt,...)
octave_refcount< int > count
virtual ~sparse_chol(void)
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
=val(i)}if ode{val(i)}occurs in table i
RowVector perm(void) const
template SparseComplexMatrix chol2inv< SparseComplexMatrix >(const SparseComplexMatrix &r)
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
octave_idx_type P(void) const
F77_RET_T const F77_INT F77_CMPLX * A
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q