Mercurial > hg > octave-lojdl
changeset 1932:682f31b20894
[project @ 1996-02-12 02:26:07 by jwe]
author | jwe |
---|---|
date | Mon, 12 Feb 1996 02:26:07 +0000 |
parents | c91f81d5f72c |
children | 333923894848 |
files | liboctave/CmplxHESS.cc liboctave/CmplxHESS.h liboctave/dbleHESS.cc liboctave/dbleHESS.h |
diffstat | 4 files changed, 120 insertions(+), 73 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxHESS.cc +++ b/liboctave/CmplxHESS.cc @@ -60,59 +60,82 @@ { int a_nr = a.rows (); int a_nc = a.cols (); - if (a_nr != a_nc) - { - (*current_liboctave_error_handler) - ("ComplexHESS requires square matrix"); - return -1; - } + + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) + ("ComplexHESS requires square matrix"); + return -1; + } - char *job = "N"; - char *side = "R"; + char job = 'N'; + char side = 'R'; + + int n = a_nc; + int lwork = 32 * n; + int info; + int ilo; + int ihi; - int n = a_nc; - int lwork = 32 * n; - int info; - int ilo; - int ihi; + hess_mat = a; + Complex *h = hess_mat.fortran_vec (); - Complex *h = dup (a.data (), a.length ()); + Array<double> scale (n); + double *pscale = scale.fortran_vec (); + + F77_XFCN (zgebal, ZGEBAL, (&job, n, h, n, ilo, ihi, pscale, info, + 1L, 1L)); - double *scale = new double [n]; - Complex *tau = new Complex [n-1]; - Complex *work = new Complex [lwork]; - Complex *z = new Complex [n*n]; + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); + else + { + Array<Complex> tau (n-1); + Complex *ptau = tau.fortran_vec (); - F77_FCN (zgebal, ZGEBAL) (job, n, h, n, ilo, ihi, scale, info, 1L, - 1L); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); - F77_FCN (zgehrd, ZGEHRD) (n, ilo, ihi, h, n, tau, work, lwork, - info, 1L, 1L); - - copy (z, h, n*n); + F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, + info, 1L, 1L)); - F77_FCN (zunghr, ZUNGHR) (n, ilo, ihi, z, n, tau, work, lwork, - info, 1L, 1L); + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgehrd"); + else + { + unitary_hess_mat = hess_mat; + Complex *z = unitary_hess_mat.fortran_vec (); - F77_FCN (zgebak, ZGEBAK) (job, side, n, ilo, ihi, scale, n, z, n, - info, 1L, 1L); - - hess_mat = ComplexMatrix (h, n, n); - unitary_hess_mat = ComplexMatrix (z, n, n); + F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info, 1L, 1L)); - // If someone thinks of a more graceful way of doing this (or faster - // for that matter :-)), please let me know! + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zunghr"); + else + { + F77_XFCN (zgebak, ZGEBAK, (&job, &side, n, ilo, ihi, + pscale, n, z, n, info, 1L, 1L)); - if (n > 2) - for (int j = 0; j < a_nc; j++) - for (int i = j+2; i < a_nr; i++) - hess_mat.elem (i, j) = 0; + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgebak"); + else + { + // If someone thinks of a more graceful way of + // doing this (or faster for that matter :-)), + // please let me know! - delete [] work; - delete [] tau; - delete [] scale; + if (n > 2) + for (int j = 0; j < a_nc; j++) + for (int i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; + } + } + } + } - return info; + return info; } /*
--- a/liboctave/CmplxHESS.h +++ b/liboctave/CmplxHESS.h @@ -56,6 +56,8 @@ return *this; } + ~ComplexHESS (void) { } + ComplexMatrix hess_matrix (void) const { return hess_mat; } ComplexMatrix unitary_hess_matrix (void) const
--- a/liboctave/dbleHESS.cc +++ b/liboctave/dbleHESS.cc @@ -59,14 +59,15 @@ { int a_nr = a.rows (); int a_nc = a.cols (); + if (a_nr != a_nc) { (*current_liboctave_error_handler) ("HESS requires square matrix"); return -1; } - char *jobbal = "N"; - char *side = "R"; + char job = 'N'; + char side = 'R'; int n = a_nc; int lwork = 32 * n; @@ -74,44 +75,63 @@ int ilo; int ihi; - double *h = dup (a.data (), a.length ()); + hess_mat = a; + double *h = hess_mat.fortran_vec (); - double *tau = new double [n+1]; - double *scale = new double [n]; - double *z = new double [n*n]; - double *work = new double [lwork]; + Array<double> scale (n); + double *pscale = scale.fortran_vec (); + + F77_XFCN (dgebal, DGEBAL, (&job, n, h, n, ilo, ihi, pscale, info, + 1L, 1L)); - F77_FCN (dgebal, DGEBAL) (jobbal, n, h, n, ilo, ihi, scale, info, - 1L, 1L); + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); + else + { + Array<double> tau (n-1); + double *ptau = tau.fortran_vec (); + + Array<double> work (lwork); + double *pwork = work.fortran_vec (); - F77_FCN (dgehrd, DGEHRD) (n, ilo, ihi, h, n, tau, work, lwork, info, - 1L, 1L); + F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork, + lwork, info, 1L, 1L)); - copy (z, h, n*n); - - F77_FCN (dorghr, DORGHR) (n, ilo, ihi, z, n, tau, work, lwork, info, - 1L, 1L); + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgehrd"); + else + { + unitary_hess_mat = hess_mat; + double *z = unitary_hess_mat.fortran_vec (); - F77_FCN (dgebak, DGEBAK) (jobbal, side, n, ilo, ihi, scale, n, z, n, - info, 1L, 1L); + F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info, 1L, 1L)); - // We need to clear out all of the area below the sub-diagonal which - // was used to store the unitary matrix. - - hess_mat = Matrix (h, n, n); - unitary_hess_mat = Matrix (z, n, n); + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dorghr"); + else + { + F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, + pscale, n, z, n, info, 1L, 1L)); - // If someone thinks of a more graceful way of doing this (or faster - // for that matter :-)), please let me know! + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgebak"); + else + { + // If someone thinks of a more graceful way of doing + // this (or faster for that matter :-)), please let + // me know! - if (n > 2) - for (int j = 0; j < a_nc; j++) - for (int i = j+2; i < a_nr; i++) - hess_mat.elem (i, j) = 0; - - delete [] tau; - delete [] work; - delete [] scale; + if (n > 2) + for (int j = 0; j < a_nc; j++) + for (int i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; + } + } + } + } return info; }