Mercurial > hg > octave-lyh
diff liboctave/CmplxCHOL.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/CmplxCHOL.cc +++ b/liboctave/CmplxCHOL.cc @@ -21,10 +21,14 @@ */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include <config.h> #endif +#include <vector> + #include "dMatrix.h" #include "dRowVector.h" #include "CmplxCHOL.h" @@ -47,6 +51,12 @@ Complex*, const octave_idx_type&, const double&, double&, Complex*, double*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); + F77_RET_T + F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, Complex*, double*); + + F77_RET_T + F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*, + octave_idx_type&); } octave_idx_type @@ -151,6 +161,55 @@ return chol2inv_internal (chol_mat); } +void +ComplexCHOL::set (const ComplexMatrix& R) +{ + if (R.is_square ()) + chol_mat = R; + else + (*current_liboctave_error_handler) ("chol2inv requires square matrix"); +} + +void +ComplexCHOL::update (const ComplexMatrix& u) +{ + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + ComplexMatrix tmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); + + F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), + tmp.fortran_vec (), w)); + } + else + (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); +} + +octave_idx_type +ComplexCHOL::downdate (const ComplexMatrix& u) +{ + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + ComplexMatrix tmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); + + F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), + tmp.fortran_vec (), w, info)); + } + else + (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + + return info; +} + ComplexMatrix chol2inv (const ComplexMatrix& r) {