Mercurial > hg > octave-lojdl
changeset 1922:8a57554f3142
[project @ 1996-02-11 02:53:20 by jwe]
author | jwe |
---|---|
date | Sun, 11 Feb 1996 02:53:20 +0000 |
parents | ce0db0ca0729 |
children | 1d4f67de06f4 |
files | liboctave/CmplxQRP.cc liboctave/CmplxQRP.h liboctave/dbleQRP.cc liboctave/dbleQRP.h |
diffstat | 4 files changed, 176 insertions(+), 101 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxQRP.cc +++ b/liboctave/CmplxQRP.cc @@ -53,6 +53,11 @@ { assert (qr_type != QR::raw); + tau = 0; + work = 0; + tmp_data = 0; + jpvt = 0; + int m = a.rows (); int n = a.cols (); @@ -63,11 +68,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 = n; - Complex *work = new Complex[lwork]; + work = new Complex[lwork]; + int info = 0; if (m > n) @@ -78,10 +84,9 @@ else tmp_data = dup (a.data (), a.length ()); - work = new Complex[n]; - double *rwork = new double[2*n]; - int *jpvt = new int[n]; + rwork = new double[2*n]; + jpvt = new int[n]; // Clear Pivot vector (code to enforce a certain permutation would // go here...) @@ -89,58 +94,85 @@ for (int i = 0; i < n; i++) jpvt[i] = 0; - F77_FCN (zgeqpf, ZGEQPF) (m, n, tmp_data, m, jpvt, tau, work, rwork, - info); + F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, jpvt, tau, work, rwork, + info)); - // Form Permutation matrix (if economy is requested, return the - // indices only!) + delete [] work; + work = 0; - if (qr_type == QR::economy && m > n) - { - p.resize (1, n, 0.0); - for (int j = 0; j < n; j++) - p.elem (0, j) = jpvt[j]; - } + delete [] rwork; + rwork = 0; + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf"); else { - p.resize (n, n, 0.0); + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy && m > n) + { + p.resize (1, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (0, j) = jpvt[j]; + } + else + { + p.resize (n, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (jpvt[j]-1, j) = 1.0; + } + + delete [] jpvt; + jpvt = 0; + + volatile int n2; + + 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++) - p.elem (jpvt[j]-1, j) = 1.0; + { + 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]; + + F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, tau, work, + lwork, info)); + + 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; + tau = 0; + + delete [] work; + work = 0; + + delete [] tmp_data; + tmp_data = 0; + delete [] jpvt; - delete [] work; - - int n2; - 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]; - - F77_FCN (zungqr, ZUNGQR) (m, m, min_mn, tmp_data, m, tau, work, - lwork, info); - - q = ComplexMatrix (tmp_data, m, m); - q.resize (m, n2); - - delete [] tau; - delete [] work; + jpvt = 0; } /*
--- a/liboctave/CmplxQRP.h +++ b/liboctave/CmplxQRP.h @@ -37,7 +37,7 @@ { public: - ComplexQRP (void) : ComplexQR (), p () { } + ComplexQRP (void) : ComplexQR (), p (), jpvt (0), rwork (0) { } ComplexQRP (const ComplexMatrix& A, QR::type qr_type = QR::std); @@ -53,6 +53,12 @@ return *this; } + ~ComplexQRP (void) + { + delete [] jpvt; + delete [] rwork; + } + Matrix P (void) const { return p; } friend ostream& operator << (ostream& os, const ComplexQRP& a); @@ -60,6 +66,9 @@ private: Matrix p; + + int *jpvt; + double *rwork; }; #endif
--- a/liboctave/dbleQRP.cc +++ b/liboctave/dbleQRP.cc @@ -53,6 +53,11 @@ { assert (qr_type != QR::raw); + tau = 0; + work = 0; + tmp_data = 0; + jpvt = 0; + int m = a.rows (); int n = a.cols (); @@ -62,11 +67,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 = 3*n; - double *work = new double[lwork]; + work = new double[lwork]; + int info = 0; if (m > n) @@ -77,7 +83,7 @@ else tmp_data = dup (a.data (), a.length ()); - int *jpvt = new int[n]; + jpvt = new int[n]; // Clear Pivot vector (code to enforce a certain permutation would // go here...) @@ -85,57 +91,81 @@ for (int i = 0; i < n; i++) jpvt[i] = 0; - F77_FCN (dgeqpf, DGEQPF) (m, n, tmp_data, m, jpvt, tau, work, info); - - // Form Permutation matrix (if economy is requested, return the - // indices only!) + F77_XFCN (dgeqpf, DGEQPF, (m, n, tmp_data, m, jpvt, tau, work, info)); - if (qr_type == QR::economy && m > n) - { - p.resize (1, n, 0.0); - for (int j = 0; j < n; j++) - p.elem (0, j) = jpvt[j]; - } + delete [] work; + work = 0; + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgeqpf"); else { - p.resize (n, n, 0.0); + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy && m > n) + { + p.resize (1, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (0, j) = jpvt[j]; + } + else + { + p.resize (n, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (jpvt[j]-1, j) = 1.0; + } + + delete [] jpvt; + jpvt = 0; + + volatile int n2; + + 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++) - p.elem (jpvt[j]-1, j) = 1.0; + { + 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]; + + F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, tau, work, + lwork, info)); + + 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; + tau = 0; + + delete [] work; + work = 0; + + delete [] tmp_data; + tmp_data = 0; + delete [] jpvt; - delete [] work; - - int n2; - 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]; - - F77_FCN (dorgqr, DORGQR) (m, m, min_mn, tmp_data, m, tau, work, - lwork, info); - - q = Matrix (tmp_data, m, m); - q.resize (m, n2); - - delete [] tau; - delete [] work; + jpvt = 0; } /*
--- a/liboctave/dbleQRP.h +++ b/liboctave/dbleQRP.h @@ -37,7 +37,7 @@ { public: - QRP (void) : QR (), p () { } + QRP (void) : QR (), p (), jpvt (0) { } QRP (const Matrix& A, QR::type qr_type = QR::std); @@ -54,13 +54,17 @@ return *this; } + ~QRP (void) { delete [] jpvt; } + Matrix P (void) const { return p; } friend ostream& operator << (ostream& os, const QRP& a); -private: +protected: Matrix p; + + int *jpvt; }; #endif