Mercurial > hg > octave-nkf
diff liboctave/dbleCHOL.cc @ 7554:40574114c514
implement Cholesky factorization updating
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 04 Mar 2008 22:25:50 -0500 |
parents | 29980c6b8604 |
children | 07522d7dcdf8 |
line wrap: on
line diff
--- a/liboctave/dbleCHOL.cc +++ b/liboctave/dbleCHOL.cc @@ -21,10 +21,14 @@ */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include <config.h> #endif +#include <vector> + #include "dRowVector.h" #include "dbleCHOL.h" #include "f77-fcn.h" @@ -47,6 +51,12 @@ double*, const octave_idx_type&, const double&, double&, double*, octave_idx_type*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); + F77_RET_T + F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*); + + F77_RET_T + F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, + int &); } octave_idx_type @@ -155,6 +165,55 @@ return chol2inv_internal (chol_mat); } +void +CHOL::set (const Matrix& R) +{ + if (R.is_square ()) + chol_mat = R; + else + (*current_liboctave_error_handler) ("chol2inv requires square matrix"); +} + +void +CHOL::update (const Matrix& u) +{ + int n = chol_mat.rows (); + + if (u.length () == n) + { + Matrix tmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); + + F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), + tmp.fortran_vec (), w)); + } + else + (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); +} + +octave_idx_type +CHOL::downdate (const Matrix& u) +{ + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + Matrix tmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); + + F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), + tmp.fortran_vec (), w, info)); + } + else + (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + + return info; +} + Matrix chol2inv (const Matrix& r) {