Mercurial > hg > octave-nkf
diff liboctave/dbleCHOL.cc @ 8547:d66c9b6e506a
imported patch qrupdate.diff
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 20 Jan 2009 21:16:42 +0100 |
parents | 25bc2d31e1bf |
children | a6edd5c23cb5 |
line wrap: on
line diff
--- a/liboctave/dbleCHOL.cc +++ b/liboctave/dbleCHOL.cc @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -52,23 +51,30 @@ double*, const octave_idx_type&, const double&, double&, double*, octave_idx_type*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*); + F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, const octave_idx_type&, + double*, double*); F77_RET_T - F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, + F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, const octave_idx_type&, + double*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, double*, double*, octave_idx_type&); F77_RET_T - F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - double*, double*, const octave_idx_type&, const octave_idx_type&); + F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, double*); F77_RET_T - F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, const double*, double*, const octave_idx_type&, - const double*, octave_idx_type&); - - F77_RET_T - F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, const double*, double*, const octave_idx_type&); + F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + double*); +#endif } octave_idx_type @@ -186,26 +192,28 @@ (*current_liboctave_error_handler) ("CHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -CHOL::update (const Matrix& u) +CHOL::update (const ColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - Matrix tmp = u; + ColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -CHOL::downdate (const Matrix& u) +CHOL::downdate (const ColumnVector& u) { octave_idx_type info = -1; @@ -213,38 +221,40 @@ if (u.length () == n) { - Matrix tmp = u; + ColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -CHOL::insert_sym (const Matrix& u, octave_idx_type j) +CHOL::insert_sym (const ColumnVector& u, octave_idx_type j) { octave_idx_type info = -1; octave_idx_type n = chol_mat.rows (); - if (u.length () != n+1) - (*current_liboctave_error_handler) ("CHOL insert dimension mismatch"); + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("CHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - Matrix chol_mat1 (n+1, n+1); + ColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), - j+1, u.data (), info)); + chol_mat.resize (n+1, n+1); - chol_mat = chol_mat1; + F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), w, info)); } return info; @@ -256,14 +266,15 @@ octave_idx_type n = chol_mat.rows (); if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - Matrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, w)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -271,14 +282,20 @@ CHOL::shift_sym (octave_idx_type i, octave_idx_type j) { octave_idx_type n = chol_mat.rows (); - double dummy; if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL shift index out of range"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (double, w, 2*n); + + F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } } +#endif + Matrix chol2inv (const Matrix& r) {