Mercurial > hg > octave-lojdl
changeset 4800:c322edde72ac
[project @ 2004-02-25 03:54:33 by jwe]
author | jwe |
---|---|
date | Wed, 25 Feb 2004 03:54:33 +0000 |
parents | e2d7d1ef5e55 |
children | b022780ac0b4 |
files | liboctave/ChangeLog liboctave/EIG.cc |
diffstat | 2 files changed, 142 insertions(+), 74 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2004-02-24 John W. Eaton <jwe@bevo.che.wisc.edu> + + * EIG.cc (EIG::init, EIG::symmetric_init): + Query Lapack for workspace size. + 2004-02-23 John W. Eaton <jwe@bevo.che.wisc.edu> * Array.cc (Array<T>::resize_and_fill (const dim_vector&, const T&)):
--- a/liboctave/EIG.cc +++ b/liboctave/EIG.cc @@ -99,12 +99,8 @@ Matrix vr (nvr, nvr); double *pvr = vr.fortran_vec (); - // XXX FIXME XXX -- it might be possible to choose a better value of - // lwork that would result in more efficient computations. - - int lwork = 8*n; - Array<double> work (lwork); - double *pwork = work.fortran_vec (); + int lwork = -1; + double dummy_work; double *dummy = 0; int idummy = 1; @@ -112,53 +108,70 @@ F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), n, tmp_data, n, pwr, pwi, dummy, - idummy, pvr, n, pwork, lwork, info + idummy, pvr, n, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); - else + if (! f77_exception_encountered && info == 0) { + lwork = static_cast<int> (dummy_work); + Array<double> work (lwork); + double *pwork = work.fortran_vec (); + + F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, tmp_data, n, pwr, pwi, dummy, + idummy, pvr, n, pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered || info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); + return info; + } + if (info > 0) - (*current_liboctave_error_handler) ("dgeev failed to converge"); - else { - lambda.resize (n); - v.resize (nvr, nvr); + (*current_liboctave_error_handler) ("dgeev failed to converge"); + return info; + } - for (int j = 0; j < n; j++) + lambda.resize (n); + v.resize (nvr, nvr); + + for (int j = 0; j < n; j++) + { + if (wi.elem (j) == 0.0) { - if (wi.elem (j) == 0.0) + lambda.elem (j) = Complex (wr.elem (j)); + for (int i = 0; i < nvr; i++) + v.elem (i, j) = vr.elem (i, j); + } + else + { + if (j+1 >= n) { - lambda.elem (j) = Complex (wr.elem (j)); - for (int i = 0; i < nvr; i++) - v.elem (i, j) = vr.elem (i, j); + (*current_liboctave_error_handler) ("EIG: internal error"); + return -1; } - else - { - if (j+1 >= n) - { - (*current_liboctave_error_handler) - ("EIG: internal error"); - return -1; - } + + lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); + lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); - lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); - lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); - - for (int i = 0; i < nvr; i++) - { - 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++; + for (int i = 0; i < nvr; i++) + { + 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++; } } } + else + (*current_liboctave_error_handler) ("dgeev workspace query failed"); return info; } @@ -182,28 +195,44 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - // XXX FIXME XXX -- it might be possible to choose a better value of - // lwork that would result in more efficient computations. - - int lwork = 8*n; - Array<double> work (lwork); - double *pwork = work.fortran_vec (); + int lwork = -1; + double dummy_work; F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, tmp_data, n, pwr, pwork, lwork, info + n, tmp_data, n, pwr, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); - else if (info > 0) - (*current_liboctave_error_handler) ("dsyev failed to converge"); - else + if (! f77_exception_encountered && info == 0) { + lwork = static_cast<int> (dummy_work); + Array<double> work (lwork); + double *pwork = work.fortran_vec (); + + F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, tmp_data, n, pwr, pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered || info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("dsyev failed to converge"); + return info; + } + lambda = ComplexColumnVector (wr); v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); } + else + (*current_liboctave_error_handler) ("dsyev workspace query failed"); return info; } @@ -234,12 +263,8 @@ ComplexMatrix vtmp (nvr, nvr); Complex *pv = vtmp.fortran_vec (); - // XXX FIXME XXX -- it might be possible to choose a better value of - // lwork that would result in more efficient computations. - - int lwork = 8*n; - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); + int lwork = -1; + Complex dummy_work; int lrwork = 2*n; Array<double> rwork (lrwork); @@ -251,19 +276,40 @@ F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), n, tmp_data, n, pw, dummy, idummy, - pv, n, pwork, lwork, prwork, info + pv, n, &dummy_work, lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); - else if (info > 0) - (*current_liboctave_error_handler) ("zgeev failed to converge"); - else + if (! f77_exception_encountered && info == 0) { + lwork = static_cast<int> (dummy_work.real ()); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); + + F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, tmp_data, n, pw, dummy, idummy, + pv, n, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered || info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("zgeev failed to converge"); + return info; + } + lambda = w; v = vtmp; } + else + (*current_liboctave_error_handler) ("zgeev workspace query failed"); return info; } @@ -287,12 +333,8 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - // XXX FIXME XXX -- it might be possible to choose a better value of - // lwork that would result in more efficient computations. - - int lwork = 8*n; - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); + int lwork = -1; + Complex dummy_work; int lrwork = 3*n; Array<double> rwork (lrwork); @@ -300,19 +342,40 @@ F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, tmp_data, n, pwr, pwork, lwork, prwork, info + n, tmp_data, n, pwr, &dummy_work, lwork, + prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in zheev"); - else if (info > 0) - (*current_liboctave_error_handler) ("zheev failed to converge"); - else + if (! f77_exception_encountered && info == 0) { + lwork = static_cast<int> (dummy_work.real ()); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); + + F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, tmp_data, n, pwr, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered || info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in zheev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("zheev failed to converge"); + return info; + } + lambda = ComplexColumnVector (wr); v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); } + else + (*current_liboctave_error_handler) ("zheev workspace query failed"); return info; }