changeset 3883:69b6bd271277

[project @ 2002-04-02 21:05:10 by jwe]
author jwe
date Tue, 02 Apr 2002 21:05:10 +0000
parents c8c1ead8474f
children fcb4931ec48a
files liboctave/ChangeLog liboctave/CmplxQR.cc liboctave/CmplxQRP.cc liboctave/dbleQR.cc liboctave/dbleQRP.cc scripts/ChangeLog scripts/linear-algebra/cross.m src/ChangeLog src/lex.h src/lex.l src/parse.y src/toplev.cc test/octave.test/linalg/linalg.exp
diffstat 13 files changed, 180 insertions(+), 141 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,11 @@
+2002-04-02  Paul Kienzle <pkienzle@users.sf.net>
+
+        * CmplxQR.cc (ComplexQR::init): Use economy QR decomposition
+	internally when the user requests it.
+	* CmplxQRP.cc (ComplexQRP::init): Ditto.
+	* dbleQR.cc (QR::init): Ditto.
+	* dbleQRP.cc (QRP::init): Ditto.
+
 2002-02-22  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* oct-fftw.cc (octave_fftw::fft2d): Avoid having to find a
--- a/liboctave/CmplxQR.cc
+++ b/liboctave/CmplxQR.cc
@@ -74,7 +74,7 @@
   int info = 0;
 
   ComplexMatrix A_fact;
-  if (m > n)
+  if (m > n && qr_type != QR::economy)
     {
       A_fact.resize (m, m);
       A_fact.insert (a, 0, 0);
@@ -106,18 +106,12 @@
 	}
       else
 	{
-	  volatile int n2;
+	  int n2 = (qr_type == QR::economy) ? min_mn : m;
 
 	  if (qr_type == QR::economy && m > n)
-	    {
-	      n2 = n;
-	      r.resize (n, n, 0.0);
-	    }
+	    r.resize (n, n, 0.0);
 	  else
-	    {
-	      n2 = m;
-	      r.resize (m, n, 0.0);
-	    }
+	    r.resize (m, n, 0.0);
 
 	  for (int j = 0; j < n; j++)
 	    {
@@ -126,11 +120,11 @@
 		r.elem (i, j) = A_fact.elem (i, j);
 	    }
 
-	  lwork = 32*m;
+	  lwork = 32 * n2;
 	  work.resize (lwork);
 	  Complex *pwork = work.fortran_vec ();
 
-	  F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, ptau,
+	  F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				     pwork, lwork, info));
 
 	  if (f77_exception_encountered)
--- a/liboctave/CmplxQRP.cc
+++ b/liboctave/CmplxQRP.cc
@@ -68,7 +68,8 @@
       return;
     }
 
-  Array<Complex> tau (m < n ? m : n);
+  int min_mn = m < n ? m : n;
+  Array<Complex> tau (min_mn);
   Complex *ptau = tau.fortran_vec ();
 
   int lwork = 3*n > 32*m ? 3*n : 32*m;
@@ -78,7 +79,7 @@
   int info = 0;
 
   ComplexMatrix A_fact = a;
-  if (m > n)
+  if (m > n && qr_type != QR::economy)
     A_fact.resize (m, m, 0.0);
 
   Complex *tmp_data = A_fact.fortran_vec ();
@@ -114,13 +115,13 @@
 	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
 	}
 
+      int n2 = (qr_type == QR::economy) ? min_mn : m;
+
       if (qr_type == QR::economy && m > n)
 	r.resize (n, n, 0.0);
       else
 	r.resize (m, n, 0.0);
 
-      int min_mn = m < n ? m : n;
-
       for (int j = 0; j < n; j++)
 	{
 	  int limit = j < min_mn-1 ? j : min_mn-1;
@@ -128,10 +129,6 @@
 	    r.elem (i, j) = A_fact.elem (i, j);
 	}
 
-      int n2 = m;
-      if (qr_type == QR::economy)
-	n2 = min_mn;
-
       F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				 pwork, lwork, info));
 
--- a/liboctave/dbleQR.cc
+++ b/liboctave/dbleQR.cc
@@ -71,14 +71,9 @@
 
   int info = 0;
 
-  Matrix A_fact;
-  if (m > n)
-    {
-      A_fact.resize (m, m);
-      A_fact.insert (a, 0, 0);
-    }
-  else
-    A_fact = a;
+  Matrix A_fact = a;
+  if (m > n && qr_type != QR::economy)
+      A_fact.resize (m, m, 0.0);
 
   double *tmp_data = A_fact.fortran_vec ();
 
