Mercurial > hg > octave-lyh
diff liboctave/CmplxCHOL.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/CmplxCHOL.cc +++ b/liboctave/CmplxCHOL.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 @@ Complex*, const octave_idx_type&, const double&, double&, Complex*, double*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, Complex*, double*); + F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, double*); F77_RET_T - F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*, + F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, Complex*, double*, octave_idx_type&); F77_RET_T - F77_FUNC (zqrshc, ZQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - Complex*, Complex*, const octave_idx_type&, const octave_idx_type&); + F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, double*); F77_RET_T - F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, const Complex*, Complex*, const octave_idx_type&, - const Complex*, octave_idx_type&); - - F77_RET_T - F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, const Complex*, Complex*, const octave_idx_type&); + F77_FUNC (zchshx, ZCHSHX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + Complex*, double*); +#endif } octave_idx_type @@ -182,26 +188,28 @@ (*current_liboctave_error_handler) ("CHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -ComplexCHOL::update (const ComplexMatrix& u) +ComplexCHOL::update (const ComplexColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - ComplexMatrix tmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -ComplexCHOL::downdate (const ComplexMatrix& u) +ComplexCHOL::downdate (const ComplexColumnVector& u) { octave_idx_type info = -1; @@ -209,38 +217,40 @@ if (u.length () == n) { - ComplexMatrix tmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw, info)); } else - (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -ComplexCHOL::insert_sym (const ComplexMatrix& u, octave_idx_type j) +ComplexCHOL::insert_sym (const ComplexColumnVector& 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 { - ComplexMatrix chol_mat1 (n+1, n+1); + ComplexColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchinx, ZCHINX, (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 (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), rw, info)); } return info; @@ -252,14 +262,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 { - ComplexMatrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, rw)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -267,14 +278,21 @@ ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) { octave_idx_type n = chol_mat.rows (); - Complex 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 (zqrshc, ZQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (Complex, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); + + F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w, rw)); + } } +#endif + ComplexMatrix chol2inv (const ComplexMatrix& r) {