Mercurial > hg > octave-jordi
diff liboctave/floatCHOL.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/floatCHOL.cc +++ b/liboctave/floatCHOL.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 @@ float*, const octave_idx_type&, const float&, float&, float*, octave_idx_type*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*, float*, float*); + F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*, const octave_idx_type&, + float*, float*); F77_RET_T - F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*, float*, float*, + F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*, const octave_idx_type&, + float*, float*, octave_idx_type&); + + F77_RET_T + F77_FUNC (schinx, SCHINX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*, float*, octave_idx_type&); F77_RET_T - F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, float*, const octave_idx_type&, const octave_idx_type&); + F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T - F77_FUNC (schinx, SCHINX) (const octave_idx_type&, const float*, float*, const octave_idx_type&, - const float*, octave_idx_type&); - - F77_RET_T - F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, const float*, float*, const octave_idx_type&); + F77_FUNC (schshx, SCHSHX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + float*); +#endif } octave_idx_type @@ -186,26 +192,28 @@ (*current_liboctave_error_handler) ("FloatCHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -FloatCHOL::update (const FloatMatrix& u) +FloatCHOL::update (const FloatColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - FloatMatrix tmp = u; + FloatColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("FloatCHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -FloatCHOL::downdate (const FloatMatrix& u) +FloatCHOL::downdate (const FloatColumnVector& u) { octave_idx_type info = -1; @@ -213,38 +221,40 @@ if (u.length () == n) { - FloatMatrix tmp = u; + FloatColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("FloatCHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -FloatCHOL::insert_sym (const FloatMatrix& u, octave_idx_type j) +FloatCHOL::insert_sym (const FloatColumnVector& 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) ("FloatCHOL 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) ("FloatCHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - FloatMatrix chol_mat1 (n+1, n+1); + FloatColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (schinx, SCHINX, (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 (schinx, SCHINX, (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) ("FloatCHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - FloatMatrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (schdex, SCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (schdex, SCHDEX, (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 @@ FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j) { octave_idx_type n = chol_mat.rows (); - float dummy; if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("FloatCHOL shift index out of range"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (sqrshc, SQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (float, w, 2*n); + + F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } } +#endif + FloatMatrix chol2inv (const FloatMatrix& r) {