@@ -104,18 +99,12 @@
 	}
       else
 	{
-	  volatile int n2;
+	  int n2 = (qr_type == QR::economy) ? min_mn : m;
 
 	  if (qr_type == QR::economy && m > n)
-	    {
-	      n2 = n;
-	      r.resize (n, n, 0.0);
-	    }
+	    r.resize (n, n, 0.0);
 	  else
-	    {
-	      n2 = m;
-	      r.resize (m, n, 0.0);
-	    }
+	    r.resize (m, n, 0.0);
 
 	  for (int j = 0; j < n; j++)
 	    {
@@ -124,11 +113,11 @@
 		r.elem (i, j) = tmp_data[m*j+i];
 	    }
 
-	  lwork = 32*m;
+	  lwork = 32 * n2;
 	  work.resize (lwork);
 	  double *pwork = work.fortran_vec ();
 
-	  F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, ptau,
+	  F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				     pwork, lwork, info));
 
 	  if (f77_exception_encountered)
--- a/liboctave/dbleQRP.cc
+++ b/liboctave/dbleQRP.cc
@@ -67,7 +67,8 @@
       return;
     }
 
-  Array<double> tau (m < n ? m : n);
+  int min_mn = m < n ? m : n;
+  Array<double> tau (min_mn);
   double *ptau = tau.fortran_vec ();
 
   int lwork = 3*n > 32*m ? 3*n : 32*m;
@@ -77,7 +78,7 @@
   int info = 0;
 
   Matrix A_fact = a;
-  if (m > n)
+  if (m > n && qr_type != QR::economy)
     A_fact.resize (m, m, 0.0);
 
   double *tmp_data = A_fact.fortran_vec ();
@@ -109,13 +110,13 @@
 	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
 	}
 
+      int n2 = (qr_type == QR::economy) ? min_mn : m;
+
       if (qr_type == QR::economy && m > n)
 	r.resize (n, n, 0.0);
       else
 	r.resize (m, n, 0.0);
 
-      int min_mn = m < n ? m : n;
-
       for (int j = 0; j < n; j++)
 	{
 	  int limit = j < min_mn-1 ? j : min_mn-1;
@@ -123,10 +124,6 @@
 	    r.elem (i, j) = A_fact.elem (i, j);
 	}
 
-      int n2 = m;
-      if (qr_type == QR::economy)
-	n2 = min_mn;
-
       F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				 pwork, lwork, info));
 
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,5 +1,9 @@
 2002-04-02  Paul Kienzle <pkienzle@users.sf.net>
 
+	* linear-algebra/cross.m: Accept nx3 and 3xn matrices, in addition
+	to vectors.  Issue a warning in the case x matches y' but return a
+	column vector as Octave currently does.
+
 	* plot/contour.m: Set default number of levels for contour(x,y,z).
 
 	* control/system/starp.m: Leave more of the documentation
--- a/scripts/linear-algebra/cross.m
+++ b/scripts/linear-algebra/cross.m
@@ -19,40 +19,53 @@
 ## Computes the vector cross product of the two 3-dimensional vectors
 ## @var{x} and @var{y}.
 ##
-## A row vector is returned if @var{x} and @var{y} are both row vectors;
-## otherwise, a column vector is returned.
-##
 ## @example
 ## @group
 ## cross ([1,1,0], [0,1,1])
 ##      @result{} [ 1; -1; 1 ]
 ## @end group
 ## @end example
+##
+## If @var{x} and @var{y} are two - dimensional matrices the
+## cross product is applied along the first dimension with 3 elements.
+##
 ## @end deftypefn
 
-## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
+## Author: Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at>
 ## Created: 15 October 1994
 ## Adapted-By: jwe
 
 function z = cross (x, y)
-
+	
   if (nargin != 2)
     usage ("cross (x, y)");
   endif
 
-  if (length (x) == 3 && length (y) == 3)
-
-    z = [x(2)*y(3) - x(3)*y(2); x(3)*y(1) - x(1)*y(3); x(1)*y(2) - x(2)*y(1)];
-
-    x_nr = rows (x);
-    y_nr = rows (y);
+  ## XXX COMPATIBILITY XXX opposite behaviour for cross(row,col)
+  ## Swap x and y in the assignments below to get the matlab behaviour.
+  ## Better yet, fix the calling code so that it uses conformant vectors.
+  if (columns(x) == 1 && rows(y) == 1)
+    warning ("cross: taking cross product of column by row");
+    y = y.';
+  elseif (rows(x) == 1 && columns(y) == 1)
+    warning ("cross: taking cross product of row by column");
+    x = x.';
+  endif
 
