Mercurial > octave
diff liboctave/numeric/chol.cc @ 22329:7f3c7a8bd131
maint: Indent namespaces in liboctave/numeric files.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 17 Aug 2016 10:55:38 -0400 |
parents | bac0d6f07a3e |
children | 4caa7b28d183 |
line wrap: on
line diff
--- a/liboctave/numeric/chol.cc Wed Aug 17 10:37:57 2016 -0400 +++ b/liboctave/numeric/chol.cc Wed Aug 17 10:55:38 2016 -0400 @@ -240,885 +240,883 @@ namespace octave { -namespace math -{ - -template <typename T> -T -chol2inv (const T& r) -{ - return chol2inv_internal (r); -} + namespace math + { + template <typename T> + T + chol2inv (const T& r) + { + return chol2inv_internal (r); + } -// Compute the inverse of a matrix using the Cholesky factorization. -template <typename T> -T -chol<T>::inverse (void) const -{ - return chol2inv_internal (chol_mat, is_upper); -} + // Compute the inverse of a matrix using the Cholesky factorization. + template <typename T> + T + chol<T>::inverse (void) const + { + return chol2inv_internal (chol_mat, is_upper); + } -template <typename T> -void -chol<T>::set (const T& R) -{ - if (! R.is_square ()) - (*current_liboctave_error_handler) ("chol: requires square matrix"); + template <typename T> + void + chol<T>::set (const T& R) + { + if (! R.is_square ()) + (*current_liboctave_error_handler) ("chol: requires square matrix"); - chol_mat = R; -} + chol_mat = R; + } #if ! defined (HAVE_QRUPDATE) -template <typename T> -void -chol<T>::update (const VT& u) -{ - warn_qrupdate_once (); + template <typename T> + void + chol<T>::update (const VT& u) + { + warn_qrupdate_once (); - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (), - true, false); -} + init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (), + true, false); + } -template <typename T> -bool -singular (const T& a) -{ - static typename T::element_type zero (0); - for (octave_idx_type i = 0; i < a.rows (); i++) - if (a(i,i) == zero) return true; - return false; -} + template <typename T> + bool + singular (const T& a) + { + static typename T::element_type zero (0); + for (octave_idx_type i = 0; i < a.rows (); i++) + if (a(i,i) == zero) return true; + return false; + } + + template <typename T> + octave_idx_type + chol<T>::downdate (const VT& u) + { + warn_qrupdate_once (); -template <typename T> -octave_idx_type -chol<T>::downdate (const VT& u) -{ - warn_qrupdate_once (); + octave_idx_type info = -1; - octave_idx_type info = -1; + octave_idx_type n = chol_mat.rows (); - octave_idx_type n = chol_mat.rows (); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (singular (chol_mat)) + info = 2; + else + { + info = init (chol_mat.hermitian () * chol_mat + - T (u) * T (u).hermitian (), true, false); + if (info) info = 1; + } - if (singular (chol_mat)) - info = 2; - else - { - info = init (chol_mat.hermitian () * chol_mat - - T (u) * T (u).hermitian (), true, false); - if (info) info = 1; + return info; } - return info; -} + template <typename T> + octave_idx_type + chol<T>::insert_sym (const VT& u, octave_idx_type j) + { + static typename T::element_type zero (0); -template <typename T> -octave_idx_type -chol<T>::insert_sym (const VT& u, octave_idx_type j) -{ - static typename T::element_type zero (0); + warn_qrupdate_once (); - warn_qrupdate_once (); + octave_idx_type info = -1; - octave_idx_type info = -1; + octave_idx_type n = chol_mat.rows (); - octave_idx_type n = chol_mat.rows (); - - if (u.numel () != n + 1) - (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); - if (j < 0 || j > n) - (*current_liboctave_error_handler) ("cholinsert: index out of range"); + if (u.numel () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); - if (singular (chol_mat)) - info = 2; - else if (octave::math::imag (u(j)) != zero) - info = 3; - else - { - T a = chol_mat.hermitian () * chol_mat; - T a1 (n+1, n+1); - for (octave_idx_type k = 0; k < n+1; k++) - for (octave_idx_type l = 0; l < n+1; l++) - { - if (l == j) - a1(k, l) = u(k); - else if (k == j) - a1(k, l) = octave::math::conj (u(l)); - else - a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); - } - info = init (a1, true, false); - if (info) info = 1; + if (singular (chol_mat)) + info = 2; + else if (octave::math::imag (u(j)) != zero) + info = 3; + else + { + T a = chol_mat.hermitian () * chol_mat; + T a1 (n+1, n+1); + for (octave_idx_type k = 0; k < n+1; k++) + for (octave_idx_type l = 0; l < n+1; l++) + { + if (l == j) + a1(k, l) = u(k); + else if (k == j) + a1(k, l) = octave::math::conj (u(l)); + else + a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); + } + info = init (a1, true, false); + if (info) info = 1; + } + + return info; } - return info; -} - -template <typename T> -void -chol<T>::delete_sym (octave_idx_type j) -{ - warn_qrupdate_once (); + template <typename T> + void + chol<T>::delete_sym (octave_idx_type j) + { + warn_qrupdate_once (); - octave_idx_type n = chol_mat.rows (); - - if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("choldelete: index out of range"); - - T a = chol_mat.hermitian () * chol_mat; - a.delete_elements (1, idx_vector (j)); - a.delete_elements (0, idx_vector (j)); - init (a, true, false); -} + octave_idx_type n = chol_mat.rows (); -template <typename T> -void -chol<T>::shift_sym (octave_idx_type i, octave_idx_type j) -{ - warn_qrupdate_once (); - - octave_idx_type n = chol_mat.rows (); - - if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("cholshift: index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); - T a = chol_mat.hermitian () * chol_mat; - Array<octave_idx_type> p (dim_vector (n, 1)); - for (octave_idx_type k = 0; k < n; k++) p(k) = k; - if (i < j) - { - for (octave_idx_type k = i; k < j; k++) p(k) = k+1; - p(j) = i; - } - else if (j < i) - { - p(j) = i; - for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + T a = chol_mat.hermitian () * chol_mat; + a.delete_elements (1, idx_vector (j)); + a.delete_elements (0, idx_vector (j)); + init (a, true, false); } - init (a.index (idx_vector (p), idx_vector (p)), true, false); -} + template <typename T> + void + chol<T>::shift_sym (octave_idx_type i, octave_idx_type j) + { + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); + + T a = chol_mat.hermitian () * chol_mat; + Array<octave_idx_type> p (dim_vector (n, 1)); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + init (a.index (idx_vector (p), idx_vector (p)), true, false); + } #endif -// Specializations. + // Specializations. -template <> -octave_idx_type -chol<Matrix>::init (const Matrix& a, bool upper, bool calc_cond) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); + template <> + octave_idx_type + chol<Matrix>::init (const Matrix& a, bool upper, bool calc_cond) + { + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("chol: requires square matrix"); + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("chol: requires square matrix"); - octave_idx_type n = a_nc; - octave_idx_type info; + octave_idx_type n = a_nc; + octave_idx_type info; - is_upper = upper; + is_upper = upper; - chol_mat.clear (n, n); - if (is_upper) - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i <= j; i++) - chol_mat.xelem (i, j) = a(i, j); - for (octave_idx_type i = j+1; i < n; i++) - chol_mat.xelem (i, j) = 0.0; - } - else - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i < j; i++) - chol_mat.xelem (i, j) = 0.0; - for (octave_idx_type i = j; i < n; i++) - chol_mat.xelem (i, j) = a(i, j); - } - double *h = chol_mat.fortran_vec (); + chol_mat.clear (n, n); + if (is_upper) + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i <= j; i++) + chol_mat.xelem (i, j) = a(i, j); + for (octave_idx_type i = j+1; i < n; i++) + chol_mat.xelem (i, j) = 0.0; + } + else + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i < j; i++) + chol_mat.xelem (i, j) = 0.0; + for (octave_idx_type i = j; i < n; i++) + chol_mat.xelem (i, j) = a(i, j); + } + double *h = chol_mat.fortran_vec (); - // Calculate the norm of the matrix, for later use. - double anorm = 0; - if (calc_cond) - anorm = xnorm (a, 1); - - if (is_upper) - F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info - F77_CHAR_ARG_LEN (1))); - else - F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info - F77_CHAR_ARG_LEN (1))); + // Calculate the norm of the matrix, for later use. + double anorm = 0; + if (calc_cond) + anorm = xnorm (a, 1); - xrcond = 0.0; - if (info > 0) - chol_mat.resize (info - 1, info - 1); - else if (calc_cond) - { - octave_idx_type dpocon_info = 0; - - // Now calculate the condition number for non-singular matrix. - Array<double> z (dim_vector (3*n, 1)); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (dim_vector (n, 1)); - octave_idx_type *piz = iz.fortran_vec (); if (is_upper) - F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, piz, dpocon_info + F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info F77_CHAR_ARG_LEN (1))); else - F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, - n, anorm, xrcond, pz, piz, dpocon_info + F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info F77_CHAR_ARG_LEN (1))); - if (dpocon_info != 0) - info = -1; + xrcond = 0.0; + if (info > 0) + chol_mat.resize (info - 1, info - 1); + else if (calc_cond) + { + octave_idx_type dpocon_info = 0; + + // Now calculate the condition number for non-singular matrix. + Array<double> z (dim_vector (3*n, 1)); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (dim_vector (n, 1)); + octave_idx_type *piz = iz.fortran_vec (); + if (is_upper) + F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, + n, anorm, xrcond, pz, piz, dpocon_info + F77_CHAR_ARG_LEN (1))); + else + F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, + n, anorm, xrcond, pz, piz, dpocon_info + F77_CHAR_ARG_LEN (1))); + + if (dpocon_info != 0) + info = -1; + } + + return info; } - return info; -} - #if defined (HAVE_QRUPDATE) -template <> -void -chol<Matrix>::update (const ColumnVector& u) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<Matrix>::update (const ColumnVector& u) + { + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - ColumnVector utmp = u; + ColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), - utmp.fortran_vec (), w)); -} + F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); + } -template <> -octave_idx_type -chol<Matrix>::downdate (const ColumnVector& u) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<Matrix>::downdate (const ColumnVector& u) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - ColumnVector utmp = u; + ColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), - utmp.fortran_vec (), w, info)); + F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); - return info; -} + return info; + } -template <> -octave_idx_type -chol<Matrix>::insert_sym (const ColumnVector& u, octave_idx_type j) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<Matrix>::insert_sym (const ColumnVector& u, octave_idx_type j) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n + 1) - (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); - if (j < 0 || j > n) - (*current_liboctave_error_handler) ("cholinsert: index out of range"); + if (u.numel () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); - ColumnVector utmp = u; + ColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, w, n); - chol_mat.resize (n+1, n+1); + chol_mat.resize (n+1, n+1); - F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - j + 1, utmp.fortran_vec (), w, info)); + F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), w, info)); - return info; -} + return info; + } -template <> -void -chol<Matrix>::delete_sym (octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<Matrix>::delete_sym (octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("choldelete: index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - j + 1, w)); + F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, w)); - chol_mat.resize (n-1, n-1); -} + chol_mat.resize (n-1, n-1); + } -template <> -void -chol<Matrix>::shift_sym (octave_idx_type i, octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<Matrix>::shift_sym (octave_idx_type i, octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("cholshift: index out of range"); + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); - OCTAVE_LOCAL_BUFFER (double, w, 2*n); + OCTAVE_LOCAL_BUFFER (double, w, 2*n); - F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - i + 1, j + 1, w)); -} + F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } #endif -template <> -octave_idx_type -chol<FloatMatrix>::init (const FloatMatrix& a, bool upper, bool calc_cond) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); + template <> + octave_idx_type + chol<FloatMatrix>::init (const FloatMatrix& a, bool upper, bool calc_cond) + { + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("chol: requires square matrix"); + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("chol: requires square matrix"); - octave_idx_type n = a_nc; - octave_idx_type info; + octave_idx_type n = a_nc; + octave_idx_type info; - is_upper = upper; + is_upper = upper; - chol_mat.clear (n, n); - if (is_upper) - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i <= j; i++) - chol_mat.xelem (i, j) = a(i, j); - for (octave_idx_type i = j+1; i < n; i++) - chol_mat.xelem (i, j) = 0.0f; - } - else - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i < j; i++) - chol_mat.xelem (i, j) = 0.0f; - for (octave_idx_type i = j; i < n; i++) - chol_mat.xelem (i, j) = a(i, j); - } - float *h = chol_mat.fortran_vec (); + chol_mat.clear (n, n); + if (is_upper) + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i <= j; i++) + chol_mat.xelem (i, j) = a(i, j); + for (octave_idx_type i = j+1; i < n; i++) + chol_mat.xelem (i, j) = 0.0f; + } + else + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i < j; i++) + chol_mat.xelem (i, j) = 0.0f; + for (octave_idx_type i = j; i < n; i++) + chol_mat.xelem (i, j) = a(i, j); + } + float *h = chol_mat.fortran_vec (); - // Calculate the norm of the matrix, for later use. - float anorm = 0; - if (calc_cond) - anorm = xnorm (a, 1); - - if (is_upper) - F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info - F77_CHAR_ARG_LEN (1))); - else - F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info - F77_CHAR_ARG_LEN (1))); + // Calculate the norm of the matrix, for later use. + float anorm = 0; + if (calc_cond) + anorm = xnorm (a, 1); - xrcond = 0.0; - if (info > 0) - chol_mat.resize (info - 1, info - 1); - else if (calc_cond) - { - octave_idx_type spocon_info = 0; - - // Now calculate the condition number for non-singular matrix. - Array<float> z (dim_vector (3*n, 1)); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (dim_vector (n, 1)); - octave_idx_type *piz = iz.fortran_vec (); if (is_upper) - F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, piz, spocon_info + F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info F77_CHAR_ARG_LEN (1))); else - F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, - n, anorm, xrcond, pz, piz, spocon_info + F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info F77_CHAR_ARG_LEN (1))); - if (spocon_info != 0) - info = -1; + xrcond = 0.0; + if (info > 0) + chol_mat.resize (info - 1, info - 1); + else if (calc_cond) + { + octave_idx_type spocon_info = 0; + + // Now calculate the condition number for non-singular matrix. + Array<float> z (dim_vector (3*n, 1)); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (dim_vector (n, 1)); + octave_idx_type *piz = iz.fortran_vec (); + if (is_upper) + F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, + n, anorm, xrcond, pz, piz, spocon_info + F77_CHAR_ARG_LEN (1))); + else + F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, + n, anorm, xrcond, pz, piz, spocon_info + F77_CHAR_ARG_LEN (1))); + + if (spocon_info != 0) + info = -1; + } + + return info; } - return info; -} - #if defined (HAVE_QRUPDATE) -template <> -void -chol<FloatMatrix>::update (const FloatColumnVector& u) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatMatrix>::update (const FloatColumnVector& u) + { + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - FloatColumnVector utmp = u; + FloatColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), - utmp.fortran_vec (), w)); -} + F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); + } -template <> -octave_idx_type -chol<FloatMatrix>::downdate (const FloatColumnVector& u) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<FloatMatrix>::downdate (const FloatColumnVector& u) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - FloatColumnVector utmp = u; + FloatColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), - utmp.fortran_vec (), w, info)); + F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); - return info; -} + return info; + } -template <> -octave_idx_type -chol<FloatMatrix>::insert_sym (const FloatColumnVector& u, octave_idx_type j) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<FloatMatrix>::insert_sym (const FloatColumnVector& u, octave_idx_type j) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n + 1) - (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); - if (j < 0 || j > n) - (*current_liboctave_error_handler) ("cholinsert: index out of range"); + if (u.numel () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); - FloatColumnVector utmp = u; + FloatColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, w, n); - chol_mat.resize (n+1, n+1); + chol_mat.resize (n+1, n+1); - F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - j + 1, utmp.fortran_vec (), w, info)); + F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), w, info)); - return info; -} + return info; + } -template <> -void -chol<FloatMatrix>::delete_sym (octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatMatrix>::delete_sym (octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("choldelete: index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - j + 1, w)); + F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, w)); - chol_mat.resize (n-1, n-1); -} + chol_mat.resize (n-1, n-1); + } -template <> -void -chol<FloatMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("cholshift: index out of range"); + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); - OCTAVE_LOCAL_BUFFER (float, w, 2*n); + OCTAVE_LOCAL_BUFFER (float, w, 2*n); - F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), - i + 1, j + 1, w)); -} + F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } #endif -template <> -octave_idx_type -chol<ComplexMatrix>::init (const ComplexMatrix& a, bool upper, bool calc_cond) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); + template <> + octave_idx_type + chol<ComplexMatrix>::init (const ComplexMatrix& a, bool upper, bool calc_cond) + { + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("chol: requires square matrix"); + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("chol: requires square matrix"); - octave_idx_type n = a_nc; - octave_idx_type info; + octave_idx_type n = a_nc; + octave_idx_type info; - is_upper = upper; + is_upper = upper; - chol_mat.clear (n, n); - if (is_upper) - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i <= j; i++) - chol_mat.xelem (i, j) = a(i, j); - for (octave_idx_type i = j+1; i < n; i++) - chol_mat.xelem (i, j) = 0.0; - } - else - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i < j; i++) - chol_mat.xelem (i, j) = 0.0; - for (octave_idx_type i = j; i < n; i++) - chol_mat.xelem (i, j) = a(i, j); - } - Complex *h = chol_mat.fortran_vec (); + chol_mat.clear (n, n); + if (is_upper) + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i <= j; i++) + chol_mat.xelem (i, j) = a(i, j); + for (octave_idx_type i = j+1; i < n; i++) + chol_mat.xelem (i, j) = 0.0; + } + else + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i < j; i++) + chol_mat.xelem (i, j) = 0.0; + for (octave_idx_type i = j; i < n; i++) + chol_mat.xelem (i, j) = a(i, j); + } + Complex *h = chol_mat.fortran_vec (); - // Calculate the norm of the matrix, for later use. - double anorm = 0; - if (calc_cond) - anorm = xnorm (a, 1); + // Calculate the norm of the matrix, for later use. + double anorm = 0; + if (calc_cond) + anorm = xnorm (a, 1); - if (is_upper) - F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (h), n, info - F77_CHAR_ARG_LEN (1))); - else - F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, F77_DBLE_CMPLX_ARG (h), n, info - F77_CHAR_ARG_LEN (1))); + if (is_upper) + F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (h), n, info + F77_CHAR_ARG_LEN (1))); + else + F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, F77_DBLE_CMPLX_ARG (h), n, info + F77_CHAR_ARG_LEN (1))); - xrcond = 0.0; - if (info > 0) - chol_mat.resize (info - 1, info - 1); - else if (calc_cond) - { - octave_idx_type zpocon_info = 0; + xrcond = 0.0; + if (info > 0) + chol_mat.resize (info - 1, info - 1); + else if (calc_cond) + { + octave_idx_type zpocon_info = 0; - // Now calculate the condition number for non-singular matrix. - Array<Complex> z (dim_vector (2*n, 1)); - Complex *pz = z.fortran_vec (); - Array<double> rz (dim_vector (n, 1)); - double *prz = rz.fortran_vec (); - F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (h), - n, anorm, xrcond, F77_DBLE_CMPLX_ARG (pz), prz, zpocon_info - F77_CHAR_ARG_LEN (1))); + // Now calculate the condition number for non-singular matrix. + Array<Complex> z (dim_vector (2*n, 1)); + Complex *pz = z.fortran_vec (); + Array<double> rz (dim_vector (n, 1)); + double *prz = rz.fortran_vec (); + F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (h), + n, anorm, xrcond, F77_DBLE_CMPLX_ARG (pz), prz, zpocon_info + F77_CHAR_ARG_LEN (1))); - if (zpocon_info != 0) - info = -1; + if (zpocon_info != 0) + info = -1; + } + + return info; } - return info; -} - #if defined (HAVE_QRUPDATE) -template <> -void -chol<ComplexMatrix>::update (const ComplexColumnVector& u) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<ComplexMatrix>::update (const ComplexColumnVector& u) + { + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - ComplexColumnVector utmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, rw, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1up, ZCH1UP, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw)); -} + F77_XFCN (zch1up, ZCH1UP, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw)); + } -template <> -octave_idx_type -chol<ComplexMatrix>::downdate (const ComplexColumnVector& u) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<ComplexMatrix>::downdate (const ComplexColumnVector& u) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - ComplexColumnVector utmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, rw, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1dn, ZCH1DN, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); + F77_XFCN (zch1dn, ZCH1DN, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); - return info; -} + return info; + } -template <> -octave_idx_type -chol<ComplexMatrix>::insert_sym (const ComplexColumnVector& u, - octave_idx_type j) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<ComplexMatrix>::insert_sym (const ComplexColumnVector& u, + octave_idx_type j) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n + 1) - (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); - if (j < 0 || j > n) - (*current_liboctave_error_handler) ("cholinsert: index out of range"); + if (u.numel () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); - ComplexColumnVector utmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, rw, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - chol_mat.resize (n+1, n+1); + chol_mat.resize (n+1, n+1); - F77_XFCN (zchinx, ZCHINX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - j + 1, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); + F77_XFCN (zchinx, ZCHINX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + j + 1, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); - return info; -} + return info; + } -template <> -void -chol<ComplexMatrix>::delete_sym (octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<ComplexMatrix>::delete_sym (octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("choldelete: index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); - OCTAVE_LOCAL_BUFFER (double, rw, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchdex, ZCHDEX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - j + 1, rw)); + F77_XFCN (zchdex, ZCHDEX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + j + 1, rw)); - chol_mat.resize (n-1, n-1); -} + chol_mat.resize (n-1, n-1); + } -template <> -void -chol<ComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<ComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("cholshift: index out of range"); + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); - OCTAVE_LOCAL_BUFFER (Complex, w, n); - OCTAVE_LOCAL_BUFFER (double, rw, n); + OCTAVE_LOCAL_BUFFER (Complex, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchshx, ZCHSHX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw)); -} + F77_XFCN (zchshx, ZCHSHX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw)); + } #endif -template <> -octave_idx_type -chol<FloatComplexMatrix>::init (const FloatComplexMatrix& a, bool upper, - bool calc_cond) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); + template <> + octave_idx_type + chol<FloatComplexMatrix>::init (const FloatComplexMatrix& a, bool upper, + bool calc_cond) + { + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("chol: requires square matrix"); + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("chol: requires square matrix"); - octave_idx_type n = a_nc; - octave_idx_type info; + octave_idx_type n = a_nc; + octave_idx_type info; - is_upper = upper; + is_upper = upper; - chol_mat.clear (n, n); - if (is_upper) - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i <= j; i++) - chol_mat.xelem (i, j) = a(i, j); - for (octave_idx_type i = j+1; i < n; i++) - chol_mat.xelem (i, j) = 0.0f; - } - else - for (octave_idx_type j = 0; j < n; j++) - { - for (octave_idx_type i = 0; i < j; i++) - chol_mat.xelem (i, j) = 0.0f; - for (octave_idx_type i = j; i < n; i++) - chol_mat.xelem (i, j) = a(i, j); - } - FloatComplex *h = chol_mat.fortran_vec (); + chol_mat.clear (n, n); + if (is_upper) + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i <= j; i++) + chol_mat.xelem (i, j) = a(i, j); + for (octave_idx_type i = j+1; i < n; i++) + chol_mat.xelem (i, j) = 0.0f; + } + else + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = 0; i < j; i++) + chol_mat.xelem (i, j) = 0.0f; + for (octave_idx_type i = j; i < n; i++) + chol_mat.xelem (i, j) = a(i, j); + } + FloatComplex *h = chol_mat.fortran_vec (); - // Calculate the norm of the matrix, for later use. - float anorm = 0; - if (calc_cond) - anorm = xnorm (a, 1); + // Calculate the norm of the matrix, for later use. + float anorm = 0; + if (calc_cond) + anorm = xnorm (a, 1); - if (is_upper) - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h), n, info - F77_CHAR_ARG_LEN (1))); - else - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, F77_CMPLX_ARG (h), n, info - F77_CHAR_ARG_LEN (1))); + if (is_upper) + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h), n, info + F77_CHAR_ARG_LEN (1))); + else + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, F77_CMPLX_ARG (h), n, info + F77_CHAR_ARG_LEN (1))); - xrcond = 0.0; - if (info > 0) - chol_mat.resize (info - 1, info - 1); - else if (calc_cond) - { - octave_idx_type cpocon_info = 0; + xrcond = 0.0; + if (info > 0) + chol_mat.resize (info - 1, info - 1); + else if (calc_cond) + { + octave_idx_type cpocon_info = 0; - // Now calculate the condition number for non-singular matrix. - Array<FloatComplex> z (dim_vector (2*n, 1)); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (dim_vector (n, 1)); - float *prz = rz.fortran_vec (); - F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h), - n, anorm, xrcond, F77_CMPLX_ARG (pz), prz, cpocon_info - F77_CHAR_ARG_LEN (1))); + // Now calculate the condition number for non-singular matrix. + Array<FloatComplex> z (dim_vector (2*n, 1)); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (dim_vector (n, 1)); + float *prz = rz.fortran_vec (); + F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h), + n, anorm, xrcond, F77_CMPLX_ARG (pz), prz, cpocon_info + F77_CHAR_ARG_LEN (1))); - if (cpocon_info != 0) - info = -1; + if (cpocon_info != 0) + info = -1; + } + + return info; } - return info; -} - #if defined (HAVE_QRUPDATE) -template <> -void -chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u) + { + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - FloatComplexColumnVector utmp = u; + FloatComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, rw, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cch1up, CCH1UP, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - F77_CMPLX_ARG (utmp.fortran_vec ()), rw)); -} + F77_XFCN (cch1up, CCH1UP, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + F77_CMPLX_ARG (utmp.fortran_vec ()), rw)); + } -template <> -octave_idx_type -chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n) - (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + if (u.numel () != n) + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); - FloatComplexColumnVector utmp = u; + FloatComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, rw, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cch1dn, CCH1DN, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); + F77_XFCN (cch1dn, CCH1DN, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); - return info; -} + return info; + } -template <> -octave_idx_type -chol<FloatComplexMatrix>::insert_sym (const FloatComplexColumnVector& u, - octave_idx_type j) -{ - octave_idx_type info = -1; + template <> + octave_idx_type + chol<FloatComplexMatrix>::insert_sym (const FloatComplexColumnVector& u, + octave_idx_type j) + { + octave_idx_type info = -1; - octave_idx_type n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); - if (u.numel () != n + 1) - (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); - if (j < 0 || j > n) - (*current_liboctave_error_handler) ("cholinsert: index out of range"); + if (u.numel () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); - FloatComplexColumnVector utmp = u; + FloatComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, rw, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - chol_mat.resize (n+1, n+1); + chol_mat.resize (n+1, n+1); - F77_XFCN (cchinx, CCHINX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - j + 1, F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); + F77_XFCN (cchinx, CCHINX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + j + 1, F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info)); - return info; -} + return info; + } -template <> -void -chol<FloatComplexMatrix>::delete_sym (octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatComplexMatrix>::delete_sym (octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("choldelete: index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); - OCTAVE_LOCAL_BUFFER (float, rw, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cchdex, CCHDEX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - j + 1, rw)); + F77_XFCN (cchdex, CCHDEX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + j + 1, rw)); - chol_mat.resize (n-1, n-1); -} + chol_mat.resize (n-1, n-1); + } -template <> -void -chol<FloatComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) -{ - octave_idx_type n = chol_mat.rows (); + template <> + void + chol<FloatComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j) + { + octave_idx_type n = chol_mat.rows (); - if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("cholshift: index out of range"); + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); - OCTAVE_LOCAL_BUFFER (FloatComplex, w, n); - OCTAVE_LOCAL_BUFFER (float, rw, n); + OCTAVE_LOCAL_BUFFER (FloatComplex, w, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cchshx, CCHSHX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), - i + 1, j + 1, F77_CMPLX_ARG (w), rw)); -} + F77_XFCN (cchshx, CCHSHX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()), chol_mat.rows (), + i + 1, j + 1, F77_CMPLX_ARG (w), rw)); + } #endif -// Instantiations we need. + // Instantiations we need. -template class chol<Matrix>; + template class chol<Matrix>; -template class chol<FloatMatrix>; + template class chol<FloatMatrix>; -template class chol<ComplexMatrix>; + template class chol<ComplexMatrix>; -template class chol<FloatComplexMatrix>; + template class chol<FloatComplexMatrix>; -template Matrix -chol2inv<Matrix> (const Matrix& r); + template Matrix + chol2inv<Matrix> (const Matrix& r); -template ComplexMatrix -chol2inv<ComplexMatrix> (const ComplexMatrix& r); + template ComplexMatrix + chol2inv<ComplexMatrix> (const ComplexMatrix& r); -template FloatMatrix -chol2inv<FloatMatrix> (const FloatMatrix& r); + template FloatMatrix + chol2inv<FloatMatrix> (const FloatMatrix& r); -template FloatComplexMatrix -chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r); - + template FloatComplexMatrix + chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r); + } } -}