Mercurial > octave-nkf
changeset 7559:07522d7dcdf8
fixes to QR and Cholesky updating code
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 05 Mar 2008 14:23:26 -0500 |
| parents | 690c91f741b8 |
| children | 0ef0f9802a37 |
| files | libcruft/ChangeLog libcruft/qrupdate/dch1dn.f libcruft/qrupdate/zch1dn.f liboctave/ChangeLog liboctave/CmplxCHOL.cc liboctave/CmplxQR.cc liboctave/dbleCHOL.cc liboctave/dbleQR.cc src/ChangeLog src/DLD-FUNCTIONS/chol.cc src/DLD-FUNCTIONS/qr.cc |
| diffstat | 11 files changed, 139 insertions(+), 94 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Wed Mar 05 04:44:49 2008 -0500 +++ b/libcruft/ChangeLog Wed Mar 05 14:23:26 2008 -0500 @@ -1,3 +1,7 @@ +2008-03-05 Jaroslav Hajek <highegg@gmail.com> + + * qrupdate/dch1dn.f, qrupdate/zch1dn.f: add "quick return" checks. + 2008-03-04 Jaroslav Hajek <highegg@gmail.com> * qrupdate/dch1dn.f, qrupdate/dch1up.f,
--- a/libcruft/qrupdate/dch1dn.f Wed Mar 05 04:44:49 2008 -0500 +++ b/libcruft/qrupdate/dch1dn.f Wed Mar 05 14:23:26 2008 -0500 @@ -41,6 +41,8 @@ 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
--- a/libcruft/qrupdate/zch1dn.f Wed Mar 05 04:44:49 2008 -0500 +++ b/libcruft/qrupdate/zch1dn.f Wed Mar 05 14:23:26 2008 -0500 @@ -41,6 +41,8 @@ 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
--- a/liboctave/ChangeLog Wed Mar 05 04:44:49 2008 -0500 +++ b/liboctave/ChangeLog Wed Mar 05 14:23:26 2008 -0500 @@ -1,3 +1,11 @@ +2008-03-05 Jaroslav Hajek <highegg@gmail.com> + + * dbleCHOL.cc: Small doc and declaration fixes. + * CmplxHOL.cc: Small doc and 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 <highegg@gmail.com> * dbleCHOL.cc (CHOL::set, CHOL::update, CHOL::downdate):
--- a/liboctave/CmplxCHOL.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/liboctave/CmplxCHOL.cc Wed Mar 05 14:23:26 2008 -0500 @@ -167,7 +167,7 @@ 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 @@ tmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); return info; }
--- a/liboctave/CmplxQR.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/liboctave/CmplxQR.cc Wed Mar 05 14:23:26 2008 -0500 @@ -155,7 +155,6 @@ } } - ComplexQR::ComplexQR (const ComplexMatrix &q, const ComplexMatrix& r) { if (q.columns () != r.rows ()) @@ -175,11 +174,11 @@ 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 @@ 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 {
--- a/liboctave/dbleCHOL.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/liboctave/dbleCHOL.cc Wed Mar 05 14:23:26 2008 -0500 @@ -56,7 +56,7 @@ 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 @@ 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 @@ tmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); return info; }
--- a/liboctave/dbleQR.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/liboctave/dbleQR.cc Wed Mar 05 14:23:26 2008 -0500 @@ -245,7 +245,7 @@ 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 {
--- a/src/ChangeLog Wed Mar 05 04:44:49 2008 -0500 +++ b/src/ChangeLog Wed Mar 05 14:23:26 2008 -0500 @@ -1,3 +1,12 @@ +2008-03-05 Jaroslav Hajek <highegg@gmail.com> + + * 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 <highegg@gmail.com> * DLD-FUNCTIONS/chol.cc (Fcholupdate): New function.
--- a/src/DLD-FUNCTIONS/chol.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/src/DLD-FUNCTIONS/chol.cc Wed Mar 05 14:23:26 2008 -0500 @@ -471,31 +471,35 @@ @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) + { + print_usage (); + return retval; + } - 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 (); + octave_value argr = args(0); + octave_value argu = args(1); - std::string op = (nargin < 3) ? "+" : argop.string_value(); + if (argr.is_matrix_type () && argu.is_matrix_type () + && (nargin < 3 || args(2).is_string ())) + { + octave_idx_type n = argr.rows (); - bool down = false; + std::string op = (nargin < 3) ? "+" : args(2).string_value (); - if (nargin < 3 || (op == "+") || (down = op == "-")) - if (argR.columns () == n - && argu.rows () == n && argu.columns () == 1) + 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 @@ else { // complex case - ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix u = argu.complex_matrix_value (); ComplexCHOL fact;
--- a/src/DLD-FUNCTIONS/qr.cc Wed Mar 05 04:44:49 2008 -0500 +++ b/src/DLD-FUNCTIONS/qr.cc Wed Mar 05 14:23:26 2008 -0500 @@ -448,30 +448,36 @@ octave_idx_type nargin = args.length (); octave_value_list retval; - octave_value argQ,argR,argu,argv; + if (nargin != 4) + { + print_usage (); + return retval; + } - 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 ())) + 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 (); + octave_idx_type m = argq.rows (); + octave_idx_type n = argr.columns (); + octave_idx_type k = argq.columns (); - if (argR.rows () == k + 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 @@ 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 @@ 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(); + 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 (); - bool row = false; + std::string orient = (nargin < 5) ? "col" : args(4).string_value (); - std::string orient = (nargin < 5) ? "col" : argor.string_value (); + bool row = orient == "row"; - if (nargin < 5 || (row = orient == "row") || (orient == "col")) - if (argR.rows () == k + 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 @@ 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 @@ \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\ -@code{columns (Q) == rows (R) = columns (R)}.\n\ +@code{columns (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 @@ @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) + { + print_usage (); + return retval; + } - 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 ()))) - { - octave_idx_type m = argQ.rows (); - octave_idx_type k = argQ.columns (); - octave_idx_type n = argR.columns (); + octave_value argq = args(0); + octave_value argr = args(1); + octave_value argj = args(2); - bool row = false; - bool colp = false; + 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" : argor.string_value (); + std::string orient = (nargin < 4) ? "col" : args(3).string_value (); - if (nargin < 4 || (row = orient == "row") - || (orient == "col") || (colp = orient == "col+")) - if (argR.rows() == k) + 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 @@ 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);