24 #if defined (HAVE_CONFIG_H)
49 #if ! defined (HAVE_QRUPDATE)
62 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
73 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
75 F77_CHAR_ARG_LEN (1)));
77 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
79 F77_CHAR_ARG_LEN (1)));
111 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
122 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
124 F77_CHAR_ARG_LEN (1)));
126 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
128 F77_CHAR_ARG_LEN (1)));
160 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
168 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
170 F77_CHAR_ARG_LEN (1)));
172 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
174 F77_CHAR_ARG_LEN (1)));
205 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
213 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
215 F77_CHAR_ARG_LEN (1)));
217 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
219 F77_CHAR_ARG_LEN (1)));
245 template <
typename T>
253 template <
typename T>
260 template <
typename T>
264 if (! R.is_square ())
270 #if ! defined (HAVE_QRUPDATE)
272 template <
typename T>
283 init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (),
287 template <
typename T>
289 singular (
const T&
a)
291 static typename T::element_type
zero (0);
293 if (
a(
i,
i) ==
zero)
return true;
297 template <
typename T>
310 if (singular (chol_mat))
314 info = init (chol_mat.hermitian () * chol_mat
315 - T (u) * T (u).hermitian (),
true,
false);
322 template <
typename T>
326 static typename T::element_type
zero (0);
334 if (u.numel () != n + 1)
337 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
339 if (singular (chol_mat))
345 T a = chol_mat.hermitian () * chol_mat;
355 a1(
k, l) =
a(
k < j ?
k :
k-1, l < j ? l : l-1);
357 info = init (a1,
true,
false);
364 template <
typename T>
372 if (j < 0 || j > n-1)
373 (*current_liboctave_error_handler) (
"choldelete: index out of range");
375 T a = chol_mat.hermitian () * chol_mat;
378 init (a,
true,
false);
381 template <
typename T>
390 (*current_liboctave_error_handler) (
"cholshift: index out of range");
392 T a = chol_mat.hermitian () * chol_mat;
421 (*current_liboctave_error_handler) (
"chol: requires square matrix");
428 chol_mat.clear (n, n);
433 chol_mat.xelem (i, j) =
a(i, j);
435 chol_mat.xelem (i, j) = 0.0;
441 chol_mat.xelem (i, j) = 0.0;
443 chol_mat.xelem (i, j) =
a(i, j);
445 double *
h = chol_mat.fortran_vec ();
450 anorm =
xnorm (a, 1);
453 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n, h, n, info
454 F77_CHAR_ARG_LEN (1)));
456 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n, h, n, info
457 F77_CHAR_ARG_LEN (1)));
461 chol_mat.resize (info - 1, info - 1);
472 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n, h,
473 n, anorm, xrcond, pz, piz, dpocon_info
474 F77_CHAR_ARG_LEN (1)));
476 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (
"L", 1), n, h,
477 n, anorm, xrcond, pz, piz, dpocon_info
478 F77_CHAR_ARG_LEN (1)));
480 if (dpocon_info != 0)
487 #if defined (HAVE_QRUPDATE)
502 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
521 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
535 if (u.
numel () != n + 1)
538 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
544 chol_mat.resize (n+1, n+1);
546 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
547 j + 1, utmp.fortran_vec (),
w, info));
558 if (j < 0 || j > n-1)
559 (*current_liboctave_error_handler) (
"choldelete: index out of range");
563 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
566 chol_mat.resize (n-1, n-1);
576 (*current_liboctave_error_handler) (
"cholshift: index out of range");
580 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
594 (*current_liboctave_error_handler) (
"chol: requires square matrix");
601 chol_mat.clear (n, n);
606 chol_mat.xelem (i, j) =
a(i, j);
608 chol_mat.xelem (i, j) = 0.0f;
614 chol_mat.xelem (i, j) = 0.0f;
616 chol_mat.xelem (i, j) =
a(i, j);
618 float *
h = chol_mat.fortran_vec ();
623 anorm =
xnorm (a, 1);
626 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n, h, n, info
627 F77_CHAR_ARG_LEN (1)));
629 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n, h, n, info
630 F77_CHAR_ARG_LEN (1)));
634 chol_mat.resize (info - 1, info - 1);
645 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n, h,
646 n, anorm, xrcond, pz, piz, spocon_info
647 F77_CHAR_ARG_LEN (1)));
649 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (
"L", 1), n, h,
650 n, anorm, xrcond, pz, piz, spocon_info
651 F77_CHAR_ARG_LEN (1)));
653 if (spocon_info != 0)
660 #if defined (HAVE_QRUPDATE)
675 F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
694 F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
708 if (u.
numel () != n + 1)
711 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
717 chol_mat.resize (n+1, n+1);
719 F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
720 j + 1, utmp.fortran_vec (),
w, info));
731 if (j < 0 || j > n-1)
732 (*current_liboctave_error_handler) (
"choldelete: index out of range");
736 F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
739 chol_mat.resize (n-1, n-1);
749 (*current_liboctave_error_handler) (
"cholshift: index out of range");
753 F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
767 (*current_liboctave_error_handler) (
"chol: requires square matrix");
774 chol_mat.clear (n, n);
779 chol_mat.xelem (i, j) =
a(i, j);
781 chol_mat.xelem (i, j) = 0.0;
787 chol_mat.xelem (i, j) = 0.0;
789 chol_mat.xelem (i, j) =
a(i, j);
791 Complex *
h = chol_mat.fortran_vec ();
796 anorm =
xnorm (a, 1);
799 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
801 F77_CHAR_ARG_LEN (1)));
803 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
805 F77_CHAR_ARG_LEN (1)));
809 chol_mat.resize (info - 1, info - 1);
819 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
822 F77_CHAR_ARG_LEN (1)));
824 if (zpocon_info != 0)
831 #if defined (HAVE_QRUPDATE)
882 if (u.
numel () != n + 1)
885 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
891 chol_mat.resize (n+1, n+1);
906 if (j < 0 || j > n-1)
907 (*current_liboctave_error_handler) (
"choldelete: index out of range");
915 chol_mat.resize (n-1, n-1);
925 (*current_liboctave_error_handler) (
"cholshift: index out of range");
946 (*current_liboctave_error_handler) (
"chol: requires square matrix");
953 chol_mat.clear (n, n);
958 chol_mat.xelem (i, j) =
a(i, j);
960 chol_mat.xelem (i, j) = 0.0f;
966 chol_mat.xelem (i, j) = 0.0f;
968 chol_mat.xelem (i, j) =
a(i, j);
975 anorm =
xnorm (a, 1);
980 F77_CHAR_ARG_LEN (1)));
984 F77_CHAR_ARG_LEN (1)));
988 chol_mat.resize (info - 1, info - 1);
1000 F77_CHAR_ARG_LEN (1)));
1002 if (cpocon_info != 0)
1009 #if defined (HAVE_QRUPDATE)
1017 if (u.
numel () != n)
1037 if (u.
numel () != n)
1060 if (u.
numel () != n + 1)
1063 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
1069 chol_mat.resize (n+1, n+1);
1084 if (j < 0 || j > n-1)
1085 (*current_liboctave_error_handler) (
"choldelete: index out of range");
1093 chol_mat.resize (n-1, n-1);
1103 (*current_liboctave_error_handler) (
"cholshift: index out of range");
template FloatComplexMatrix chol2inv< FloatComplexMatrix >(const FloatComplexMatrix &r)
Octave interface to the compression and uncompression libraries.
octave_idx_type init(const T &a, bool upper, bool calc_cond)
octave_idx_type numel(void) const
Number of elements in the array.
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define F77_DBLE_CMPLX_ARG(x)
template Matrix chol2inv< Matrix >(const Matrix &r)
octave_idx_type insert_sym(const VT &u, octave_idx_type j)
#define F77_XFCN(f, F, args)
void warn_qrupdate_once(void)
octave_idx_type rows(void) const
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
ComplexColumnVector conj(const ComplexColumnVector &a)
std::complex< double > w(std::complex< double > z, double relerr=0)
octave_idx_type downdate(const VT &u)
static Matrix chol2inv_internal(const Matrix &r, bool is_upper=true)
void delete_sym(octave_idx_type j)
OCTAVE_API double xnorm(const ColumnVector &x, double p)
template ComplexMatrix chol2inv< ComplexMatrix >(const ComplexMatrix &r)
T & xelem(octave_idx_type n)
=val(i)}if ode{val(i)}occurs in table i
void shift_sym(octave_idx_type i, octave_idx_type j)
template FloatMatrix chol2inv< FloatMatrix >(const FloatMatrix &r)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
std::complex< float > FloatComplex
std::complex< double > Complex
const T * fortran_vec(void) const
octave_idx_type cols(void) const
Vector representing the dimensions (size) of an Array.