24 #if defined (HAVE_CONFIG_H) 44 #if ! defined (HAVE_QRUPDATE) 57 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
59 F77_INT n = octave::to_f77_int (r_nc);
63 double *v =
tmp.fortran_vec ();
68 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
70 F77_CHAR_ARG_LEN (1)));
72 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
74 F77_CHAR_ARG_LEN (1)));
106 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
108 F77_INT n = octave::to_f77_int (r_nc);
112 float *v =
tmp.fortran_vec ();
117 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
119 F77_CHAR_ARG_LEN (1)));
121 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
123 F77_CHAR_ARG_LEN (1)));
133 tmp.xelem (
i, j) =
tmp.xelem (j,
i);
137 tmp.xelem (j,
i) =
tmp.xelem (
i, j);
155 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
157 F77_INT n = octave::to_f77_int (r_nc);
163 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
165 F77_CHAR_ARG_LEN (1)));
167 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
169 F77_CHAR_ARG_LEN (1)));
200 (*current_liboctave_error_handler) (
"chol2inv requires square matrix");
202 F77_INT n = octave::to_f77_int (r_nc);
208 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
210 F77_CHAR_ARG_LEN (1)));
212 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
214 F77_CHAR_ARG_LEN (1)));
240 template <
typename T>
248 template <
typename T>
255 template <
typename T>
265 #if ! defined (HAVE_QRUPDATE) 267 template <
typename T>
278 init (chol_mat.hermitian () * chol_mat + T (
u) * T (
u).hermitian (),
282 template <
typename T>
284 singular (
const T&
a)
286 static typename T::element_type
zero (0);
288 if (
a(
i,
i) ==
zero)
return true;
292 template <
typename T>
305 if (singular (chol_mat))
309 info = init (chol_mat.hermitian () * chol_mat
310 - T (
u) * T (
u).hermitian (),
true,
false);
317 template <
typename T>
321 static typename T::element_type
zero (0);
329 if (
u.numel () != n + 1)
332 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
334 if (singular (chol_mat))
340 T
a = chol_mat.hermitian () * chol_mat;
350 a1(
k, l) =
a(
k < j ?
k :
k-1, l < j ? l : l-1);
352 info = init (a1,
true,
false);
359 template <
typename T>
367 if (j < 0 || j > n-1)
368 (*current_liboctave_error_handler) (
"choldelete: index out of range");
370 T
a = chol_mat.hermitian () * chol_mat;
373 init (
a,
true,
false);
376 template <
typename T>
384 if (i < 0 || i > n-1 || j < 0 || j > n-1)
385 (*current_liboctave_error_handler) (
"cholshift: index out of range");
387 T
a = chol_mat.hermitian () * chol_mat;
416 (*current_liboctave_error_handler) (
"chol: requires square matrix");
423 chol_mat.clear (n, n);
428 chol_mat.xelem (
i, j) =
a(
i, j);
430 chol_mat.xelem (
i, j) = 0.0;
436 chol_mat.xelem (
i, j) = 0.0;
438 chol_mat.xelem (
i, j) =
a(
i, j);
440 double *
h = chol_mat.fortran_vec ();
448 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
h, n, info
449 F77_CHAR_ARG_LEN (1)));
451 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
h, n, info
452 F77_CHAR_ARG_LEN (1)));
456 chol_mat.resize (info - 1, info - 1);
466 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
h,
467 n, anorm, xrcond, pz, iz, dpocon_info
468 F77_CHAR_ARG_LEN (1)));
470 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
h,
471 n, anorm, xrcond, pz, iz, dpocon_info
472 F77_CHAR_ARG_LEN (1)));
474 if (dpocon_info != 0)
481 #if defined (HAVE_QRUPDATE) 487 F77_INT n = to_f77_int (chol_mat.rows ());
496 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), n,
506 F77_INT n = to_f77_int (chol_mat.rows ());
515 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), n,
527 F77_INT n = to_f77_int (chol_mat.rows ());
528 F77_INT j = to_f77_int (j_arg);
530 if (
u.numel () != n + 1)
533 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
539 chol_mat.resize (n+1, n+1);
540 F77_INT ldcm = to_f77_int (chol_mat.rows ());
542 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), ldcm,
543 j + 1, utmp.fortran_vec (),
w, info));
552 F77_INT n = to_f77_int (chol_mat.rows ());
553 F77_INT j = to_f77_int (j_arg);
555 if (j < 0 || j > n-1)
556 (*current_liboctave_error_handler) (
"choldelete: index out of range");
560 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), n, j + 1,
w));
562 chol_mat.resize (n-1, n-1);
569 F77_INT n = to_f77_int (chol_mat.rows ());
571 F77_INT j = to_f77_int (j_arg);
573 if (i < 0 || i > n-1 || j < 0 || j > n-1)
574 (*current_liboctave_error_handler) (
"cholshift: index out of range");
578 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), n,
592 (*current_liboctave_error_handler) (
"chol: requires square matrix");
599 chol_mat.clear (n, n);
604 chol_mat.xelem (
i, j) =
a(
i, j);
606 chol_mat.xelem (
i, j) = 0.0f;
612 chol_mat.xelem (
i, j) = 0.0f;
614 chol_mat.xelem (
i, j) =
a(
i, j);
616 float *
h = chol_mat.fortran_vec ();
624 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
h, n, info
625 F77_CHAR_ARG_LEN (1)));
627 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
h, n, info
628 F77_CHAR_ARG_LEN (1)));
632 chol_mat.resize (info - 1, info - 1);
642 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
h,
643 n, anorm, xrcond, pz, iz, spocon_info
644 F77_CHAR_ARG_LEN (1)));
646 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
h,
647 n, anorm, xrcond, pz, iz, spocon_info
648 F77_CHAR_ARG_LEN (1)));
650 if (spocon_info != 0)
657 #if defined (HAVE_QRUPDATE) 663 F77_INT n = to_f77_int (chol_mat.rows ());
672 F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), n,
682 F77_INT n = to_f77_int (chol_mat.rows ());
691 F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), n,
704 F77_INT n = to_f77_int (chol_mat.rows ());
705 F77_INT j = to_f77_int (j_arg);
707 if (
u.numel () != n + 1)
710 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
716 chol_mat.resize (n+1, n+1);
717 F77_INT ldcm = to_f77_int (chol_mat.rows ());
719 F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), ldcm,
720 j + 1, utmp.fortran_vec (),
w, info));
729 F77_INT n = to_f77_int (chol_mat.rows ());
730 F77_INT j = to_f77_int (j_arg);
732 if (j < 0 || j > n-1)
733 (*current_liboctave_error_handler) (
"choldelete: index out of range");
737 F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), n,
740 chol_mat.resize (n-1, n-1);
747 F77_INT n = to_f77_int (chol_mat.rows ());
749 F77_INT j = to_f77_int (j_arg);
751 if (i < 0 || i > n-1 || j < 0 || j > n-1)
752 (*current_liboctave_error_handler) (
"cholshift: index out of range");
756 F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), n,
770 (*current_liboctave_error_handler) (
"chol: requires square matrix");
777 chol_mat.clear (n, n);
782 chol_mat.xelem (
i, j) =
a(
i, j);
784 chol_mat.xelem (
i, j) = 0.0;
790 chol_mat.xelem (
i, j) = 0.0;
792 chol_mat.xelem (
i, j) =
a(
i, j);
794 Complex *
h = chol_mat.fortran_vec ();
802 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
804 F77_CHAR_ARG_LEN (1)));
806 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1), n,
808 F77_CHAR_ARG_LEN (1)));
812 chol_mat.resize (info - 1, info - 1);
822 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
825 F77_CHAR_ARG_LEN (1)));
827 if (zpocon_info != 0)
834 #if defined (HAVE_QRUPDATE) 840 F77_INT n = to_f77_int (chol_mat.rows ());
862 F77_INT n = to_f77_int (chol_mat.rows ());
887 F77_INT n = to_f77_int (chol_mat.rows ());
888 F77_INT j = to_f77_int (j_arg);
890 if (
u.numel () != n + 1)
893 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
899 chol_mat.resize (n+1, n+1);
900 F77_INT ldcm = to_f77_int (chol_mat.rows ());
915 F77_INT n = to_f77_int (chol_mat.rows ());
916 F77_INT j = to_f77_int (j_arg);
918 if (j < 0 || j > n-1)
919 (*current_liboctave_error_handler) (
"choldelete: index out of range");
927 chol_mat.resize (n-1, n-1);
935 F77_INT n = to_f77_int (chol_mat.rows ());
937 F77_INT j = to_f77_int (j_arg);
939 if (i < 0 || i > n-1 || j < 0 || j > n-1)
940 (*current_liboctave_error_handler) (
"cholshift: index out of range");
962 (*current_liboctave_error_handler) (
"chol: requires square matrix");
969 chol_mat.clear (n, n);
974 chol_mat.xelem (
i, j) =
a(
i, j);
976 chol_mat.xelem (
i, j) = 0.0f;
982 chol_mat.xelem (
i, j) = 0.0f;
984 chol_mat.xelem (
i, j) =
a(
i, j);
994 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (
"U", 1),
996 F77_CHAR_ARG_LEN (1)));
998 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1),
1000 F77_CHAR_ARG_LEN (1)));
1004 chol_mat.resize (info - 1, info - 1);
1014 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (
"U", 1), n,
1017 F77_CHAR_ARG_LEN (1)));
1019 if (cpocon_info != 0)
1026 #if defined (HAVE_QRUPDATE) 1032 F77_INT n = to_f77_int (chol_mat.rows ());
1034 if (
u.numel () != n)
1051 F77_INT n = to_f77_int (chol_mat.rows ());
1053 if (
u.numel () != n)
1073 F77_INT j = to_f77_int (j_arg);
1075 F77_INT n = to_f77_int (chol_mat.rows ());
1077 if (
u.numel () != n + 1)
1080 (*current_liboctave_error_handler) (
"cholinsert: index out of range");
1086 chol_mat.resize (n+1, n+1);
1087 F77_INT ldcm = to_f77_int (chol_mat.rows ());
1101 F77_INT n = to_f77_int (chol_mat.rows ());
1102 F77_INT j = to_f77_int (j_arg);
1104 if (j < 0 || j > n-1)
1105 (*current_liboctave_error_handler) (
"choldelete: index out of range");
1112 chol_mat.resize (n-1, n-1);
1120 F77_INT n = to_f77_int (chol_mat.rows ());
1122 F77_INT j = to_f77_int (j_arg);
1124 if (i < 0 || i > n-1 || j < 0 || j > n-1)
1125 (*current_liboctave_error_handler) (
"cholshift: index out of range");
octave_idx_type rows(void) const
template FloatComplexMatrix chol2inv< FloatComplexMatrix >(const FloatComplexMatrix &r)
octave_idx_type init(const T &a, bool upper, bool calc_cond)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by zero($0/0$)
#define F77_DBLE_CMPLX_ARG(x)
const T * fortran_vec(void) const
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 cols(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)
octave_f77_int_type F77_INT
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)
ColumnVector imag(const ComplexColumnVector &a)
std::complex< float > FloatComplex
std::complex< double > Complex
Vector representing the dimensions (size) of an Array.