# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1235549438 -3600
# Node ID c3b743b1b1c6ffd8047f9ddbdc462681f5c8a7c0
# Parent  4d812facab0e6bca6aacb7a8a42d634335f4b28a
preset triangular type if possible for lu and qr outputs

diff -r 4d812facab0e -r c3b743b1b1c6 src/ChangeLog
--- 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
diff -r 4d812facab0e -r c3b743b1b1c6 src/DLD-FUNCTIONS/lu.cc
--- 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;
 		    }
diff -r 4d812facab0e -r c3b743b1b1c6 src/DLD-FUNCTIONS/qr.cc
--- 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;