-    if (x_nr == y_nr && x_nr == 1)
-      z = z.';
+  if (size(x) == size(y))
+    if (rows(x) == 3)
+      z = [ ( x (2,:) .* y (3,:) - x (3,:) .* y (2,:) ) ;
+            ( x (3,:) .* y (1,:) - x (1,:) .* y (3,:) ) ;
+            ( x (1,:) .* y (2,:) - x (2,:) .* y (1,:) ) ];
+    elseif (columns(x) == 3)
+      z = [ ( x (:,2) .* y (:,3) - x (:,3) .* y (:,2) ) , \
+            ( x (:,3) .* y (:,1) - x (:,1) .* y (:,3) ) , \
+            ( x (:,1) .* y (:,2) - x (:,2) .* y (:,1) ) ];
+    else
+      error ("cross: x,y must have dimension nx3 or 3xn");
     endif
-
   else
-    error ("cross: both x and y must be 3-dimensional vectors");
+    error ("cross: x and y must have the same dimensions");
   endif
 
 endfunction
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,7 +1,19 @@
 2002-04-02  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* lex.l, lex.h (parser_end_of_input): New global variable.
+	(reset_parser): Reset it here.
+	* parse.y (input): Set it to TRUE on EOF.
+	(parse_and_executed): Save and restore it, check it to correctly
+	break out of parse-execute loop.
+	(parse_fcn_file): Likewise.
+	(eval_string): Likewise.
+	* toplev.cc (main_loop): Likewise.
+	
+
 	* parse.y (input): Call YYACCEPT for END_OF_INPUT.
+	Return no-op command for bare newline.
 	(parse_and_execute): Handle change in yyparse return value semantics.
+
 	* toplev.cc (main_loop): Likewise.
 
 2002-03-25  John W. Eaton  <jwe@bevo.che.wisc.edu>
--- a/src/lex.h
+++ b/src/lex.h
@@ -197,6 +197,9 @@
   lexical_feedback& operator = (const lexical_feedback&);
 };
 
+// TRUE means that we have encountered EOF on the input stream.
+extern bool parser_end_of_input;
+
 // Flags that need to be shared between the lexer and parser.
 extern lexical_feedback lexer_flags;
 
--- a/src/lex.l
+++ b/src/lex.l
@@ -66,6 +66,9 @@
 #error lex.l requires flex version 2.5.4 or later
 #endif
 
+// TRUE means that we have encountered EOF on the input stream.
+bool parser_end_of_input = false;
+
 // Flags that need to be shared between the lexer and parser.
 lexical_feedback lexer_flags;
 
@@ -761,6 +764,7 @@
   BEGIN 0;
   error_state = 0;
   warning_state = 0;
+  parser_end_of_input = false;
 
   // We do want a prompt by default.
   promptflag = 1;
--- a/src/parse.y
+++ b/src/parse.y
@@ -462,12 +462,6 @@
 		    promptflag = 1;
 		    YYACCEPT;
 		  }
-		| END_OF_INPUT
-		  {
-		    global_command = 0;
-		    promptflag = 1;
-		    YYACCEPT;
-		  }
 		| simple_list parse_error
 		  { ABORT_PARSE; }
 		| parse_error
@@ -476,6 +470,11 @@
 
 input1		: '\n'
 		  { $$ = 0; }
+		| END_OF_INPUT
+		  {
+		    parser_end_of_input = 1;
+		    $$ = 0;
+		  }
 		| simple_list
 		  { $$ = $1; }
 		| simple_list '\n'
@@ -2732,10 +2731,12 @@
   unwind_protect_bool (line_editing);
   unwind_protect_bool (input_from_command_line_file);
   unwind_protect_bool (get_input_from_eval_string);
+  unwind_protect_bool (parser_end_of_input);
 
   line_editing = false;
   input_from_command_line_file = false;
   get_input_from_eval_string = false;
+  parser_end_of_input = false;
 
   unwind_protect_ptr (curr_sym_tab);
 
@@ -2746,36 +2747,39 @@
 
       retval = yyparse ();
 
-      if (retval == 0 && global_command)
-	{
-	  global_command->eval ();
-
-	  delete global_command;
-
-	  global_command = 0;
-
-	  bool quit = (tree_return_command::returning
-		       || tree_break_command::breaking);
-
-	  if (tree_return_command::returning)
-	    tree_return_command::returning = 0;
-
-	  if (tree_break_command::breaking)
-	    tree_break_command::breaking--;
-
-	  if (error_state)
+      if (retval == 0)
+        {
+          if (global_command)
 	    {
-	      error ("near line %d of file `%s'", input_line_number,
-		     curr_fcn_file_full_name.c_str ());
-
-	      break;
+	      global_command->eval ();
+
+	      delete global_command;
+
+	      global_command = 0;
+
+	      bool quit = (tree_return_command::returning
+			   || tree_break_command::breaking);
+
+	      if (tree_return_command::returning)
+		tree_return_command::returning = 0;
+
+	      if (tree_break_command::breaking)
+		tree_break_command::breaking--;
+
+	      if (error_state)
+		{
+		  error ("near line %d of file `%s'", input_line_number,
+			 curr_fcn_file_full_name.c_str ());
+
+		  break;
+		}
+
+	      if (quit)
+		break;
 	    }
-
-	  if (quit)
+	  else if (parser_end_of_input)
 	    break;
-	}
-      else
-	break;
+        }
     }
   while (retval == 0);
 
