Mercurial > hg > octave-lojdl
changeset 2763:d9d00d7e271e
[project @ 1997-03-01 02:14:33 by jwe]
author | jwe |
---|---|
date | Sat, 01 Mar 1997 02:14:33 +0000 |
parents | ac1427f5c9e6 |
children | 2c0f259cf83d |
files | liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/CmplxQRP.cc liboctave/CmplxQRP.h liboctave/dbleQR.cc liboctave/dbleQR.h liboctave/dbleQRP.cc liboctave/dbleQRP.h |
diffstat | 8 files changed, 73 insertions(+), 60 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxQR.cc +++ b/liboctave/CmplxQR.cc @@ -45,6 +45,13 @@ } ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type) + : q (), r () +{ + init (a, qr_type); +} + +void +ComplexQR::init (const ComplexMatrix& a, QR::type qr_type) { int m = a.rows (); int n = a.cols ();
--- a/liboctave/CmplxQR.h +++ b/liboctave/CmplxQR.h @@ -39,7 +39,7 @@ ComplexQR (void) : q (), r () { } - ComplexQR (const ComplexMatrix& A, QR::type qr_type = QR::std); + ComplexQR (const ComplexMatrix&, QR::type = QR::std); ComplexQR (const ComplexQR& a) : q (a.q), r (a.r) { } @@ -55,10 +55,13 @@ ~ComplexQR (void) { } + void init (const ComplexMatrix&, QR::type = QR::std); + ComplexMatrix Q (void) const { return q; } + ComplexMatrix R (void) const { return r; } - friend ostream& operator << (ostream& os, const ComplexQR& a); + friend ostream& operator << (ostream&, const ComplexQR&); protected:
--- a/liboctave/CmplxQRP.cc +++ b/liboctave/CmplxQRP.cc @@ -49,6 +49,13 @@ // It would be best to share some of this code with ComplexQR class... ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type) + : ComplexQR (), p () +{ + init (a, qr_type); +} + +void +ComplexQRP::init (const ComplexMatrix& a, QR::type qr_type) { assert (qr_type != QR::raw); @@ -62,24 +69,18 @@ return; } - int min_mn = m < n ? m : n; - Array<Complex> tau (min_mn); + Array<Complex> tau (m < n ? m : n); Complex *ptau = tau.fortran_vec (); - int lwork = n; + int lwork = 3*n > 32*m ? 3*n : 32*m; Array<Complex> work (lwork); Complex *pwork = work.fortran_vec (); int info = 0; - ComplexMatrix A_fact; + ComplexMatrix A_fact = a; if (m > n) - { - A_fact.resize (m, m); - A_fact.insert (a, 0, 0); - } - else - A_fact = a; + A_fact.resize (m, m, 0.0); Complex *tmp_data = A_fact.fortran_vec (); @@ -101,7 +102,7 @@ // Form Permutation matrix (if economy is requested, return the // indices only!) - if (qr_type == QR::economy && m > n) + if (qr_type == QR::economy) { p.resize (1, n, 0.0); for (int j = 0; j < n; j++) @@ -114,18 +115,12 @@ p.elem (jpvt.elem (j) - 1, j) = 1.0; } - volatile int n2; + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); - if (qr_type == QR::economy && m > n) - { - n2 = n; - r.resize (n, n, 0.0); - } - else - { - n2 = m; - r.resize (m, n, 0.0); - } + int min_mn = m < n ? m : n; for (int j = 0; j < n; j++) { @@ -134,11 +129,11 @@ r.elem (i, j) = A_fact.elem (i, j); } - lwork = 32*m; - work.resize (lwork); - Complex *pwork = work.fortran_vec (); + int n2 = m; + if (qr_type == QR::economy) + n2 = min_mn; - 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.h +++ b/liboctave/CmplxQRP.h @@ -38,7 +38,7 @@ ComplexQRP (void) : ComplexQR (), p () { } - ComplexQRP (const ComplexMatrix& A, QR::type qr_type = QR::std); + ComplexQRP (const ComplexMatrix&, QR::type = QR::std); ComplexQRP (const ComplexQRP& a) : ComplexQR (a), p (a.p) { } @@ -54,9 +54,11 @@ ~ComplexQRP (void) { } + void init (const ComplexMatrix&, QR::type = QR::std); + Matrix P (void) const { return p; } - friend ostream& operator << (ostream& os, const ComplexQRP& a); + friend ostream& operator << (ostream&, const ComplexQRP&); private:
--- a/liboctave/dbleQR.cc +++ b/liboctave/dbleQR.cc @@ -45,6 +45,13 @@ } QR::QR (const Matrix& a, QR::type qr_type) + : q (), r () +{ + init (a, qr_type); +} + +void +QR::init (const Matrix& a, QR::type qr_type) { int m = a.rows (); int n = a.cols ();
--- a/liboctave/dbleQR.h +++ b/liboctave/dbleQR.h @@ -45,7 +45,7 @@ QR (void) : q (), r () { } - QR (const Matrix& A, type qr_type = QR::std); + QR (const Matrix&, QR::type = QR::std); QR (const QR& a) : q (a.q), r (a.r) { } @@ -61,11 +61,13 @@ ~QR (void) { } + void init (const Matrix&, QR::type); + Matrix Q (void) const { return q; } Matrix R (void) const { return r; } - friend ostream& operator << (ostream& os, const QR& a); + friend ostream& operator << (ostream&, const QR&); protected:
--- a/liboctave/dbleQRP.cc +++ b/liboctave/dbleQRP.cc @@ -49,6 +49,13 @@ // It would be best to share some of this code with QR class... QRP::QRP (const Matrix& a, QR::type qr_type) + : QR (), p () +{ + init (a, qr_type); +} + +void +QRP::init (const Matrix& a, QR::type qr_type) { assert (qr_type != QR::raw); @@ -61,24 +68,18 @@ return; } - int min_mn = m < n ? m : n; - Array<double> tau (min_mn); + Array<double> tau (m < n ? m : n); double *ptau = tau.fortran_vec (); - int lwork = 3*n; + int lwork = 3*n > 32*m ? 3*n : 32*m; Array<double> work (lwork); double *pwork = work.fortran_vec (); int info = 0; - Matrix A_fact; + Matrix A_fact = a; if (m > n) - { - A_fact.resize (m, m); - A_fact.insert (a, 0, 0); - } - else - A_fact = a; + A_fact.resize (m, m, 0.0); double *tmp_data = A_fact.fortran_vec (); @@ -96,7 +97,7 @@ // Form Permutation matrix (if economy is requested, return the // indices only!) - if (qr_type == QR::economy && m > n) + if (qr_type == QR::economy) { p.resize (1, n, 0.0); for (int j = 0; j < n; j++) @@ -109,18 +110,12 @@ p.elem (jpvt.elem (j) - 1, j) = 1.0; } - volatile int n2; + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); - if (qr_type == QR::economy && m > n) - { - n2 = n; - r.resize (n, n, 0.0); - } - else - { - n2 = m; - r.resize (m, n, 0.0); - } + int min_mn = m < n ? m : n; for (int j = 0; j < n; j++) { @@ -129,11 +124,11 @@ r.elem (i, j) = A_fact.elem (i, j); } - lwork = 32*m; - work.resize (lwork); - double *pwork = work.fortran_vec (); + int n2 = m; + if (qr_type == QR::economy) + n2 = min_mn; - 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.h +++ b/liboctave/dbleQRP.h @@ -38,7 +38,7 @@ QRP (void) : QR (), p () { } - QRP (const Matrix& A, QR::type qr_type = QR::std); + QRP (const Matrix&, QR::type = QR::std); QRP (const QRP& a) : QR (a), p (a.p) { } @@ -55,9 +55,11 @@ ~QRP (void) { } + void init (const Matrix&, QR::type = QR::std); + Matrix P (void) const { return p; } - friend ostream& operator << (ostream& os, const QRP& a); + friend ostream& operator << (ostream&, const QRP&); protected: