# HG changeset patch # User jwe # Date 774830529 0 # Node ID 5ec10a984241bc21734f6307494160c81ef41251 # Parent 8e134d3b21c92f69c23005212024aaca61ea9ae1 [project @ 1994-07-21 22:40:04 by jwe] diff -r 8e134d3b21c9 -r 5ec10a984241 liboctave/CmplxQR.cc --- 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; + } } /* diff -r 8e134d3b21c9 -r 5ec10a984241 liboctave/CmplxQR.h --- 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; diff -r 8e134d3b21c9 -r 5ec10a984241 liboctave/dbleQR.cc --- 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; + } } /* diff -r 8e134d3b21c9 -r 5ec10a984241 liboctave/dbleQR.h --- 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; diff -r 8e134d3b21c9 -r 5ec10a984241 src/qr.cc --- 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: diff -r 8e134d3b21c9 -r 5ec10a984241 src/svd.cc --- 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 ();