Mercurial > octave-nkf

changeset 539:5ec10a984241

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-07-21 22:40:04 by jwe]
author jwe
date Thu, 21 Jul 1994 22:42:09 +0000
parents 8e134d3b21c9
children c07674bbc3b1
files liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleQR.cc liboctave/dbleQR.h src/qr.cc src/svd.cc
diffstat 6 files changed, 167 insertions(+), 53 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxQR.cc	Thu Jul 21 22:40:04 1994 +0000
+++ b/liboctave/CmplxQR.cc	Thu Jul 21 22:42:09 1994 +0000
@@ -43,7 +43,7 @@
 			const int*, Complex*, Complex*, const int*, int*);
 }
 
-ComplexQR::ComplexQR (const ComplexMatrix& a)
+ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type)
 {
   int m = a.rows ();
   int n = a.cols ();
@@ -74,23 +74,47 @@
 
   delete [] work;
 
-  r.resize (m, n, 0.0);
-  for (int j = 0; j < n; j++)
+  if (qr_type == QR::raw)
+    {
+      for (int j = 0; j < min_mn; j++)
+	{
+	  int limit = j < min_mn - 1 ? j : min_mn - 1;
+	  for (int i = limit + 1; i < m; i++)
+	    tmp_data[m*j+i] *= tau[j];
+	}
+    }
+  else
     {
-      int limit = j < min_mn-1 ? j : min_mn-1;
-      for (int i = 0; i <= limit; i++)
-	r.elem (i, j) = tmp_data[m*j+i];
-    }
+      int m2;
+      if (qr_type == QR::economy && m > n)
+	{
+	  m2 = n;
+	  r.resize (n, n, 0.0);
+	}
+      else
+	{
+	  m2 = m;
+	  r.resize (m, n, 0.0);
+	}
 
-  lwork = 32*m;
-  work = new Complex[lwork];
-
-  F77_FCN (zungqr) (&m, &m, &min_mn, tmp_data, &m, tau, work, &lwork, &info);
+      for (int j = 0; j < n; j++)
+	{
+	  int limit = j < min_mn-1 ? j : min_mn-1;
+	  for (int i = 0; i <= limit; i++)
+	    r.elem (i, j) = tmp_data[m*j+i];
+	}
 
-  q = ComplexMatrix (tmp_data, m, m);
+      lwork = 32*m;
+      work = new Complex[lwork];
+
+      F77_FCN (zungqr) (&m, &m, &min_mn, tmp_data, &m, tau, work,
+			&lwork, &info);
 
-  delete [] tau;
-  delete [] work;
+      q = ComplexMatrix (tmp_data, m, m);
+
+      delete [] tau;
+      delete [] work;
+    }
 }
 
 /*
--- a/liboctave/CmplxQR.h	Thu Jul 21 22:40:04 1994 +0000
+++ b/liboctave/CmplxQR.h	Thu Jul 21 22:42:09 1994 +0000
@@ -31,6 +31,7 @@
 class ostream;
 
 #include "CMatrix.h"
+#include "dbleQR.h"
 
 extern "C++" {
 
@@ -40,7 +41,7 @@
 
   ComplexQR (void) {}
 
-  ComplexQR (const ComplexMatrix& A);
+  ComplexQR (const ComplexMatrix& A, QR::type qr_type = QR::std);
 
   ComplexQR (const ComplexQR& a);
 
@@ -51,7 +52,7 @@
 
   friend ostream&  operator << (ostream& os, const ComplexQR& a);
 
-private:
+protected:
 
   ComplexMatrix q;
   ComplexMatrix r;
--- a/liboctave/dbleQR.cc	Thu Jul 21 22:40:04 1994 +0000
+++ b/liboctave/dbleQR.cc	Thu Jul 21 22:42:09 1994 +0000
@@ -41,9 +41,12 @@
 
   int F77_FCN (dorgqr) (const int*, const int*, const int*, double*,
 			const int*, double*, double*, const int*, int*);
+
+  int F77_FCN (dgeqpf) (const int*, const int*, double*, const int*,
+			int*, double*, double*, int*);
 }
 
-QR::QR (const Matrix& a)
+QR::QR (const Matrix& a, QR::type qr_type)
 {
   int m = a.rows ();
   int n = a.cols ();
@@ -73,23 +76,47 @@
 
   delete [] work;
 
-  r.resize (m, n, 0.0);
-  for (int j = 0; j < n; j++)
+  if (qr_type == QR::raw)
+    {
+      for (int j = 0; j < min_mn; j++)
+	{
+	  int limit = j < min_mn - 1 ? j : min_mn - 1;
+	  for (int i = limit + 1; i < m; i++)
+	    tmp_data[m*j+i] *= tau[j];
+	}
+    }
+  else
     {
-      int limit = j < min_mn-1 ? j : min_mn-1;
-      for (int i = 0; i <= limit; i++)
-	r.elem (i, j) = tmp_data[m*j+i];
-    }
+      int m2;
+      if (qr_type == QR::economy && m > n)
+	{
+	  m2 = n;
+	  r.resize (n, n, 0.0);
+	}
+      else
+	{
+	  m2 = m;
+	  r.resize (m, n, 0.0);
+	}
 
-  lwork = 32*m;
-  work = new double[lwork];
-
-  F77_FCN (dorgqr) (&m, &m, &min_mn, tmp_data, &m, tau, work, &lwork, &info);
+      for (int j = 0; j < n; j++)
+	{
+	  int limit = j < min_mn-1 ? j : min_mn-1;
+	  for (int i = 0; i <= limit; i++)
+	    r.elem (i, j) = tmp_data[m*j+i];
+	}
 
-  q = Matrix (tmp_data, m, m);
+      lwork = 32*m;
+      work = new double[lwork];
+
+      F77_FCN (dorgqr) (&m, &m, &min_mn, tmp_data, &m, tau, work,
+			&lwork, &info);
 
-  delete [] tau;
-  delete [] work;
+      q = Matrix (tmp_data, m, m);
+
+      delete [] tau;
+      delete [] work;
+    }
 }
 
 /*
--- a/liboctave/dbleQR.h	Thu Jul 21 22:40:04 1994 +0000
+++ b/liboctave/dbleQR.h	Thu Jul 21 22:42:09 1994 +0000
@@ -38,9 +38,16 @@
 {
 public:
 
+  enum type
+    {
+      std,
+      raw,
+      economy,
+    };
+
   QR (void) {}
 
-  QR (const Matrix& A);
+  QR (const Matrix& A, type qr_type = QR::std);
 
   QR (const QR& a);
 
@@ -51,7 +58,7 @@
 
   friend ostream&  operator << (ostream& os, const QR& a);
 
-private:
+protected:
 
   Matrix q;
   Matrix r;
--- a/src/qr.cc	Thu Jul 21 22:40:04 1994 +0000
+++ b/src/qr.cc	Thu Jul 21 22:42:09 1994 +0000
@@ -28,19 +28,37 @@
 #include "dbleQR.h"
 #include "CmplxQR.h"
 
+#include "dbleQRP.h"
+#include "CmplxQRP.h"
+
 #include "tree-const.h"
 #include "user-prefs.h"
 #include "gripes.h"
 #include "defun-dld.h"
 
 DEFUN_DLD ("qr", Fqr, Sqr, 2, 2,
-  "[Q, R] = qr (X): form QR factorization of X")
+  "[Q, R] = qr (X):      form Q unitary and R upper triangular such\n\
+                       that Q * R = X\n\
+\n\
+[Q, R] = qr (X, 0):    form the economy decomposition such that if X is\n\
+                       if X is m by n then only the first n columns of Q\n\
+                       are computed.\n\
+\n\
+[Q, R, P] = qr (X):    form QRP factorization of X where\n\
+                       P is a permutation matrix such that\n\
+                       A * P = Q * R\n\
+\n\
+[Q, R, P] = qr (X, 0): form the economy decomposition with \n\
+                       permutation vector P such that Q * R = X (:, P)\n\
+\n\
+qr (X) alone returns the output of the LAPACK routine dgeqrf, such\n\
+that R = triu (qr (X))")
 {
   Octave_object retval;
 
   int nargin = args.length ();
 
-  if (nargin != 2 || nargout > 2)
+  if (nargin != 2 && nargin != 3 || nargout > 3)
     {
       print_usage ("qr");
       return retval;
@@ -48,10 +66,7 @@
 
   tree_constant tmp = args(1).make_numeric ();
     
-  int nr = tmp.rows ();
-  int nc = tmp.columns ();
-
-  if (nr == 0 || nc == 0)
+  if (tmp.rows () == 0 || tmp.columns () == 0)
     {
       int flag = user_pref.propagate_empty_matrices;
       if (flag != 0)
@@ -59,6 +74,7 @@
 	  if (flag < 0)
 	    gripe_empty_arg ("qr", 0);
 	  Matrix m;
+	  retval(2) = m;
 	  retval(1) = m;
 	  retval(0) = m;
 	}
@@ -68,36 +84,71 @@
       return retval;
     }
 
+  QR::type type = nargout == 1 ? QR::raw
+    : (nargin == 3 ? QR::economy : QR::std);
+
   switch (tmp.const_type ())
     {
     case tree_constant_rep::matrix_constant:
       {
 	Matrix m = tmp.matrix_value ();
-	QR fact (m);
-	retval(1) = fact.R ();
-	retval(0) = fact.Q ();
+	if (nargout < 3)
+	  {
+	    QR fact (m, type);
+	    retval(1) = fact.R ();
+	    retval(0) = fact.Q ();
+	  }
+	else
+	  {
+	    QRP fact (m, type);
+	    retval(2) = fact.P ();
+	    retval(1) = fact.R ();
+	    retval(0) = fact.Q ();
+	  }
       }
       break;
     case tree_constant_rep::complex_matrix_constant:
       {
 	ComplexMatrix m = tmp.complex_matrix_value ();
-	ComplexQR fact (m);
-	retval(1) = fact.R ();
-	retval(0) = fact.Q ();
+	if (nargout < 3)
+	  {
+	    ComplexQR fact (m, type);
+	    retval(1) = fact.R ();
+	    retval(0) = fact.Q ();
+	  }
+	else
+	  {
+	    ComplexQRP fact (m, type);
+	    retval(2) = fact.P ();
+	    retval(1) = fact.R ();
+	    retval(0) = fact.Q ();
+	  }
       }
       break;
     case tree_constant_rep::scalar_constant:
       {
 	double d = tmp.double_value ();
-	retval(1) = d;
-	retval(0) = 1.0;
+	if (nargout == 1)
+	  retval(0) = d;
+	else
+	  {
+	    retval(2) = 1.0;
+	    retval(1) = d;
+	    retval(0) = 1.0;
+	  }
       }
       break;
     case tree_constant_rep::complex_scalar_constant:
       {
 	Complex c = tmp.complex_value ();
-	retval(1) = c;
-	retval(0) = 1.0;
+	if (nargout == 1)
+	  retval(0) = c;
+	else
+	  {
+	    retval(2) = 1.0;
+	    retval(1) = c;
+	    retval(0) = 1.0;
+	  }
       }
       break;
     default:
--- a/src/svd.cc	Thu Jul 21 22:40:04 1994 +0000
+++ b/src/svd.cc	Thu Jul 21 22:42:09 1994 +0000
@@ -35,15 +35,17 @@
 #include "defun-dld.h"
 
 DEFUN_DLD ("svd", Fsvd, Ssvd, 2, 3,
-  "S = svd (X) or [U, S, V] = svd (X)\n\
+  "S = svd (X) or [U, S, V] = svd (X [, 0])\n\
 \n\
-compute the singular value decomposition of X")
+Compute the singular value decomposition of X.  Given a second input\n\
+argument, an `economy' sized factorization is computed that omits\n\
+unnecessary rows and columns of U and V")
 {
   Octave_object retval;
 
   int nargin = args.length ();
 
-  if (nargin != 2 || nargout == 2 || nargout > 3)
+  if (nargin < 2 || nargin > 3 || nargout == 2 || nargout > 3)
     {
       print_usage ("svd");
       return retval;
@@ -94,12 +96,14 @@
       break;
     }
 
+  SVD::type type = (nargin == 3) ? SVD::economy : SVD::std;
+
   switch (arg.const_type ())
     {
     case tree_constant_rep::scalar_constant:
     case tree_constant_rep::matrix_constant:
       {
-	SVD result (tmp);
+	SVD result (tmp, type);
 
 	DiagMatrix sigma = result.singular_values ();
 
@@ -120,7 +124,7 @@
     case tree_constant_rep::complex_scalar_constant:
     case tree_constant_rep::complex_matrix_constant:
       {
-	ComplexSVD result (ctmp);
+	ComplexSVD result (ctmp, type);
 
 	DiagMatrix sigma = result.singular_values ();