@@ -3141,12 +3145,14 @@
 	  unwind_protect_bool (reading_fcn_file);
 	  unwind_protect_bool (input_from_command_line_file);
 	  unwind_protect_bool (get_input_from_eval_string);
+	  unwind_protect_bool (parser_end_of_input);
 
 	  Vecho_executing_commands = ECHO_OFF;
 	  Vsaving_history = false;
 	  reading_fcn_file = true;
 	  input_from_command_line_file = false;
 	  get_input_from_eval_string = false;
+	  parser_end_of_input = false;
 
 	  YY_BUFFER_STATE old_buf = current_buffer ();
 	  YY_BUFFER_STATE new_buf = create_buffer (ffile);
@@ -3406,11 +3412,13 @@
   unwind_protect_bool (get_input_from_eval_string);
   unwind_protect_bool (input_from_eval_string_pending);
   unwind_protect_bool (input_from_command_line_file);
+  unwind_protect_bool (parser_end_of_input);
   unwind_protect_str (current_eval_string);
 
   get_input_from_eval_string = true;
   input_from_eval_string_pending = true;
   input_from_command_line_file = false;
+  parser_end_of_input = false;
   current_eval_string = s;
 
   unwind_protect_ptr (global_command);
@@ -3431,22 +3439,25 @@
 
       tree_statement_list *command = global_command;
 
-      if (parse_status == 0 && command)
+      if (parse_status == 0)
         {
-	  retval = command->eval (silent, nargout);
-
-	  delete command;
-
-	  command = 0;
-
-	  if (error_state
-	      || tree_return_command::returning
-	      || tree_break_command::breaking
-	      || tree_continue_command::continuing)
+	  if (command)
+	    {
+	      retval = command->eval (silent, nargout);
+
+	      delete command;
+
+	      command = 0;
+
+	      if (error_state
+		  || tree_return_command::returning
+		  || tree_break_command::breaking
+		  || tree_continue_command::continuing)
+		break;
+	    }
+	  else if (parser_end_of_input)
 	    break;
-	}
-      else
-	break;
+        }
     }
   while (parse_status == 0);
 
--- a/src/toplev.cc
+++ b/src/toplev.cc
@@ -125,48 +125,51 @@
 
       retval = yyparse ();
 
-      if (retval == 0 && global_command)
+      if (retval == 0)
 	{
-	  global_command->eval ();
-
-	  delete global_command;
-
-	  global_command = 0;
-
-	  if (! (interactive || forced_interactive))
+	  if (global_command)
 	    {
-	      bool quit = (tree_return_command::returning
-			   || tree_break_command::breaking);
+	      global_command->eval ();
 
-	      if (tree_return_command::returning)
-		tree_return_command::returning = 0;
-
-	      if (tree_break_command::breaking)
-		tree_break_command::breaking--;
+	      delete global_command;
 
-	      if (quit)
-		break;
-	    }
+	      global_command = 0;
 
-	  if (error_state)
-	    {
 	      if (! (interactive || forced_interactive))
 		{
-		  // We should exit with a non-zero status.
-		  retval = 1;
-		  break;
+		  bool quit = (tree_return_command::returning
+			       || tree_break_command::breaking);
+
+		  if (tree_return_command::returning)
+		    tree_return_command::returning = 0;
+
+		  if (tree_break_command::breaking)
+		    tree_break_command::breaking--;
+
+		  if (quit)
+		    break;
+		}
+
+	      if (error_state)
+		{
+		  if (! (interactive || forced_interactive))
+		    {
+		      // We should exit with a non-zero status.
+		      retval = 1;
+		      break;
+		    }
+		}
+	      else
+		{
+		  if (octave_completion_matches_called)
+		    octave_completion_matches_called = false;	    
+		  else
+		    command_editor::increment_current_command_number ();
 		}
 	    }
-	  else
-	    {
-	      if (octave_completion_matches_called)
-		octave_completion_matches_called = false;	    
-	      else
-		command_editor::increment_current_command_number ();
-	    }
+	  else if (parser_end_of_input)
+	    break;
 	}
-      else
-	break;
     }
   while (retval == 0);
 
--- a/test/octave.test/linalg/linalg.exp
+++ b/test/octave.test/linalg/linalg.exp
@@ -186,6 +186,10 @@
 set prog_output "\n... qr:.*"
 do_test qr-6.m
 
+set test qr-7
+set prog_output "ans = 1";
+do_test qr-7.m
+
 set test schur-1
 set prog_output "ans = 1"
 do_test schur-1.m