Mercurial > hg > octave-lojdl
changeset 1921:ce0db0ca0729
[project @ 1996-02-11 02:31:58 by jwe]
author | jwe |
---|---|
date | Sun, 11 Feb 1996 02:31:58 +0000 |
parents | 6fdb56b2cb33 |
children | 8a57554f3142 |
files | liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleQR.cc liboctave/dbleQR.h |
diffstat | 4 files changed, 144 insertions(+), 73 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxQR.cc +++ b/liboctave/CmplxQR.cc @@ -47,6 +47,10 @@ ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type) { + tau = 0; + work = 0; + tmp_data = 0; + int m = a.rows (); int n = a.cols (); @@ -57,11 +61,12 @@ return; } - Complex *tmp_data; int min_mn = m < n ? m : n; - Complex *tau = new Complex[min_mn]; + tau = new Complex[min_mn]; + int lwork = 32*n; - Complex *work = new Complex[lwork]; + work = new Complex[lwork]; + int info = 0; if (m > n) @@ -72,53 +77,72 @@ else tmp_data = dup (a.data (), a.length ()); - F77_FCN (zgeqrf, ZGEQRF) (m, n, tmp_data, m, tau, work, lwork, - info); + F77_XFCN (zgeqrf, ZGEQRF, (m, n, tmp_data, m, tau, work, lwork, info)); delete [] work; + work = 0; - 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]; - } - } + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgeqrf"); else { - int n2; - if (qr_type == QR::economy && m > n) + if (qr_type == QR::raw) { - n2 = n; - r.resize (n, n, 0.0); + 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 { - n2 = m; - r.resize (m, n, 0.0); - } + volatile int n2; - 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]; - } + if (qr_type == QR::economy && m > n) + { + n2 = n; + r.resize (n, n, 0.0); + } + else + { + n2 = m; + r.resize (m, n, 0.0); + } + + 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]; + } + + lwork = 32*m; + work = new Complex[lwork]; - lwork = 32*m; - work = new Complex[lwork]; + F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, tau, work, + lwork, info)); - F77_FCN (zungqr, ZUNGQR) (m, m, min_mn, tmp_data, m, tau, work, - lwork, info); + delete [] work; + work = 0; - q = ComplexMatrix (tmp_data, m, m); - q.resize (m, n2); + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zungqr"); + else + { + q = ComplexMatrix (tmp_data, m, m); + tmp_data = 0; + q.resize (m, n2); + } + } + } - delete [] tau; - delete [] work; - } + delete [] tau; + tau = 0; + + delete [] tmp_data; + tmp_data = 0; } /*
--- a/liboctave/CmplxQR.h +++ b/liboctave/CmplxQR.h @@ -38,7 +38,7 @@ { public: - ComplexQR (void) : q (), r () { } + ComplexQR (void) : q (), r (), tau (0), work (0), tmp_data (0) { } ComplexQR (const ComplexMatrix& A, QR::type qr_type = QR::std); @@ -54,6 +54,13 @@ return *this; } + ~ComplexQR (void) + { + delete [] tau; + delete [] work; + delete [] tmp_data; + } + ComplexMatrix Q (void) const { return q; } ComplexMatrix R (void) const { return r; } @@ -63,6 +70,10 @@ ComplexMatrix q; ComplexMatrix r; + + Complex *tau; + Complex *work; + Complex *tmp_data; }; #endif
--- a/liboctave/dbleQR.cc +++ b/liboctave/dbleQR.cc @@ -47,6 +47,10 @@ QR::QR (const Matrix& a, QR::type qr_type) { + tau = 0; + work = 0; + tmp_data = 0; + int m = a.rows (); int n = a.cols (); @@ -56,11 +60,12 @@ return; } - double *tmp_data; int min_mn = m < n ? m : n; - double *tau = new double[min_mn]; + tau = new double[min_mn]; + int lwork = 32*n; - double *work = new double[lwork]; + work = new double[lwork]; + int info = 0; if (m > n) @@ -71,52 +76,72 @@ else tmp_data = dup (a.data (), a.length ()); - F77_FCN (dgeqrf, DGEQRF) (m, n, tmp_data, m, tau, work, lwork, info); + F77_XFCN (dgeqrf, DGEQRF, (m, n, tmp_data, m, tau, work, lwork, info)); delete [] work; + work = 0; - 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]; - } - } + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgeqrf"); else { - int n2; - if (qr_type == QR::economy && m > n) + if (qr_type == QR::raw) { - n2 = n; - r.resize (n, n, 0.0); + 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 { - n2 = m; - r.resize (m, n, 0.0); - } + volatile int n2; - 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]; - } + if (qr_type == QR::economy && m > n) + { + n2 = n; + r.resize (n, n, 0.0); + } + else + { + n2 = m; + r.resize (m, n, 0.0); + } + + 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]; + } + + lwork = 32*m; + work = new double[lwork]; - lwork = 32*m; - work = new double[lwork]; + F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, tau, work, + lwork, info)); - F77_FCN (dorgqr, DORGQR) (m, m, min_mn, tmp_data, m, tau, work, - lwork, info); + delete [] work; + work = 0; - q = Matrix (tmp_data, m, m); - q.resize (m, n2); + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dorgqr"); + else + { + q = Matrix (tmp_data, m, m); + tmp_data = 0; + q.resize (m, n2); + } + } + } - delete [] tau; - delete [] work; - } + delete [] tau; + tau = 0; + + delete [] tmp_data; + tmp_data = 0; } /*
--- a/liboctave/dbleQR.h +++ b/liboctave/dbleQR.h @@ -44,7 +44,7 @@ economy, }; - QR (void) : q (), r () { } + QR (void) : q (), r (), tau (0), work (0), tmp_data (0) { } QR (const Matrix& A, type qr_type = QR::std); @@ -60,6 +60,13 @@ return *this; } + ~QR (void) + { + delete [] tau; + delete [] work; + delete [] tmp_data; + } + Matrix Q (void) const { return q; } Matrix R (void) const { return r; } @@ -70,6 +77,10 @@ Matrix q; Matrix r; + + double *tau; + double *work; + double *tmp_data; }; #endif