Mercurial > octave-nkf
changeset 8869:c3b743b1b1c6
preset triangular type if possible for lu and qr outputs
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 25 Feb 2009 09:10:38 +0100 |
| parents | 4d812facab0e |
| children | eea0e1b45ec0 |
| files | src/ChangeLog src/DLD-FUNCTIONS/lu.cc src/DLD-FUNCTIONS/qr.cc |
| diffstat | 3 files changed, 73 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog Wed Feb 25 03:02:28 2009 -0500 +++ b/src/ChangeLog Wed Feb 25 09:10:38 2009 +0100 @@ -1,3 +1,10 @@ +2009-02-25 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/lu.cc (maybe_set_triangular): New function. + (Flu): Use it. + * DLD-FUNCTIONS/qr.cc (maybe_set_triangular): New function. + (Fqr): Use it. + 2009-02-25 John W. Eaton <jwe@octave.org> * input.cc (get_debug_input): Write debugging location info
--- a/src/DLD-FUNCTIONS/lu.cc Wed Feb 25 03:02:28 2009 -0500 +++ b/src/DLD-FUNCTIONS/lu.cc Wed Feb 25 09:10:38 2009 +0100 @@ -40,6 +40,29 @@ #include "ov-re-sparse.h" #include "ov-cx-sparse.h" +template <class MT> +static octave_value +maybe_set_triangular (const MT& m, MatrixType::matrix_type t = MatrixType::Upper) +{ + typedef typename MT::element_type T; + octave_value retval; + octave_idx_type r = m.rows (), c = m.columns (); + if (r == c) + { + const T zero = T(); + octave_idx_type i = 0; + for (;i != r && m(i,i) != zero; i++) ; + if (i == r) + retval = octave_value (m, MatrixType (t)); + else + retval = m; + } + else + retval = m; + + return retval; +} + DEFUN_DLD (lu, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =} lu (@var{a})\n\ @@ -358,7 +381,7 @@ { PermMatrix P = fact.P (); FloatMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); retval(0) = L; } break; @@ -370,8 +393,8 @@ retval(2) = fact.P_vec (); else retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); + retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower); } break; } @@ -396,7 +419,7 @@ { PermMatrix P = fact.P (); Matrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); retval(0) = L; } break; @@ -408,8 +431,8 @@ retval(2) = fact.P_vec (); else retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); + retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower); } break; } @@ -437,7 +460,7 @@ { PermMatrix P = fact.P (); FloatComplexMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); retval(0) = L; } break; @@ -449,8 +472,8 @@ retval(2) = fact.P_vec (); else retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); + retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower); } break; } @@ -475,7 +498,7 @@ { PermMatrix P = fact.P (); ComplexMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); retval(0) = L; } break; @@ -487,8 +510,8 @@ retval(2) = fact.P_vec (); else retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper); + retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower); } break; }
--- a/src/DLD-FUNCTIONS/qr.cc Wed Feb 25 03:02:28 2009 -0500 +++ b/src/DLD-FUNCTIONS/qr.cc Wed Feb 25 09:10:38 2009 +0100 @@ -44,6 +44,29 @@ #include "oct-obj.h" #include "utils.h" +template <class MT> +static octave_value +maybe_set_triangular (const MT& m, MatrixType::matrix_type t = MatrixType::Upper) +{ + typedef typename MT::element_type T; + octave_value retval; + octave_idx_type r = m.rows (), c = m.columns (); + if (r == c) + { + const T zero = T(); + octave_idx_type i = 0; + for (;i != r && m(i,i) != zero; i++) ; + if (i == r) + retval = octave_value (m, MatrixType (t)); + else + retval = m; + } + else + retval = m; + + return retval; +} + // [Q, R] = qr (X): form Q unitary and R upper triangular such // that Q * R = X // @@ -296,7 +319,7 @@ case 2: { FloatQR fact (m, type); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -308,7 +331,7 @@ retval(2) = fact.Pvec (); else retval(2) = fact.P (); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -334,7 +357,7 @@ case 2: { FloatComplexQR fact (m, type); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -346,7 +369,7 @@ retval(2) = fact.Pvec (); else retval(2) = fact.P (); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -375,7 +398,7 @@ case 2: { QR fact (m, type); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -387,7 +410,7 @@ retval(2) = fact.Pvec (); else retval(2) = fact.P (); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -413,7 +436,7 @@ case 2: { ComplexQR fact (m, type); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break; @@ -425,7 +448,7 @@ retval(2) = fact.Pvec (); else retval(2) = fact.P (); - retval(1) = fact.R (); + retval(1) = maybe_set_triangular (fact.R ()); retval(0) = fact.Q (); } break;