diff src/DLD-FUNCTIONS/qr.cc @ 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 40574114c514
children 0ef0f9802a37
line wrap: on
line diff
--- 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);