Mercurial > hg > octave-nkf
diff liboctave/EIG.cc @ 1934:0e591d443ff0
[project @ 1996-02-12 03:25:10 by jwe]
author | jwe |
---|---|
date | Mon, 12 Feb 1996 03:25:10 +0000 |
parents | 1281a23a34dd |
children | 1b57120c997b |
line wrap: on
line diff
--- a/liboctave/EIG.cc +++ b/liboctave/EIG.cc @@ -52,106 +52,117 @@ int EIG::init (const Matrix& a) { - int a_nr = a.rows (); - if (a_nr != a.cols ()) + int n = a.rows (); + + if (n != a.cols ()) { (*current_liboctave_error_handler) ("EIG requires square matrix"); return -1; } - int n = a_nr; - int info; - double *tmp_data = dup (a.data (), a.length ()); - double *wr = new double[n]; - double *wi = new double[n]; + Matrix atmp = a; + double *tmp_data = atmp.fortran_vec (); + + Array<double> wr (n); + double *pwr = wr.fortran_vec (); + + Array<double> wi (n); + double *pwi = wi.fortran_vec (); + Matrix vr (n, n); double *pvr = vr.fortran_vec (); + int lwork = 8*n; - double *work = new double[lwork]; + + Array<double> work (lwork); + double *pwork = work.fortran_vec (); double *dummy = 0; int idummy = 1; - F77_FCN (dgeev, DGEEV) ("N", "V", n, tmp_data, n, wr, wi, dummy, - idummy, pvr, n, work, lwork, info, 1L, 1L); + F77_XFCN (dgeev, DGEEV, ("N", "V", n, tmp_data, n, pwr, pwi, dummy, + idummy, pvr, n, pwork, lwork, info, 1L, 1L)); - lambda.resize (n); - v.resize (n, n); - - for (int j = 0; j < n; j++) + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); + else { - if (wi[j] == 0.0) - { - lambda.elem (j) = Complex (wr[j]); - for (int i = 0; i < n; i++) - v.elem (i, j) = vr.elem (i, j); - } - else + lambda.resize (n); + v.resize (n, n); + + for (int j = 0; j < n; j++) { - if (j+1 >= n) + if (wi.elem (j) == 0.0) { - (*current_liboctave_error_handler) ("EIG: internal error"); - return -1; + lambda.elem (j) = Complex (wr.elem (j)); + for (int i = 0; i < n; i++) + v.elem (i, j) = vr.elem (i, j); } - - for (int i = 0; i < n; i++) + else { - lambda.elem (j) = Complex (wr[j], wi[j]); - lambda.elem (j+1) = Complex (wr[j+1], wi[j+1]); - double real_part = vr.elem (i, j); - double imag_part = vr.elem (i, j+1); - v.elem (i, j) = Complex (real_part, imag_part); - v.elem (i, j+1) = Complex (real_part, -imag_part); + if (j+1 >= n) + { + (*current_liboctave_error_handler) ("EIG: internal error"); + return -1; + } + + for (int i = 0; i < n; i++) + { + lambda.elem (j) = Complex (wr.elem (j), wi.elem (j)); + lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1)); + double real_part = vr.elem (i, j); + double imag_part = vr.elem (i, j+1); + v.elem (i, j) = Complex (real_part, imag_part); + v.elem (i, j+1) = Complex (real_part, -imag_part); + } + j++; } - j++; } } - delete [] tmp_data; - delete [] wr; - delete [] wi; - delete [] work; - return info; } int EIG::init (const ComplexMatrix& a) { - int a_nr = a.rows (); - if (a_nr != a.cols ()) + int n = a.rows (); + + if (n != a.cols ()) { (*current_liboctave_error_handler) ("EIG requires square matrix"); return -1; } - int n = a_nr; - int info; lambda.resize (n); v.resize (n, n); - + Complex *pw = lambda.fortran_vec (); Complex *pvr = v.fortran_vec (); - Complex *tmp_data = dup (a.data (), a.length ()); + ComplexMatrix atmp = a; + Complex *tmp_data = atmp.fortran_vec (); int lwork = 8*n; - Complex *work = new Complex[lwork]; - double *rwork = new double[4*n]; + + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); + + Array<double> rwork (4*n); + double *prwork = rwork.fortran_vec (); Complex *dummy = 0; int idummy = 1; - F77_FCN (zgeev, ZGEEV) ("N", "V", n, tmp_data, n, pw, dummy, idummy, - pvr, n, work, lwork, rwork, info, 1L, 1L); + F77_XFCN (zgeev, ZGEEV, ("N", "V", n, tmp_data, n, pw, dummy, idummy, + pvr, n, pwork, lwork, prwork, info, 1L, 1L)); - delete [] tmp_data; - delete [] work; - delete [] rwork; + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); return info; }