Mercurial > hg > octave-lojdl
changeset 1919:b582f7f765e0
[project @ 1996-02-11 01:58:47 by jwe]
author | jwe |
---|---|
date | Sun, 11 Feb 1996 02:02:55 +0000 |
parents | adde795927b5 |
children | 6fdb56b2cb33 |
files | liboctave/CmplxLU.cc liboctave/CmplxLU.h liboctave/dbleLU.cc liboctave/dbleLU.h |
diffstat | 4 files changed, 112 insertions(+), 68 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxLU.cc +++ b/liboctave/CmplxLU.cc @@ -43,8 +43,12 @@ ComplexLU::ComplexLU (const ComplexMatrix& a) { + ipvt = 0; + pvt = 0; + int a_nr = a.rows (); int a_nc = a.cols (); + if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) { (*current_liboctave_error_handler) ("ComplexLU requires square matrix"); @@ -53,51 +57,60 @@ int n = a_nr; - int *ipvt = new int [n]; - int *pvt = new int [n]; - Complex *tmp_data = dup (a.data (), a.length ()); + ipvt = new int [n]; + pvt = new int [n]; + + ComplexMatrix A_fact = a; + Complex *tmp_data = A_fact.fortran_vec (); + int info = 0; Complex *dummy = 0; - F77_FCN (zgesv, ZGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info); + F77_XFCN (zgesv, ZGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info); - ComplexMatrix A_fact (tmp_data, n, n); - - for (int i = 0; i < n; i++) + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgesv"); + else { - ipvt[i] -= 1; - pvt[i] = i; - } + for (int i = 0; i < n; i++) + { + ipvt[i] -= 1; + pvt[i] = i; + } - for (int i = 0; i < n - 1; i++) - { - int k = ipvt[i]; - if (k != i) + for (int i = 0; i < n - 1; i++) { - int tmp = pvt[k]; - pvt[k] = pvt[i]; - pvt[i] = tmp; + int k = ipvt[i]; + if (k != i) + { + int tmp = pvt[k]; + pvt[k] = pvt[i]; + pvt[i] = tmp; + } + } + + l.resize (n, n, 0.0); + u.resize (n, n, 0.0); + p.resize (n, n, 0.0); + + for (int i = 0; i < n; i++) + { + p.elem (i, pvt[i]) = 1.0; + + l.elem (i, i) = 1.0; + for (int j = 0; j < i; j++) + l.elem (i, j) = A_fact.elem (i, j); + + for (int j = i; j < n; j++) + u.elem (i, j) = A_fact.elem (i, j); } } - l.resize (n, n, 0.0); - u.resize (n, n, 0.0); - p.resize (n, n, 0.0); - - for (int i = 0; i < n; i++) - { - p.elem (i, pvt[i]) = 1.0; + delete [] ipvt; + ipvt = 0; - l.elem (i, i) = 1.0; - for (int j = 0; j < i; j++) - l.elem (i, j) = A_fact.elem (i, j); - - for (int j = i; j < n; j++) - u.elem (i, j) = A_fact.elem (i, j); - } - - delete [] ipvt; delete [] pvt; + pvt = 0; } /*
--- a/liboctave/CmplxLU.h +++ b/liboctave/CmplxLU.h @@ -38,7 +38,7 @@ { public: - ComplexLU (void) : l (), u (), p () { } + ComplexLU (void) : l (), u (), p (), ipvt (0), pvt (0) { } ComplexLU (const ComplexMatrix& a); @@ -55,6 +55,12 @@ return *this; } + ~ComplexLU (void) + { + delete [] ipvt; + delete [] pvt; + } + ComplexMatrix L (void) const { return l; } ComplexMatrix U (void) const { return u; } @@ -67,6 +73,9 @@ ComplexMatrix l; ComplexMatrix u; Matrix p; + + int *ipvt; + int *pvt; }; #endif
--- a/liboctave/dbleLU.cc +++ b/liboctave/dbleLU.cc @@ -43,8 +43,12 @@ LU::LU (const Matrix& a) { + ipvt = 0; + pvt = 0; + int a_nr = a.rows (); int a_nc = a.cols (); + if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) { (*current_liboctave_error_handler) ("LU requires square matrix"); @@ -53,51 +57,60 @@ int n = a_nr; - int *ipvt = new int [n]; - int *pvt = new int [n]; - double *tmp_data = dup (a.data (), a.length ()); + ipvt = new int [n]; + pvt = new int [n]; + + Matrix A_fact = a; + double *tmp_data = A_fact.fortran_vec (); + int info = 0; double dummy = 0; - F77_FCN (dgesv, DGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info); + F77_XFCN (dgesv, DGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info); - Matrix A_fact (tmp_data, n, n); - - for (int i = 0; i < n; i++) + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgesv"); + else { - ipvt[i] -= 1; - pvt[i] = i; - } + for (int i = 0; i < n; i++) + { + ipvt[i] -= 1; + pvt[i] = i; + } - for (int i = 0; i < n - 1; i++) - { - int k = ipvt[i]; - if (k != i) + for (int i = 0; i < n - 1; i++) { - int tmp = pvt[k]; - pvt[k] = pvt[i]; - pvt[i] = tmp; + int k = ipvt[i]; + if (k != i) + { + int tmp = pvt[k]; + pvt[k] = pvt[i]; + pvt[i] = tmp; + } + } + + l.resize (n, n, 0.0); + u.resize (n, n, 0.0); + p.resize (n, n, 0.0); + + for (int i = 0; i < n; i++) + { + p.elem (i, pvt[i]) = 1.0; + + l.elem (i, i) = 1.0; + for (int j = 0; j < i; j++) + l.elem (i, j) = A_fact.elem (i, j); + + for (int j = i; j < n; j++) + u.elem (i, j) = A_fact.elem (i, j); } } - l.resize (n, n, 0.0); - u.resize (n, n, 0.0); - p.resize (n, n, 0.0); - - for (int i = 0; i < n; i++) - { - p.elem (i, pvt[i]) = 1.0; + delete [] ipvt; + ipvt = 0; - l.elem (i, i) = 1.0; - for (int j = 0; j < i; j++) - l.elem (i, j) = A_fact.elem (i, j); - - for (int j = i; j < n; j++) - u.elem (i, j) = A_fact.elem (i, j); - } - - delete [] ipvt; delete [] pvt; + pvt = 0; } /*
--- a/liboctave/dbleLU.h +++ b/liboctave/dbleLU.h @@ -36,7 +36,7 @@ { public: - LU (void) : l (), u (), p () { } + LU (void) : l (), u (), p (), ipvt (0), pvt (0) { } LU (const Matrix& a); @@ -53,6 +53,12 @@ return *this; } + ~LU (void) + { + delete [] ipvt; + delete [] pvt; + } + Matrix L (void) const { return l; } Matrix U (void) const { return u; } @@ -66,6 +72,9 @@ Matrix l; Matrix u; Matrix p; + + int *ipvt; + int *pvt; }; #endif