# HG changeset patch # User Jaroslav Hajek # Date 1232619047 -3600 # Node ID a6edd5c23cb5e9c70c4ce78552d67ba0ccd3cfe2 # Parent 66165de2cc4285857c07812f481b3d6d1c2dbea3 use replacement methods if qrupdate is not available diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2009-01-22 Jaroslav Hajek + + * configure.in: Fix qrupdate warning message. + 2009-01-21 John W. Eaton * Makeconf.in: Substitute X11_INCFLAGS and X11_LIBS. diff --git a/configure.in b/configure.in --- a/configure.in +++ b/configure.in @@ -885,7 +885,7 @@ [don't use qrupdate, disable QR & Cholesky updating functions])], with_qrupdate=$withval, with_qrupdate=yes) -warn_qrupdate="qrupdate not found. The QR & Cholesky updating function will not be available." +warn_qrupdate="qrupdate not found. The QR & Cholesky updating functions will be slow." if test "$with_qrupdate" = yes; then with_qrupdate=no if $have_fortran_compiler; then diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,16 @@ +2009-01-22 Jaroslav Hajek + + * dbleQR.h: Optionally declare warn_qrupdate_once. + * dbleQR.cc: Define it. + * (CmplxQR.h, dbleQR.h, fCmplxQR.h, floatQR.h): Declare replacement + methods unconditionally. + * (CmplxQR.cc, dbleQR.cc, fCmplxQR.cc, floatQR.cc): Define + updating replacement methods. + * (CmplxCHOL.h, dbleCHOL.h, fCmplxCHOL.h, floatCHOL.h): Declare + replacement methods unconditionally. + * (CmplxCHOL.cc, dbleCHOL.cc, fCmplxCHOL.cc, floatCHOL.cc): Define + updating replacement methods. + 2009-01-21 Jaroslav Hajek * Range.cc ( operator + (double x, const Range& r), diff --git a/liboctave/CmplxCHOL.cc b/liboctave/CmplxCHOL.cc --- a/liboctave/CmplxCHOL.cc +++ b/liboctave/CmplxCHOL.cc @@ -34,6 +34,9 @@ #include "f77-fcn.h" #include "lo-error.h" #include "oct-locbuf.h" +#ifndef HAVE_QRUPDATE +#include "dbleQR.h" +#endif extern "C" { @@ -291,6 +294,146 @@ } } +#else + +void +ComplexCHOL::update (const ComplexColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + init (chol_mat.hermitian () * chol_mat + + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false); + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); +} + +static bool +singular (const ComplexMatrix& a) +{ + for (octave_idx_type i = 0; i < a.rows (); i++) + if (a(i,i) == 0.0) return true; + return false; +} + +octave_idx_type +ComplexCHOL::downdate (const ComplexColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + if (singular (chol_mat)) + info = 2; + else + { + info = init (chol_mat.hermitian () * chol_mat + - ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false); + if (info) info = 1; + } + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + + return info; +} + +octave_idx_type +ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); + else + { + if (singular (chol_mat)) + info = 2; + else if (u(j).imag () != 0.0) + info = 3; + else + { + ComplexMatrix a = chol_mat.hermitian () * chol_mat; + ComplexMatrix a1 (n+1, n+1); + for (octave_idx_type k = 0; k < n+1; k++) + for (octave_idx_type l = 0; l < n+1; l++) + { + if (l == j) + a1(k, l) = u(k); + else if (k == j) + a1(k, l) = std::conj (u(l)); + else + a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); + } + info = init (a1, false); + if (info) info = 1; + } + } + + return info; +} + +void +ComplexCHOL::delete_sym (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); + else + { + ComplexMatrix a = chol_mat.hermitian () * chol_mat; + a.delete_elements (1, idx_vector (j)); + a.delete_elements (0, idx_vector (j)); + init (a, false); + } +} + +void +ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); + else + { + ComplexMatrix a = chol_mat.hermitian () * chol_mat; + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + init (a.index (idx_vector (p), idx_vector (p)), false); + } +} + #endif ComplexMatrix diff --git a/liboctave/CmplxCHOL.h b/liboctave/CmplxCHOL.h --- a/liboctave/CmplxCHOL.h +++ b/liboctave/CmplxCHOL.h @@ -67,8 +67,6 @@ void set (const ComplexMatrix& R); -#ifdef HAVE_QRUPDATE - void update (const ComplexColumnVector& u); octave_idx_type downdate (const ComplexColumnVector& u); @@ -79,8 +77,6 @@ void shift_sym (octave_idx_type i, octave_idx_type j); -#endif - friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a); private: diff --git a/liboctave/CmplxQR.cc b/liboctave/CmplxQR.cc --- a/liboctave/CmplxQR.cc +++ b/liboctave/CmplxQR.cc @@ -168,14 +168,28 @@ ComplexQR::ComplexQR (const ComplexMatrix& q_arg, const ComplexMatrix& r_arg) { - if (q_arg.columns () != r_arg.rows ()) + octave_idx_type qr = q_arg.rows (), qc = q_arg.columns (); + octave_idx_type rr = r_arg.rows (), rc = r_arg.columns (); + if (qc == rr && (qr == qc || (qr > qc && rr == rc))) { - (*current_liboctave_error_handler) ("QR dimensions mismatch"); - return; + q = q_arg; + r = r_arg; } + else + (*current_liboctave_error_handler) ("QR dimensions mismatch"); +} - this->q = q_arg; - this->r = r_arg; +QR::type +ComplexQR::get_type (void) const +{ + QR::type retval; + if (!q.is_empty () && q.is_square ()) + retval = QR::std; + else if (q.rows () > q.columns () && r.is_square ()) + retval = QR::economy; + else + retval = QR::raw; + return retval; } #ifdef HAVE_QRUPDATE @@ -196,7 +210,7 @@ utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void @@ -431,6 +445,249 @@ } } +#else + +// Replacement update methods. + +void +ComplexQR::update (const ComplexColumnVector& u, const ComplexColumnVector& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () == m && v.length () == n) + { + init(q*r + ComplexMatrix (u) * ComplexMatrix (v).hermitian (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +void +ComplexQR::update (const ComplexMatrix& u, const ComplexMatrix& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + init(q*r + u * v.hermitian (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +static +ComplexMatrix insert_col (const ComplexMatrix& a, octave_idx_type i, + const ComplexColumnVector& x) +{ + ComplexMatrix retval (a.rows (), a.columns () + 1); + retval.assign (idx_vector::colon, idx_vector (0, i), + a.index (idx_vector::colon, idx_vector (0, i))); + retval.assign (idx_vector::colon, idx_vector (i), x); + retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()), + a.index (idx_vector::colon, idx_vector (i, a.columns ()))); + return retval; +} + +static +ComplexMatrix insert_row (const ComplexMatrix& a, octave_idx_type i, + const ComplexRowVector& x) +{ + ComplexMatrix retval (a.rows () + 1, a.columns ()); + retval.assign (idx_vector (0, i), idx_vector::colon, + a.index (idx_vector (0, i), idx_vector::colon)); + retval.assign (idx_vector (i), idx_vector::colon, x); + retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon, + a.index (idx_vector (i, a.rows ()), idx_vector::colon)); + return retval; +} + +static +ComplexMatrix delete_col (const ComplexMatrix& a, octave_idx_type i) +{ + ComplexMatrix retval = a; + retval.delete_elements (1, idx_vector (i)); + return retval; +} + +static +ComplexMatrix delete_row (const ComplexMatrix& a, octave_idx_type i) +{ + ComplexMatrix retval = a; + retval.delete_elements (0, idx_vector (i)); + return retval; +} + +static +ComplexMatrix shift_cols (const ComplexMatrix& a, + octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = a.columns (); + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + return a.index (idx_vector::colon, idx_vector (p)); +} + +void +ComplexQR::insert_col (const ComplexColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_col (q*r, j, u), get_type ()); + } +} + +void +ComplexQR::insert_col (const ComplexMatrix& u, const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + ComplexMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::insert_col (a, js(i), u.column (i)); + init (a, get_type ()); + } +} + +void +ComplexQR::delete_col (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_col (q*r, j), get_type ()); + } +} + +void +ComplexQR::delete_col (const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + ComplexMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::delete_col (a, js(i)); + init (a, get_type ()); + } +} + +void +ComplexQR::insert_row (const ComplexRowVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_row (q*r, j, u), get_type ()); + } +} + +void +ComplexQR::delete_row (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_row (q*r, j), get_type ()); + } +} + +void +ComplexQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrshift: index out of range"); + else + { + init (::shift_cols (q*r, i, j), get_type ()); + } +} + #endif /* diff --git a/liboctave/CmplxQR.h b/liboctave/CmplxQR.h --- a/liboctave/CmplxQR.h +++ b/liboctave/CmplxQR.h @@ -65,7 +65,7 @@ ComplexMatrix R (void) const { return r; } -#ifdef HAVE_QRUPDATE + QR::type get_type (void) const; void update (const ComplexColumnVector& u, const ComplexColumnVector& v); @@ -85,8 +85,6 @@ void shift_cols (octave_idx_type i, octave_idx_type j); -#endif - friend std::ostream& operator << (std::ostream&, const ComplexQR&); protected: diff --git a/liboctave/dbleCHOL.cc b/liboctave/dbleCHOL.cc --- a/liboctave/dbleCHOL.cc +++ b/liboctave/dbleCHOL.cc @@ -33,6 +33,9 @@ #include "f77-fcn.h" #include "lo-error.h" #include "oct-locbuf.h" +#ifndef HAVE_QRUPDATE +#include "dbleQR.h" +#endif extern "C" { @@ -294,6 +297,144 @@ } } +#else + +void +CHOL::update (const ColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + init (chol_mat.transpose () * chol_mat + + Matrix (u) * Matrix (u).transpose (), false); + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); +} + +static bool +singular (const Matrix& a) +{ + for (octave_idx_type i = 0; i < a.rows (); i++) + if (a(i,i) == 0.0) return true; + return false; +} + +octave_idx_type +CHOL::downdate (const ColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + if (singular (chol_mat)) + info = 2; + else + { + info = init (chol_mat.transpose () * chol_mat + - Matrix (u) * Matrix (u).transpose (), false); + if (info) info = 1; + } + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + + return info; +} + +octave_idx_type +CHOL::insert_sym (const ColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); + else + { + if (singular (chol_mat)) + info = 2; + else + { + Matrix a = chol_mat.transpose () * chol_mat; + Matrix a1 (n+1, n+1); + for (octave_idx_type k = 0; k < n+1; k++) + for (octave_idx_type l = 0; l < n+1; l++) + { + if (l == j) + a1(k, l) = u(k); + else if (k == j) + a1(k, l) = u(l); + else + a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); + } + info = init (a1, false); + if (info) info = 1; + } + } + + return info; +} + +void +CHOL::delete_sym (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); + else + { + Matrix a = chol_mat.transpose () * chol_mat; + a.delete_elements (1, idx_vector (j)); + a.delete_elements (0, idx_vector (j)); + init (a, false); + } +} + +void +CHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); + else + { + Matrix a = chol_mat.transpose () * chol_mat; + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + init (a.index (idx_vector (p), idx_vector (p)), false); + } +} + #endif Matrix diff --git a/liboctave/dbleCHOL.h b/liboctave/dbleCHOL.h --- a/liboctave/dbleCHOL.h +++ b/liboctave/dbleCHOL.h @@ -64,8 +64,6 @@ void set (const Matrix& R); -#ifdef HAVE_QRUPDATE - void update (const ColumnVector& u); octave_idx_type downdate (const ColumnVector& u); @@ -76,8 +74,6 @@ void shift_sym (octave_idx_type i, octave_idx_type j); -#endif - friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a); private: diff --git a/liboctave/dbleQR.cc b/liboctave/dbleQR.cc --- a/liboctave/dbleQR.cc +++ b/liboctave/dbleQR.cc @@ -159,14 +159,28 @@ QR::QR (const Matrix& q_arg, const Matrix& r_arg) { - if (q_arg.columns () != r_arg.rows ()) + octave_idx_type qr = q_arg.rows (), qc = q_arg.columns (); + octave_idx_type rr = r_arg.rows (), rc = r_arg.columns (); + if (qc == rr && (qr == qc || (qr > qc && rr == rc))) { - (*current_liboctave_error_handler) ("QR dimensions mismatch"); - return; + q = q_arg; + r = r_arg; } + else + (*current_liboctave_error_handler) ("QR dimensions mismatch"); +} - this->q = q_arg; - this->r = r_arg; +QR::type +QR::get_type (void) const +{ + QR::type retval; + if (!q.is_empty () && q.is_square ()) + retval = QR::std; + else if (q.rows () > q.columns () && r.is_square ()) + retval = QR::economy; + else + retval = QR::raw; + return retval; } #ifdef HAVE_QRUPDATE @@ -186,7 +200,7 @@ utmp.fortran_vec (), vtmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void @@ -418,6 +432,261 @@ } } +#else + +// Replacement update methods. + +void +QR::update (const ColumnVector& u, const ColumnVector& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () == m && v.length () == n) + { + init(q*r + Matrix (u) * Matrix (v).transpose (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +void +QR::update (const Matrix& u, const Matrix& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + init(q*r + u * v.transpose (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +static +Matrix insert_col (const Matrix& a, octave_idx_type i, + const ColumnVector& x) +{ + Matrix retval (a.rows (), a.columns () + 1); + retval.assign (idx_vector::colon, idx_vector (0, i), + a.index (idx_vector::colon, idx_vector (0, i))); + retval.assign (idx_vector::colon, idx_vector (i), x); + retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()), + a.index (idx_vector::colon, idx_vector (i, a.columns ()))); + return retval; +} + +static +Matrix insert_row (const Matrix& a, octave_idx_type i, + const RowVector& x) +{ + Matrix retval (a.rows () + 1, a.columns ()); + retval.assign (idx_vector (0, i), idx_vector::colon, + a.index (idx_vector (0, i), idx_vector::colon)); + retval.assign (idx_vector (i), idx_vector::colon, x); + retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon, + a.index (idx_vector (i, a.rows ()), idx_vector::colon)); + return retval; +} + +static +Matrix delete_col (const Matrix& a, octave_idx_type i) +{ + Matrix retval = a; + retval.delete_elements (1, idx_vector (i)); + return retval; +} + +static +Matrix delete_row (const Matrix& a, octave_idx_type i) +{ + Matrix retval = a; + retval.delete_elements (0, idx_vector (i)); + return retval; +} + +static +Matrix shift_cols (const Matrix& a, + octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = a.columns (); + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + return a.index (idx_vector::colon, idx_vector (p)); +} + +void +QR::insert_col (const ColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_col (q*r, j, u), get_type ()); + } +} + +void +QR::insert_col (const Matrix& u, const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + Matrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::insert_col (a, js(i), u.column (i)); + init (a, get_type ()); + } +} + +void +QR::delete_col (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_col (q*r, j), get_type ()); + } +} + +void +QR::delete_col (const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + Matrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::delete_col (a, js(i)); + init (a, get_type ()); + } +} + +void +QR::insert_row (const RowVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_row (q*r, j, u), get_type ()); + } +} + +void +QR::delete_row (octave_idx_type j) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_row (q*r, j), get_type ()); + } +} + +void +QR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrshift: index out of range"); + else + { + init (::shift_cols (q*r, i, j), get_type ()); + } +} + +void warn_qrupdate_once (void) +{ + static bool warned = false; + if (! warned) + { + (*current_liboctave_warning_handler) + ("In this version of Octave, QR & Cholesky updating routines\n" + "simply update the matrix and recalculate factorizations.\n" + "To use fast algorithms, link Octave with the qrupdate library.\n" + "See .\n"); + warned = true; + } +} + #endif /* diff --git a/liboctave/dbleQR.h b/liboctave/dbleQR.h --- a/liboctave/dbleQR.h +++ b/liboctave/dbleQR.h @@ -70,7 +70,7 @@ Matrix R (void) const { return r; } -#ifdef HAVE_QRUPDATE + QR::type get_type (void) const; void update (const ColumnVector& u, const ColumnVector& v); @@ -90,8 +90,6 @@ void shift_cols (octave_idx_type i, octave_idx_type j); -#endif - friend std::ostream& operator << (std::ostream&, const QR&); protected: @@ -100,6 +98,10 @@ Matrix r; }; +#ifndef HAVE_QRUPDATE +void warn_qrupdate_once (void); +#endif + #endif /* diff --git a/liboctave/fCmplxCHOL.cc b/liboctave/fCmplxCHOL.cc --- a/liboctave/fCmplxCHOL.cc +++ b/liboctave/fCmplxCHOL.cc @@ -34,6 +34,9 @@ #include "f77-fcn.h" #include "lo-error.h" #include "oct-locbuf.h" +#ifndef HAVE_QRUPDATE +#include "dbleQR.h" +#endif extern "C" { @@ -291,6 +294,146 @@ } } +#else + +void +FloatComplexCHOL::update (const FloatComplexColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + init (chol_mat.hermitian () * chol_mat + + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false); + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); +} + +static bool +singular (const FloatComplexMatrix& a) +{ + for (octave_idx_type i = 0; i < a.rows (); i++) + if (a(i,i) == 0.0f) return true; + return false; +} + +octave_idx_type +FloatComplexCHOL::downdate (const FloatComplexColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + if (singular (chol_mat)) + info = 2; + else + { + info = init (chol_mat.hermitian () * chol_mat + - FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false); + if (info) info = 1; + } + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + + return info; +} + +octave_idx_type +FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); + else + { + if (singular (chol_mat)) + info = 2; + else if (u(j).imag () != 0.0) + info = 3; + else + { + FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; + FloatComplexMatrix a1 (n+1, n+1); + for (octave_idx_type k = 0; k < n+1; k++) + for (octave_idx_type l = 0; l < n+1; l++) + { + if (l == j) + a1(k, l) = u(k); + else if (k == j) + a1(k, l) = std::conj (u(l)); + else + a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); + } + info = init (a1, false); + if (info) info = 1; + } + } + + return info; +} + +void +FloatComplexCHOL::delete_sym (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); + else + { + FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; + a.delete_elements (1, idx_vector (j)); + a.delete_elements (0, idx_vector (j)); + init (a, false); + } +} + +void +FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); + else + { + FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + init (a.index (idx_vector (p), idx_vector (p)), false); + } +} + #endif FloatComplexMatrix diff --git a/liboctave/fCmplxCHOL.h b/liboctave/fCmplxCHOL.h --- a/liboctave/fCmplxCHOL.h +++ b/liboctave/fCmplxCHOL.h @@ -67,8 +67,6 @@ void set (const FloatComplexMatrix& R); -#ifdef HAVE_QRUPDATE - void update (const FloatComplexColumnVector& u); octave_idx_type downdate (const FloatComplexColumnVector& u); @@ -79,8 +77,6 @@ void shift_sym (octave_idx_type i, octave_idx_type j); -#endif - friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatComplexCHOL& a); private: diff --git a/liboctave/fCmplxQR.cc b/liboctave/fCmplxQR.cc --- a/liboctave/fCmplxQR.cc +++ b/liboctave/fCmplxQR.cc @@ -168,14 +168,28 @@ FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& q_arg, const FloatComplexMatrix& r_arg) { - if (q_arg.columns () != r_arg.rows ()) + octave_idx_type qr = q_arg.rows (), qc = q_arg.columns (); + octave_idx_type rr = r_arg.rows (), rc = r_arg.columns (); + if (qc == rr && (qr == qc || (qr > qc && rr == rc))) { - (*current_liboctave_error_handler) ("QR dimensions mismatch"); - return; + q = q_arg; + r = r_arg; } + else + (*current_liboctave_error_handler) ("QR dimensions mismatch"); +} - this->q = q_arg; - this->r = r_arg; +QR::type +FloatComplexQR::get_type (void) const +{ + QR::type retval; + if (!q.is_empty () && q.is_square ()) + retval = QR::std; + else if (q.rows () > q.columns () && r.is_square ()) + retval = QR::economy; + else + retval = QR::raw; + return retval; } #ifdef HAVE_QRUPDATE @@ -196,7 +210,7 @@ utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void @@ -430,6 +444,249 @@ } } +#else + +// Replacement update methods. + +void +FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () == m && v.length () == n) + { + init(q*r + FloatComplexMatrix (u) * FloatComplexMatrix (v).hermitian (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +void +FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + init(q*r + u * v.hermitian (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +static +FloatComplexMatrix insert_col (const FloatComplexMatrix& a, octave_idx_type i, + const FloatComplexColumnVector& x) +{ + FloatComplexMatrix retval (a.rows (), a.columns () + 1); + retval.assign (idx_vector::colon, idx_vector (0, i), + a.index (idx_vector::colon, idx_vector (0, i))); + retval.assign (idx_vector::colon, idx_vector (i), x); + retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()), + a.index (idx_vector::colon, idx_vector (i, a.columns ()))); + return retval; +} + +static +FloatComplexMatrix insert_row (const FloatComplexMatrix& a, octave_idx_type i, + const FloatComplexRowVector& x) +{ + FloatComplexMatrix retval (a.rows () + 1, a.columns ()); + retval.assign (idx_vector (0, i), idx_vector::colon, + a.index (idx_vector (0, i), idx_vector::colon)); + retval.assign (idx_vector (i), idx_vector::colon, x); + retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon, + a.index (idx_vector (i, a.rows ()), idx_vector::colon)); + return retval; +} + +static +FloatComplexMatrix delete_col (const FloatComplexMatrix& a, octave_idx_type i) +{ + FloatComplexMatrix retval = a; + retval.delete_elements (1, idx_vector (i)); + return retval; +} + +static +FloatComplexMatrix delete_row (const FloatComplexMatrix& a, octave_idx_type i) +{ + FloatComplexMatrix retval = a; + retval.delete_elements (0, idx_vector (i)); + return retval; +} + +static +FloatComplexMatrix shift_cols (const FloatComplexMatrix& a, + octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = a.columns (); + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + return a.index (idx_vector::colon, idx_vector (p)); +} + +void +FloatComplexQR::insert_col (const FloatComplexColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_col (q*r, j, u), get_type ()); + } +} + +void +FloatComplexQR::insert_col (const FloatComplexMatrix& u, const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + FloatComplexMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::insert_col (a, js(i), u.column (i)); + init (a, get_type ()); + } +} + +void +FloatComplexQR::delete_col (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_col (q*r, j), get_type ()); + } +} + +void +FloatComplexQR::delete_col (const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + FloatComplexMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::delete_col (a, js(i)); + init (a, get_type ()); + } +} + +void +FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_row (q*r, j, u), get_type ()); + } +} + +void +FloatComplexQR::delete_row (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_row (q*r, j), get_type ()); + } +} + +void +FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrshift: index out of range"); + else + { + init (::shift_cols (q*r, i, j), get_type ()); + } +} + #endif /* diff --git a/liboctave/fCmplxQR.h b/liboctave/fCmplxQR.h --- a/liboctave/fCmplxQR.h +++ b/liboctave/fCmplxQR.h @@ -66,7 +66,7 @@ FloatComplexMatrix R (void) const { return r; } -#ifdef HAVE_QRUPDATE + QR::type get_type (void) const; void update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v); @@ -86,8 +86,6 @@ void shift_cols (octave_idx_type i, octave_idx_type j); -#endif - friend std::ostream& operator << (std::ostream&, const FloatComplexQR&); protected: diff --git a/liboctave/floatCHOL.cc b/liboctave/floatCHOL.cc --- a/liboctave/floatCHOL.cc +++ b/liboctave/floatCHOL.cc @@ -33,6 +33,9 @@ #include "f77-fcn.h" #include "lo-error.h" #include "oct-locbuf.h" +#ifndef HAVE_QRUPDATE +#include "dbleQR.h" +#endif extern "C" { @@ -294,6 +297,144 @@ } } +#else + +void +FloatCHOL::update (const FloatColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + init (chol_mat.transpose () * chol_mat + + FloatMatrix (u) * FloatMatrix (u).transpose (), false); + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); +} + +static bool +singular (const FloatMatrix& a) +{ + for (octave_idx_type i = 0; i < a.rows (); i++) + if (a(i,i) == 0.0f) return true; + return false; +} + +octave_idx_type +FloatCHOL::downdate (const FloatColumnVector& u) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () == n) + { + if (singular (chol_mat)) + info = 2; + else + { + info = init (chol_mat.transpose () * chol_mat + - FloatMatrix (u) * FloatMatrix (u).transpose (), false); + if (info) info = 1; + } + } + else + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); + + return info; +} + +octave_idx_type +FloatCHOL::insert_sym (const FloatColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("cholinsert: index out of range"); + else + { + if (singular (chol_mat)) + info = 2; + else + { + FloatMatrix a = chol_mat.transpose () * chol_mat; + FloatMatrix a1 (n+1, n+1); + for (octave_idx_type k = 0; k < n+1; k++) + for (octave_idx_type l = 0; l < n+1; l++) + { + if (l == j) + a1(k, l) = u(k); + else if (k == j) + a1(k, l) = u(l); + else + a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); + } + info = init (a1, false); + if (info) info = 1; + } + } + + return info; +} + +void +FloatCHOL::delete_sym (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("choldelete: index out of range"); + else + { + FloatMatrix a = chol_mat.transpose () * chol_mat; + a.delete_elements (1, idx_vector (j)); + a.delete_elements (0, idx_vector (j)); + init (a, false); + } +} + +void +FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type n = chol_mat.rows (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("cholshift: index out of range"); + else + { + FloatMatrix a = chol_mat.transpose () * chol_mat; + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + init (a.index (idx_vector (p), idx_vector (p)), false); + } +} + #endif FloatMatrix diff --git a/liboctave/floatCHOL.h b/liboctave/floatCHOL.h --- a/liboctave/floatCHOL.h +++ b/liboctave/floatCHOL.h @@ -64,8 +64,6 @@ void set (const FloatMatrix& R); -#ifdef HAVE_QRUPDATE - void update (const FloatColumnVector& u); octave_idx_type downdate (const FloatColumnVector& u); @@ -76,8 +74,6 @@ void shift_sym (octave_idx_type i, octave_idx_type j); -#endif - friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatCHOL& a); private: diff --git a/liboctave/floatQR.cc b/liboctave/floatQR.cc --- a/liboctave/floatQR.cc +++ b/liboctave/floatQR.cc @@ -159,14 +159,28 @@ FloatQR::FloatQR (const FloatMatrix& q_arg, const FloatMatrix& r_arg) { - if (q_arg.columns () != r_arg.rows ()) + octave_idx_type qr = q_arg.rows (), qc = q_arg.columns (); + octave_idx_type rr = r_arg.rows (), rc = r_arg.columns (); + if (qc == rr && (qr == qc || (qr > qc && rr == rc))) { - (*current_liboctave_error_handler) ("QR dimensions mismatch"); - return; + q = q_arg; + r = r_arg; } + else + (*current_liboctave_error_handler) ("QR dimensions mismatch"); +} - this->q = q_arg; - this->r = r_arg; +QR::type +FloatQR::get_type (void) const +{ + QR::type retval; + if (!q.is_empty () && q.is_square ()) + retval = QR::std; + else if (q.rows () > q.columns () && r.is_square ()) + retval = QR::economy; + else + retval = QR::raw; + return retval; } #ifdef HAVE_QRUPDATE @@ -186,7 +200,7 @@ utmp.fortran_vec (), vtmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void @@ -418,6 +432,249 @@ } } +#else + +// Replacement update methods. + +void +FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () == m && v.length () == n) + { + init(q*r + FloatMatrix (u) * FloatMatrix (v).transpose (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +void +FloatQR::update (const FloatMatrix& u, const FloatMatrix& v) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + init(q*r + u * v.transpose (), get_type ()); + } + else + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); +} + +static +FloatMatrix insert_col (const FloatMatrix& a, octave_idx_type i, + const FloatColumnVector& x) +{ + FloatMatrix retval (a.rows (), a.columns () + 1); + retval.assign (idx_vector::colon, idx_vector (0, i), + a.index (idx_vector::colon, idx_vector (0, i))); + retval.assign (idx_vector::colon, idx_vector (i), x); + retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()), + a.index (idx_vector::colon, idx_vector (i, a.columns ()))); + return retval; +} + +static +FloatMatrix insert_row (const FloatMatrix& a, octave_idx_type i, + const FloatRowVector& x) +{ + FloatMatrix retval (a.rows () + 1, a.columns ()); + retval.assign (idx_vector (0, i), idx_vector::colon, + a.index (idx_vector (0, i), idx_vector::colon)); + retval.assign (idx_vector (i), idx_vector::colon, x); + retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon, + a.index (idx_vector (i, a.rows ()), idx_vector::colon)); + return retval; +} + +static +FloatMatrix delete_col (const FloatMatrix& a, octave_idx_type i) +{ + FloatMatrix retval = a; + retval.delete_elements (1, idx_vector (i)); + return retval; +} + +static +FloatMatrix delete_row (const FloatMatrix& a, octave_idx_type i) +{ + FloatMatrix retval = a; + retval.delete_elements (0, idx_vector (i)); + return retval; +} + +static +FloatMatrix shift_cols (const FloatMatrix& a, + octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = a.columns (); + Array p (n); + for (octave_idx_type k = 0; k < n; k++) p(k) = k; + if (i < j) + { + for (octave_idx_type k = i; k < j; k++) p(k) = k+1; + p(j) = i; + } + else if (j < i) + { + p(j) = i; + for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; + } + + return a.index (idx_vector::colon, idx_vector (p)); +} + +void +FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_col (q*r, j, u), get_type ()); + } +} + +void +FloatQR::insert_col (const FloatMatrix& u, const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + FloatMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::insert_col (a, js(i), u.column (i)); + init (a, get_type ()); + } +} + +void +FloatQR::delete_col (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_col (q*r, j), get_type ()); + } +} + +void +FloatQR::delete_col (const Array& j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + FloatMatrix a = q*r; + for (octave_idx_type i = 0; i < js.length (); i++) + a = ::delete_col (a, js(i)); + init (a, get_type ()); + } +} + +void +FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else + { + init (::insert_row (q*r, j, u), get_type ()); + } +} + +void +FloatQR::delete_row (octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); + else + { + init (::delete_row (q*r, j), get_type ()); + } +} + +void +FloatQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + warn_qrupdate_once (); + + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrshift: index out of range"); + else + { + init (::shift_cols (q*r, i, j), get_type ()); + } +} + #endif diff --git a/liboctave/floatQR.h b/liboctave/floatQR.h --- a/liboctave/floatQR.h +++ b/liboctave/floatQR.h @@ -64,7 +64,7 @@ FloatMatrix R (void) const { return r; } -#ifdef HAVE_QRUPDATE + QR::type get_type (void) const; void update (const FloatColumnVector& u, const FloatColumnVector& v); @@ -84,8 +84,6 @@ void shift_cols (octave_idx_type i, octave_idx_type j); -#endif - friend std::ostream& operator << (std::ostream&, const FloatQR&); protected: diff --git a/scripts/ChangeLog b/scripts/ChangeLog --- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,7 @@ +2009-01-22 Jaroslav Hajek + + * optimization/fsolve.m: Undo the last change. + 2009-01-18 Thorsten Meyer * miscellaneous/doc.m: Add test for existence of info file. diff --git a/scripts/optimization/fsolve.m b/scripts/optimization/fsolve.m --- a/scripts/optimization/fsolve.m +++ b/scripts/optimization/fsolve.m @@ -67,8 +67,6 @@ function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options) - persistent have_qrupdate = exist ('qrupdate') == 5; - if (nargin < 3) options = struct (); endif @@ -268,11 +266,7 @@ endif ## Update the QR factorization. - if (have_qrupdate) - [q, r] = qrupdate (q, r, u, v); - else - [q, r] = qr (q*r + u*v'); - endif + [q, r] = qrupdate (q, r, u, v); endwhile endwhile diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2009-01-22 Jaroslav Hajek + + * DLD-FUNCTIONS/qr.cc: Remove HAVE_QRUPDATE check. + * DLD-FUNCTIONS/chol.cc: Dtto. + 2009-01-21 John W. Eaton * Makefile.in (display.o): Add X11_INCFLAGS to CPPFLAGS. diff --git a/src/DLD-FUNCTIONS/chol.cc b/src/DLD-FUNCTIONS/chol.cc --- a/src/DLD-FUNCTIONS/chol.cc +++ b/src/DLD-FUNCTIONS/chol.cc @@ -576,8 +576,6 @@ return retval; } -#ifdef HAVE_QRUPDATE - DEFUN_DLD (cholupdate, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\ @@ -1276,8 +1274,6 @@ %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) */ -#endif - /* ;;; Local Variables: *** ;;; mode: C++ *** diff --git a/src/DLD-FUNCTIONS/qr.cc b/src/DLD-FUNCTIONS/qr.cc --- a/src/DLD-FUNCTIONS/qr.cc +++ b/src/DLD-FUNCTIONS/qr.cc @@ -739,8 +739,6 @@ */ -#ifdef HAVE_QRUPDATE - static bool check_qr_dims (const octave_value& q, const octave_value& r, bool allow_ecf = false) @@ -1580,8 +1578,6 @@ %! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single')) */ -#endif - /* ;;; Local Variables: *** ;;; mode: C++ ***