# HG changeset patch # User Jaroslav Hajek # Date 1204726701 -3600 # Node ID d917a7de88c7513db6d347320174d925caa0317f # Parent 2ba84879f961fbc023aab0c72c593c70968fe2b2 fixes to QR and Cholesky updating code diff -r 2ba84879f961 -r d917a7de88c7 libcruft/ChangeLog --- a/libcruft/ChangeLog Tue Mar 04 22:58:05 2008 -0500 +++ b/libcruft/ChangeLog Wed Mar 05 15:18:21 2008 +0100 @@ -1,3 +1,7 @@ 2008-03-04 Jaroslav Hajek + + * qrupdate/dch1dn.f, qrupdate/zch1dn.f: add "quick return" checks. + 2008-03-04 Jaroslav Hajek * qrupdate/dch1dn.f, qrupdate/dch1up.f, diff -r 2ba84879f961 -r d917a7de88c7 libcruft/qrupdate/dch1dn.f --- a/libcruft/qrupdate/dch1dn.f Tue Mar 04 22:58:05 2008 -0500 +++ b/libcruft/qrupdate/dch1dn.f Wed Mar 05 15:18:21 2008 +0100 @@ -41,6 +41,8 @@ c double precision rr,ui,t integer i,j +c quick return if possible + if (n <= 0) return c check for singularity of R do i = 1,n if (R(i,i) == 0d0) then diff -r 2ba84879f961 -r d917a7de88c7 libcruft/qrupdate/zch1dn.f --- a/libcruft/qrupdate/zch1dn.f Tue Mar 04 22:58:05 2008 -0500 +++ b/libcruft/qrupdate/zch1dn.f Wed Mar 05 15:18:21 2008 +0100 @@ -41,6 +41,8 @@ c double complex crho,rr,ui,t integer i,j +c quick return if possible + if (n <= 0) return c check for singularity of R do i = 1,n if (R(i,i) == 0d0) then diff -r 2ba84879f961 -r d917a7de88c7 liboctave/ChangeLog --- a/liboctave/ChangeLog Tue Mar 04 22:58:05 2008 -0500 +++ b/liboctave/ChangeLog Wed Mar 05 15:18:21 2008 +0100 @@ -1,3 +1,11 @@ 2008-03-04 Jaroslav Hajek + + * dbleCHOL.cc: Small doc & declaration fixes. + * CmplxHOL.cc: Small doc & declaration fixes. + * CmplxQR.cc (ComplexQR::ComplexQR): adjust code to match dbleQR.cc. + * dbleQR.cc (QR::delete_row): fix incorrect test. + * CmplxQR.cc (ComplexQR::delete_row): fix incorrect test. + 2008-03-04 Jaroslav Hajek * dbleCHOL.cc (CHOL::set, CHOL::update, CHOL::downdate): diff -r 2ba84879f961 -r d917a7de88c7 liboctave/CmplxCHOL.cc --- a/liboctave/CmplxCHOL.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/liboctave/CmplxCHOL.cc Wed Mar 05 15:18:21 2008 +0100 @@ -167,7 +167,7 @@ ComplexCHOL::set (const ComplexMatrix& R if (R.is_square ()) chol_mat = R; else - (*current_liboctave_error_handler) ("chol2inv requires square matrix"); + (*current_liboctave_error_handler) ("CHOL requires square matrix"); } void @@ -205,7 +205,7 @@ ComplexCHOL::downdate (const ComplexMatr tmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); return info; } diff -r 2ba84879f961 -r d917a7de88c7 liboctave/CmplxQR.cc --- a/liboctave/CmplxQR.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/liboctave/CmplxQR.cc Wed Mar 05 15:18:21 2008 +0100 @@ -155,7 +155,6 @@ ComplexQR::init (const ComplexMatrix& a, } } - ComplexQR::ComplexQR (const ComplexMatrix &q, const ComplexMatrix& r) { if (q.columns () != r.rows ()) @@ -175,11 +174,11 @@ ComplexQR::update (const ComplexMatrix& octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); - if (u.length () != m || v.length () != n) - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); - else + if (u.length () == m && v.length () == n) F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), u.data (), v.data ())); + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); } void @@ -257,7 +256,7 @@ ComplexQR::delete_row (octave_idx_type j if (! q.is_square ()) (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); - else if (j < 1 || j > n) + else if (j < 1 || j > m) (*current_liboctave_error_handler) ("QR delete index out of range"); else { diff -r 2ba84879f961 -r d917a7de88c7 liboctave/dbleCHOL.cc --- a/liboctave/dbleCHOL.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/liboctave/dbleCHOL.cc Wed Mar 05 15:18:21 2008 +0100 @@ -56,7 +56,7 @@ extern "C" F77_RET_T F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, - int &); + octave_idx_type&); } octave_idx_type @@ -171,13 +171,13 @@ CHOL::set (const Matrix& R) if (R.is_square ()) chol_mat = R; else - (*current_liboctave_error_handler) ("chol2inv requires square matrix"); + (*current_liboctave_error_handler) ("CHOL requires square matrix"); } void CHOL::update (const Matrix& u) { - int n = chol_mat.rows (); + octave_idx_type n = chol_mat.rows (); if (u.length () == n) { @@ -209,7 +209,7 @@ CHOL::downdate (const Matrix& u) tmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); return info; } diff -r 2ba84879f961 -r d917a7de88c7 liboctave/dbleQR.cc --- a/liboctave/dbleQR.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/liboctave/dbleQR.cc Wed Mar 05 15:18:21 2008 +0100 @@ -245,7 +245,7 @@ QR::delete_row (octave_idx_type j) if (! q.is_square ()) (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); - else if (j < 1 || j > n) + else if (j < 1 || j > m) (*current_liboctave_error_handler) ("QR delete index out of range"); else { diff -r 2ba84879f961 -r d917a7de88c7 src/ChangeLog --- a/src/ChangeLog Tue Mar 04 22:58:05 2008 -0500 +++ b/src/ChangeLog Wed Mar 05 15:18:21 2008 +0100 @@ -1,3 +1,12 @@ 2008-03-04 Jaroslav Hajek + + * DLD-FUNCTIONS/chol.cc (Fcholupdate): Adjust code to meet + Octave's coding guidelines. + + * DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrdelete): Adjust + code to meet Octave's coding guidelines. + * DLD-FUNCTIONS/qr.cc (Fqrdelete): Fix incorrect test. + 2008-03-04 Jaroslav Hajek * DLD-FUNCTIONS/chol.cc (Fcholupdate): New function. diff -r 2ba84879f961 -r d917a7de88c7 src/DLD-FUNCTIONS/chol.cc --- a/src/DLD-FUNCTIONS/chol.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/src/DLD-FUNCTIONS/chol.cc Wed Mar 05 15:18:21 2008 +0100 @@ -471,31 +471,35 @@ If @var{info} is not present, an error m @seealso{chol, qrupdate}\n\ @end deftypefn") { - int nargin = args.length (); + octave_idx_type nargin = args.length (); octave_value_list retval; - octave_value argR,argu,argop; - - if ((nargin == 3 || nargin == 2) - && (argR = args(0), argR.is_matrix_type ()) - && (argu = args(1), argu.is_matrix_type ()) - && (nargin < 3 || (argop = args(2), argop.is_string ()))) - { - octave_idx_type n = argR.rows (); - - std::string op = (nargin < 3) ? "+" : argop.string_value(); - - bool down = false; - - if (nargin < 3 || (op == "+") || (down = op == "-")) - if (argR.columns () == n - && argu.rows () == n && argu.columns () == 1) + if (nargin > 3 || nargin < 2) + { + print_usage (); + return retval; + } + + octave_value argr = args(0); + octave_value argu = args(1); + + if (argr.is_matrix_type () && argu.is_matrix_type () + && (nargin < 3 || args(2).is_string ())) + { + octave_idx_type n = argr.rows (); + + std::string op = (nargin < 3) ? "+" : args(2).string_value (); + + bool down = op == "-"; + + if (down || op == "+") + if (argr.columns () == n && argu.rows () == n && argu.columns () == 1) { - if (argR.is_real_matrix () && argu.is_real_matrix ()) + if (argr.is_real_matrix () && argu.is_real_matrix ()) { // real case - Matrix R = argR.matrix_value (); + Matrix R = argr.matrix_value (); Matrix u = argu.matrix_value (); CHOL fact; @@ -517,7 +521,7 @@ If @var{info} is not present, an error m else { // complex case - ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix u = argu.complex_matrix_value (); ComplexCHOL fact; diff -r 2ba84879f961 -r d917a7de88c7 src/DLD-FUNCTIONS/qr.cc --- a/src/DLD-FUNCTIONS/qr.cc Tue Mar 04 22:58:05 2008 -0500 +++ b/src/DLD-FUNCTIONS/qr.cc Wed Mar 05 15:18:21 2008 +0100 @@ -448,30 +448,36 @@ Q*Q'*u*v' instead of u*v'.\n\ octave_idx_type nargin = args.length (); octave_value_list retval; - octave_value argQ,argR,argu,argv; - - if (nargin == 4 - && (argQ = args(0),argQ.is_matrix_type ()) - && (argR = args(1),argR.is_matrix_type ()) - && (argu = args(2),argu.is_matrix_type ()) - && (argv = args(3),argv.is_matrix_type ())) + if (nargin != 4) { - octave_idx_type m = argQ.rows (); - octave_idx_type n = argR.columns (); - octave_idx_type k = argQ.columns (); - - if (argR.rows () == k + print_usage (); + return retval; + } + + octave_value argq = args(0); + octave_value argr = args(1); + octave_value argu = args(2); + octave_value argv = args(3); + + if (argq.is_matrix_type () && argr.is_matrix_type () + && argu.is_matrix_type () && argv.is_matrix_type ()) + { + octave_idx_type m = argq.rows (); + octave_idx_type n = argr.columns (); + octave_idx_type k = argq.columns (); + + if (argr.rows () == k && argu.rows () == m && argu.columns () == 1 && argv.rows () == n && argv.columns () == 1) { - if (argQ.is_real_matrix () - && argR.is_real_matrix () + if (argq.is_real_matrix () + && argr.is_real_matrix () && argu.is_real_matrix () && argv.is_real_matrix ()) { // all real case - Matrix Q = argQ.matrix_value (); - Matrix R = argR.matrix_value (); + Matrix Q = argq.matrix_value (); + Matrix R = argr.matrix_value (); Matrix u = argu.matrix_value (); Matrix v = argv.matrix_value (); @@ -484,8 +490,8 @@ Q*Q'*u*v' instead of u*v'.\n\ else { // complex case - ComplexMatrix Q = argQ.complex_matrix_value (); - ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix Q = argq.complex_matrix_value (); + ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix u = argu.complex_matrix_value (); ComplexMatrix v = argv.complex_matrix_value (); @@ -577,48 +583,54 @@ If @var{orient} is @code{\"row\"}, @var{ octave_idx_type nargin = args.length (); octave_value_list retval; - octave_value argQ,argR,argj,argx,argor; - - if ((nargin == 4 || nargin == 5) - && (argQ = args(0), argQ.is_matrix_type ()) - && (argR = args(1), argR.is_matrix_type ()) - && (argj = args(2), argj.is_scalar_type ()) - && (argx = args(3), argx.is_matrix_type ()) - && (nargin < 5 || (argor = args (4), argor.is_string ()))) + if (nargin < 4 || nargin > 5) { - octave_idx_type m = argQ.rows (); - octave_idx_type n = argR.columns (); - octave_idx_type k = argQ.columns(); - - bool row = false; - - std::string orient = (nargin < 5) ? "col" : argor.string_value (); - - if (nargin < 5 || (row = orient == "row") || (orient == "col")) - if (argR.rows () == k + print_usage (); + return retval; + } + + octave_value argq = args(0); + octave_value argr = args(1); + octave_value argj = args(2); + octave_value argx = args(3); + + if (argq.is_matrix_type () && argr.is_matrix_type () + && argj.is_scalar_type () && argx.is_matrix_type () + && (nargin < 5 || args(4).is_string ())) + { + octave_idx_type m = argq.rows (); + octave_idx_type n = argr.columns (); + octave_idx_type k = argq.columns (); + + std::string orient = (nargin < 5) ? "col" : args(4).string_value (); + + bool row = orient == "row"; + + if (row || orient == "col") + if (argr.rows () == k && (! row || m == k) && argx.rows () == (row ? 1 : m) && argx.columns () == (row ? n : 1)) { - int j = argj.int_value (); + octave_idx_type j = argj.int_value (); if (j >= 1 && j <= (row ? n : m)+1) { - if (argQ.is_real_matrix () - && argR.is_real_matrix () + if (argq.is_real_matrix () + && argr.is_real_matrix () && argx.is_real_matrix ()) { // real case - Matrix Q = argQ.matrix_value (); - Matrix R = argR.matrix_value (); + Matrix Q = argq.matrix_value (); + Matrix R = argr.matrix_value (); Matrix x = argx.matrix_value (); QR fact (Q, R); if (row) - fact.insert_row(x, j); + fact.insert_row (x, j); else - fact.insert_col(x, j); + fact.insert_col (x, j); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -626,8 +638,8 @@ If @var{orient} is @code{\"row\"}, @var{ else { // complex case - ComplexMatrix Q = argQ.complex_matrix_value (); - ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix Q = argq.complex_matrix_value (); + ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix x = argx.complex_matrix_value (); ComplexQR fact (Q, R); @@ -745,7 +757,7 @@ be square.\n\ \n\ For @sc{Matlab} compatibility, if @var{Q} is nonsquare on input, the\n\ updated factorization is always stripped to the economical form, i.e.\n\ address@hidden (Q) == rows (R) = columns (R)}.\n\ address@hidden (Q) == rows (R) <= columns (R)}.\n\ \n\ To get the less intelligent but more natural behaviour when @var{Q}\n\ retains it shape and @var{R} loses one column, set @var{orient} to\n\ @@ -755,39 +767,44 @@ If @var{orient} is \"row\", @var{Q} must @seealso{qr, qrinsert, qrupdate}\n\ @end deftypefn") { - int nargin = args.length (); + octave_idx_type nargin = args.length (); octave_value_list retval; - octave_value argQ,argR,argj,argor; - - if ((nargin == 3 || nargin == 4) - && (argQ = args(0), argQ.is_matrix_type ()) - && (argR = args(1), argR.is_matrix_type ()) - && (argj = args(2), argj.is_scalar_type ()) - && (nargin < 4 || (argor = args (3), argor.is_string ()))) + if (nargin < 3 || nargin > 4) { - octave_idx_type m = argQ.rows (); - octave_idx_type k = argQ.columns (); - octave_idx_type n = argR.columns (); - - bool row = false; - bool colp = false; - - std::string orient = (nargin < 4) ? "col" : argor.string_value (); - - if (nargin < 4 || (row = orient == "row") - || (orient == "col") || (colp = orient == "col+")) - if (argR.rows() == k) + print_usage (); + return retval; + } + + octave_value argq = args(0); + octave_value argr = args(1); + octave_value argj = args(2); + + if (argq.is_matrix_type () && argr.is_matrix_type () && argj.is_scalar_type () + && (nargin < 4 || args(3).is_string ())) + { + octave_idx_type m = argq.rows (); + octave_idx_type k = argq.columns (); + octave_idx_type n = argr.columns (); + + std::string orient = (nargin < 4) ? "col" : args(3).string_value (); + + bool row = orient == "row"; + bool colp = orient == "col+"; + + if (row || colp || orient == "col") + if (argr.rows () == k + && (! row || m == k)) { - int j = argj.scalar_value (); - if (j >= 1 && j <= (row ? n : m)+1) + octave_idx_type j = argj.scalar_value (); + if (j >= 1 && j <= (row ? n : m)) { - if (argQ.is_real_matrix () - && argR.is_real_matrix ()) + if (argq.is_real_matrix () + && argr.is_real_matrix ()) { // real case - Matrix Q = argQ.matrix_value (); - Matrix R = argR.matrix_value (); + Matrix Q = argq.matrix_value (); + Matrix R = argr.matrix_value (); QR fact (Q, R); @@ -807,8 +824,8 @@ If @var{orient} is \"row\", @var{Q} must else { // complex case - ComplexMatrix Q = argQ.complex_matrix_value (); - ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix Q = argq.complex_matrix_value (); + ComplexMatrix R = argr.complex_matrix_value (); ComplexQR fact (Q, R);