Mercurial > hg > octave-lojdl
changeset 8339:18c4ded8612a
Add generalized eigenvalue functions
author | Jarkko Kaleva <d3roga@gmail.com> |
---|---|
date | Mon, 24 Nov 2008 10:55:50 +0100 |
parents | a35bf360b919 |
children | fa9e6619fa99 |
files | liboctave/ChangeLog liboctave/EIG.cc liboctave/EIG.h liboctave/fEIG.cc liboctave/fEIG.h src/ChangeLog src/DLD-FUNCTIONS/eig.cc |
diffstat | 7 files changed, 1196 insertions(+), 39 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,41 @@ +2008-11-21 Jarkko Kaleva <d3roga@gmail.com> + + * EIG.h (EIG::EIG (const Matrix& a, const Matrix& b, + bool calc_eigenvectors = true)): New constructor. + (EIG::EIG (const Matrix& a, const Matrix& b, octave_idx_type& info, + bool calc_eigenvectors = true)): New constructor. + (EIG::EIG (const ComplexMatrix& a, const ComplexMatrix& b, + bool calc_eigenvectors = true)): New constructor. + (EIG::EIG (const ComplexMatrix& a, const ComplexMatrix& b, + octave_idx_type& info, bool calc_eigenvectors = true)): New + constructor. + * EIG.cc (EIG::init (const Matrix& a, const Matrix& b, + bool calc_eigenvectors)): New function. + (EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, + bool calc_eigenvectors)): New function. + (EIG::symmetric_init (const Matrix& a, const Matrix& b, + bool calc_eigenvectors)): New function. + (EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, + bool calc_eigenvectors)): New function. + * fEIG.h (fEIG::fEIG (const FloatMatrix& a, const FloatMatrix& b, + bool calc_eigenvectors = true)): New constructor. + (fEIG::fEIG (const FloatMatrix& a, const FloatMatrix& b, + octave_idx_type& info, bool calc_eigenvectors = true)): New + constructor. + (fEIG::fEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, + bool calc_eigenvectors = true)): New constructor. + (fEIG::fEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, + octave_idx_type& info, bool calc_eigenvectors = true)): New + constructor. + (fEIG::init (const FloatMatrix& a, const FloatMatrix& b, + bool calc_eigenvectors)): New function. + (fEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, + bool calc_eigenvectors)): New function. + (fEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, + bool calc_eigenvectors)): New function. + (fEIG::hermitian_init (const FloatComplexMatrix& a, + const FloatComplexMatrix& b, bool calc_eigenvectors)): New function. + 2008-11-19 Jaroslav Hajek <highegg@gmail.com> * dMatrix.cc (Matrix::determinant),
--- a/liboctave/EIG.cc +++ b/liboctave/EIG.cc @@ -65,6 +65,68 @@ Complex*, const octave_idx_type&, double*, octave_idx_type& F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + double*, const octave_idx_type&, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + Complex*, const octave_idx_type&, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + double*, const octave_idx_type&, + double*, const octave_idx_type&, + double*, double*, double *, + double*, const octave_idx_type&, double*, const octave_idx_type&, + double*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dsygv, DSYGV) (const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + double*, const octave_idx_type&, + double*, const octave_idx_type&, + double*, double*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + Complex*, const octave_idx_type&, + Complex*, const octave_idx_type&, + Complex*, Complex*, + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, const octave_idx_type&, double*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + Complex*, const octave_idx_type&, + Complex*, const octave_idx_type&, + double*, Complex*, const octave_idx_type&, double*, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); } octave_idx_type @@ -391,6 +453,414 @@ return info; } +octave_idx_type +EIG::init (const Matrix& a, const Matrix& b, bool calc_ev) +{ + if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) + { + (*current_liboctave_error_handler) + ("EIG: matrix contains Inf or NaN values"); + return -1; + } + + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + Matrix tmp = b; + double *tmp_data = tmp.fortran_vec (); + + F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (a.is_symmetric () && b.is_symmetric () && info == 0) + return symmetric_init (a, b, calc_ev); + + Matrix atmp = a; + double *atmp_data = atmp.fortran_vec (); + + Matrix btmp = b; + double *btmp_data = btmp.fortran_vec (); + + Array<double> ar (n); + double *par = ar.fortran_vec (); + + Array<double> ai (n); + double *pai = ai.fortran_vec (); + + Array<double> beta (n); + double *pbeta = beta.fortran_vec (); + + volatile octave_idx_type nvr = calc_ev ? n : 0; + Matrix vr (nvr, nvr); + double *pvr = vr.fortran_vec (); + + octave_idx_type lwork = -1; + double dummy_work; + + double *dummy = 0; + octave_idx_type idummy = 1; + + F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + par, pai, pbeta, + dummy, idummy, pvr, n, + &dummy_work, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work); + Array<double> work (lwork); + double *pwork = work.fortran_vec (); + + F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + par, pai, pbeta, + dummy, idummy, pvr, n, + pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in dggev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("dggev failed to converge"); + return info; + } + + lambda.resize (n); + v.resize (nvr, nvr); + + for (octave_idx_type j = 0; j < n; j++) + { + if (ai.elem (j) == 0.0) + { + lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); + for (octave_idx_type i = 0; i < nvr; i++) + v.elem (i, j) = vr.elem (i, j); + } + else + { + if (j+1 >= n) + { + (*current_liboctave_error_handler) ("EIG: internal error"); + return -1; + } + + lambda.elem(j) = Complex (ar.elem(j) / beta.elem (j), + ai.elem(j) / beta.elem (j)); + lambda.elem(j+1) = Complex (ar.elem(j+1) / beta.elem (j+1), + ai.elem(j+1) / beta.elem (j+1)); + + for (octave_idx_type 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) ("dggev workspace query failed"); + + return info; +} + +octave_idx_type +EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev) +{ + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + Matrix atmp = a; + double *atmp_data = atmp.fortran_vec (); + + Matrix btmp = b; + double *btmp_data = btmp.fortran_vec (); + + ColumnVector wr (n); + double *pwr = wr.fortran_vec (); + + octave_idx_type lwork = -1; + double dummy_work; + + F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, &dummy_work, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work); + Array<double> work (lwork); + double *pwork = work.fortran_vec (); + + F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in dsygv"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("dsygv failed to converge"); + return info; + } + + lambda = ComplexColumnVector (wr); + v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + } + else + (*current_liboctave_error_handler) ("dsygv workspace query failed"); + + return info; +} + +octave_idx_type +EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) +{ + if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) + { + (*current_liboctave_error_handler) + ("EIG: matrix contains Inf or NaN values"); + return -1; + } + + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + ComplexMatrix tmp = b; + Complex*tmp_data = tmp.fortran_vec (); + + F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (a.is_hermitian () && b.is_hermitian () && info == 0) + return hermitian_init (a, calc_ev); + + ComplexMatrix atmp = a; + Complex *atmp_data = atmp.fortran_vec (); + + ComplexMatrix btmp = b; + Complex *btmp_data = btmp.fortran_vec (); + + ComplexColumnVector alpha (n); + Complex *palpha = alpha.fortran_vec (); + + ComplexColumnVector beta (n); + Complex *pbeta = beta.fortran_vec (); + + octave_idx_type nvr = calc_ev ? n : 0; + ComplexMatrix vtmp (nvr, nvr); + Complex *pv = vtmp.fortran_vec (); + + octave_idx_type lwork = -1; + Complex dummy_work; + + octave_idx_type lrwork = 8*n; + Array<double> rwork (lrwork); + double *prwork = rwork.fortran_vec (); + + Complex *dummy = 0; + octave_idx_type idummy = 1; + + F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + palpha, pbeta, dummy, idummy, + pv, n, &dummy_work, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work.real ()); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); + + F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + palpha, pbeta, dummy, idummy, + pv, n, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in zggev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("zggev failed to converge"); + return info; + } + + lambda.resize (n); + + for (octave_idx_type j = 0; j < n; j++) + lambda.elem (j) = alpha.elem (j) / beta.elem(j); + + v = vtmp; + } + else + (*current_liboctave_error_handler) ("zggev workspace query failed"); + + return info; +} + +octave_idx_type +EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) +{ + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + ComplexMatrix atmp = a; + Complex *atmp_data = atmp.fortran_vec (); + + ComplexMatrix btmp = b; + Complex *btmp_data = btmp.fortran_vec (); + + ColumnVector wr (n); + double *pwr = wr.fortran_vec (); + + octave_idx_type lwork = -1; + Complex dummy_work; + + octave_idx_type lrwork = 3*n; + Array<double> rwork (lrwork); + double *prwork = rwork.fortran_vec (); + + F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, &dummy_work, lwork, + prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work.real ()); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); + + F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("zhegv failed to converge"); + return info; + } + + lambda = ComplexColumnVector (wr); + v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + } + else + (*current_liboctave_error_handler) ("zhegv workspace query failed"); + + return info; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/EIG.h +++ b/liboctave/EIG.h @@ -48,12 +48,24 @@ EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true) { info = init (a, calc_eigenvectors); } + EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true) + { init (a, b, calc_eigenvectors); } + + EIG (const Matrix& a, const Matrix& b, octave_idx_type& info, bool calc_eigenvectors = true) + { info = init (a, b, calc_eigenvectors); } + EIG (const ComplexMatrix& a, bool calc_eigenvectors = true) { init (a, calc_eigenvectors); } EIG (const ComplexMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true) { info = init (a, calc_eigenvectors); } + EIG (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors = true) + { init (a, b, calc_eigenvectors); } + + EIG (const ComplexMatrix& a, const ComplexMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true) + { info = init (a, b, calc_eigenvectors); } + EIG (const EIG& a) : lambda (a.lambda), v (a.v) { } @@ -81,10 +93,14 @@ ComplexMatrix v; octave_idx_type init (const Matrix& a, bool calc_eigenvectors); + octave_idx_type init (const Matrix& a, const Matrix& b, bool calc_eigenvectors); octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors); + octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors); octave_idx_type symmetric_init (const Matrix& a, bool calc_eigenvectors); + octave_idx_type symmetric_init (const Matrix& a, const Matrix& b, bool calc_eigenvectors); octave_idx_type hermitian_init (const ComplexMatrix& a, bool calc_eigenvectors); + octave_idx_type hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors); }; #endif
--- a/liboctave/fEIG.cc +++ b/liboctave/fEIG.cc @@ -65,6 +65,68 @@ FloatComplex*, const octave_idx_type&, float*, octave_idx_type& F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + float*, const octave_idx_type&, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sggev, SGGEV) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + float*, const octave_idx_type&, + float*, const octave_idx_type&, + float*, float*, float*, + float*, const octave_idx_type&, float*, const octave_idx_type&, + float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (ssygv, SSYGV) (const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + float*, const octave_idx_type&, + float*, const octave_idx_type&, + float*, float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cggev, CGGEV) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + FloatComplex*, FloatComplex*, + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (chegv, CHEGV) (const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + float*, FloatComplex*, const octave_idx_type&, float*, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); } octave_idx_type @@ -391,6 +453,414 @@ return info; } +octave_idx_type +FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev) +{ + if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) + { + (*current_liboctave_error_handler) + ("EIG: matrix contains Inf or NaN values"); + return -1; + } + + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + FloatMatrix tmp = b; + float *tmp_data = tmp.fortran_vec (); + + F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (a.is_symmetric () && b.is_symmetric () && info == 0) + return symmetric_init (a, b, calc_ev); + + FloatMatrix atmp = a; + float *atmp_data = atmp.fortran_vec (); + + FloatMatrix btmp = b; + float *btmp_data = btmp.fortran_vec (); + + Array<float> ar (n); + float *par = ar.fortran_vec (); + + Array<float> ai (n); + float *pai = ai.fortran_vec (); + + Array<float> beta (n); + float *pbeta = beta.fortran_vec (); + + volatile octave_idx_type nvr = calc_ev ? n : 0; + FloatMatrix vr (nvr, nvr); + float *pvr = vr.fortran_vec (); + + octave_idx_type lwork = -1; + float dummy_work; + + float *dummy = 0; + octave_idx_type idummy = 1; + + F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + par, pai, pbeta, + dummy, idummy, pvr, n, + &dummy_work, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work); + Array<float> work (lwork); + float *pwork = work.fortran_vec (); + + F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + par, pai, pbeta, + dummy, idummy, pvr, n, + pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in sggev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("sggev failed to converge"); + return info; + } + + lambda.resize (n); + v.resize (nvr, nvr); + + for (octave_idx_type j = 0; j < n; j++) + { + if (ai.elem (j) == 0.0) + { + lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j)); + for (octave_idx_type i = 0; i < nvr; i++) + v.elem (i, j) = vr.elem (i, j); + } + else + { + if (j+1 >= n) + { + (*current_liboctave_error_handler) ("EIG: internal error"); + return -1; + } + + lambda.elem(j) = FloatComplex (ar.elem(j) / beta.elem (j), + ai.elem(j) / beta.elem (j)); + lambda.elem(j+1) = FloatComplex (ar.elem(j+1) / beta.elem (j+1), + ai.elem(j+1) / beta.elem (j+1)); + + for (octave_idx_type i = 0; i < nvr; i++) + { + float real_part = vr.elem (i, j); + float imag_part = vr.elem (i, j+1); + v.elem (i, j) = FloatComplex (real_part, imag_part); + v.elem (i, j+1) = FloatComplex (real_part, -imag_part); + } + j++; + } + } + } + else + (*current_liboctave_error_handler) ("sggev workspace query failed"); + + return info; +} + +octave_idx_type +FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev) +{ + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + FloatMatrix atmp = a; + float *atmp_data = atmp.fortran_vec (); + + FloatMatrix btmp = b; + float *btmp_data = btmp.fortran_vec (); + + FloatColumnVector wr (n); + float *pwr = wr.fortran_vec (); + + octave_idx_type lwork = -1; + float dummy_work; + + F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, &dummy_work, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work); + Array<float> work (lwork); + float *pwork = work.fortran_vec (); + + F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, pwork, lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in ssygv"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("ssygv failed to converge"); + return info; + } + + lambda = FloatComplexColumnVector (wr); + v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + } + else + (*current_liboctave_error_handler) ("ssygv workspace query failed"); + + return info; +} + +octave_idx_type +FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_ev) +{ + if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) + { + (*current_liboctave_error_handler) + ("EIG: matrix contains Inf or NaN values"); + return -1; + } + + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + FloatComplexMatrix tmp = b; + FloatComplex *tmp_data = tmp.fortran_vec (); + + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (a.is_hermitian () && b.is_hermitian () && info == 0) + return hermitian_init (a, calc_ev); + + FloatComplexMatrix atmp = a; + FloatComplex *atmp_data = atmp.fortran_vec (); + + FloatComplexMatrix btmp = b; + FloatComplex *btmp_data = btmp.fortran_vec (); + + FloatComplexColumnVector alpha (n); + FloatComplex *palpha = alpha.fortran_vec (); + + FloatComplexColumnVector beta (n); + FloatComplex *pbeta = beta.fortran_vec (); + + octave_idx_type nvr = calc_ev ? n : 0; + FloatComplexMatrix vtmp (nvr, nvr); + FloatComplex *pv = vtmp.fortran_vec (); + + octave_idx_type lwork = -1; + FloatComplex dummy_work; + + octave_idx_type lrwork = 8*n; + Array<float> rwork (lrwork); + float *prwork = rwork.fortran_vec (); + + FloatComplex *dummy = 0; + octave_idx_type idummy = 1; + + F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + palpha, pbeta, dummy, idummy, + pv, n, &dummy_work, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work.real ()); + Array<FloatComplex> work (lwork); + FloatComplex *pwork = work.fortran_vec (); + + F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + n, atmp_data, n, btmp_data, n, + palpha, pbeta, dummy, idummy, + pv, n, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in cggev"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("cggev failed to converge"); + return info; + } + + lambda.resize (n); + + for (octave_idx_type j = 0; j < n; j++) + lambda.elem (j) = alpha.elem (j) / beta.elem(j); + + v = vtmp; + } + else + (*current_liboctave_error_handler) ("cggev workspace query failed"); + + return info; +} + +octave_idx_type +FloatEIG::hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_ev) +{ + octave_idx_type n = a.rows (); + octave_idx_type nb = b.rows (); + + if (n != a.cols () || nb != b.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + if (n != nb) + { + (*current_liboctave_error_handler) ("EIG requires same size matrices"); + return -1; + } + + octave_idx_type info = 0; + + FloatComplexMatrix atmp = a; + FloatComplex *atmp_data = atmp.fortran_vec (); + + FloatComplexMatrix btmp = b; + FloatComplex *btmp_data = btmp.fortran_vec (); + + FloatColumnVector wr (n); + float *pwr = wr.fortran_vec (); + + octave_idx_type lwork = -1; + FloatComplex dummy_work; + + octave_idx_type lrwork = 3*n; + Array<float> rwork (lrwork); + float *prwork = rwork.fortran_vec (); + + F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, &dummy_work, lwork, + prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info == 0) + { + lwork = static_cast<octave_idx_type> (dummy_work.real ()); + Array<FloatComplex> work (lwork); + FloatComplex *pwork = work.fortran_vec (); + + F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("U", 1), + n, atmp_data, n, + btmp_data, n, + pwr, pwork, lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info < 0) + { + (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); + return info; + } + + if (info > 0) + { + (*current_liboctave_error_handler) ("zhegv failed to converge"); + return info; + } + + lambda = FloatComplexColumnVector (wr); + v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + } + else + (*current_liboctave_error_handler) ("zhegv workspace query failed"); + + return info; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/fEIG.h +++ b/liboctave/fEIG.h @@ -48,12 +48,24 @@ FloatEIG (const FloatMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true) { info = init (a, calc_eigenvectors); } + FloatEIG (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors = true) + { init (a, b, calc_eigenvectors); } + + FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true) + { info = init (a, b, calc_eigenvectors); } + FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true) { init (a, calc_eigenvectors); } FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true) { info = init (a, calc_eigenvectors); } + FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors = true) + { init (a, b, calc_eigenvectors); } + + FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true) + { info = init (a, b, calc_eigenvectors); } + FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v) { } @@ -81,10 +93,14 @@ FloatComplexMatrix v; octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors); + octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors); octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors); + octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors); octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_eigenvectors); + octave_idx_type symmetric_init (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors); octave_idx_type hermitian_init (const FloatComplexMatrix& a, bool calc_eigenvectors); + octave_idx_type hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors); }; #endif
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2008-11-21 Jarkko Kaleva <d3roga@gmail.com> + + * DLD-FUNCTIONS/eig.cc (Feig): Handle generalized eigenvalues and + eigenvectors. + 2008-11-19 Jaroslav Hajek <highegg@gmail.com> * DLD_FUNCTIONS/det.cc: Include only DET.h. Retrieve & matrix type &
--- a/src/DLD-FUNCTIONS/eig.cc +++ b/src/DLD-FUNCTIONS/eig.cc @@ -37,7 +37,9 @@ DEFUN_DLD (eig, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{lambda} =} eig (@var{a})\n\ +@deftypefnx {Loadable Function} {@var{lambda} =} eig (@var{a}, @var{b})\n\ @deftypefnx {Loadable Function} {[@var{v}, @var{lambda}] =} eig (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{v}, @var{lambda}] =} eig (@var{a}, @var{b})\n\ The eigenvalues (and eigenvectors) of a matrix are computed in a several\n\ step process which begins with a Hessenberg decomposition, followed by a\n\ Schur decomposition, from which the eigenvalues are apparent. The\n\ @@ -51,55 +53,116 @@ int nargin = args.length (); - if (nargin != 1 || nargout > 2) + if (nargin > 2 || nargin == 0 || nargout > 2) { print_usage (); return retval; } - octave_value arg = args(0); + octave_value arg_a, arg_b; + + octave_idx_type nr_a = 0, nr_b = 0; + octave_idx_type nc_a = 0, nc_b = 0; - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); + arg_a = args(0); + nr_a = arg_a.rows (); + nc_a = arg_a.columns (); - int arg_is_empty = empty_arg ("eig", nr, nc); + int arg_is_empty = empty_arg ("eig", nr_a, nc_a); if (arg_is_empty < 0) return retval; else if (arg_is_empty > 0) return octave_value_list (2, Matrix ()); - if (nr != nc) + if (!(arg_a.is_single_type () || arg_a.is_double_type ())) + { + gripe_wrong_type_arg ("eig", arg_a); + return retval; + } + + if (nargin == 2) + { + arg_b = args(1); + nr_b = arg_b.rows (); + nc_b = arg_b.columns (); + + arg_is_empty = empty_arg ("eig", nr_b, nc_b); + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (2, Matrix ()); + + if (!(arg_b.is_single_type() || arg_b.is_double_type ())) + { + gripe_wrong_type_arg ("eig", arg_b); + return retval; + } + } + + if (nr_a != nc_a) { gripe_square_matrix_required ("eig"); return retval; } - Matrix tmp; - ComplexMatrix ctmp; - FloatMatrix ftmp; - FloatComplexMatrix fctmp; + if (nargin == 2 && nr_b != nc_b) + { + gripe_square_matrix_required ("eig"); + return retval; + } - if (arg.is_single_type ()) + Matrix tmp_a, tmp_b; + ComplexMatrix ctmp_a, ctmp_b; + FloatMatrix ftmp_a, ftmp_b; + FloatComplexMatrix fctmp_a, fctmp_b; + + if (arg_a.is_single_type ()) { FloatEIG result; - if (arg.is_real_type ()) + if (nargin == 1) { - ftmp = arg.float_matrix_value (); + if (arg_a.is_real_type ()) + { + ftmp_a = arg_a.float_matrix_value (); - if (error_state) - return retval; + if (error_state) + return retval; + else + result = FloatEIG (ftmp_a, nargout > 1); + } else - result = FloatEIG (ftmp, nargout > 1); + { + fctmp_a = arg_a.float_complex_matrix_value (); + + if (error_state) + return retval; + else + result = FloatEIG (fctmp_a, nargout > 1); + } } - else if (arg.is_complex_type ()) + else if (nargin == 2) { - fctmp = arg.float_complex_matrix_value (); + if (arg_a.is_real_type () && arg_b.is_real_type ()) + { + ftmp_a = arg_a.float_matrix_value (); + ftmp_b = arg_b.float_matrix_value (); - if (error_state) - return retval; + if (error_state) + return retval; + else + result = FloatEIG (ftmp_a, ftmp_b, nargout > 1); + } else - result = FloatEIG (fctmp, nargout > 1); + { + fctmp_a = arg_a.float_complex_matrix_value (); + fctmp_b = arg_b.float_complex_matrix_value (); + + if (error_state) + return retval; + else + result = FloatEIG (fctmp_a, fctmp_b, nargout > 1); + } } if (! error_state) @@ -123,28 +186,49 @@ { EIG result; - if (arg.is_real_type ()) - { - tmp = arg.matrix_value (); - - if (error_state) - return retval; - else - result = EIG (tmp, nargout > 1); - } - else if (arg.is_complex_type ()) + if (nargin == 1) { - ctmp = arg.complex_matrix_value (); + if (arg_a.is_real_type ()) + { + tmp_a = arg_a.matrix_value (); - if (error_state) - return retval; + if (error_state) + return retval; + else + result = EIG (tmp_a, nargout > 1); + } else - result = EIG (ctmp, nargout > 1); + { + ctmp_a = arg_a.complex_matrix_value (); + + if (error_state) + return retval; + else + result = EIG (ctmp_a, nargout > 1); + } } - else + else if (nargin == 2) { - gripe_wrong_type_arg ("eig", tmp); - return retval; + if (arg_a.is_real_type () && arg_b.is_real_type ()) + { + tmp_a = arg_a.matrix_value (); + tmp_b = arg_b.matrix_value (); + + if (error_state) + return retval; + else + result = EIG (tmp_a, tmp_b, nargout > 1); + } + else + { + ctmp_a = arg_a.complex_matrix_value (); + ctmp_b = arg_b.complex_matrix_value (); + + if (error_state) + return retval; + else + result = EIG (ctmp_a, ctmp_b, nargout > 1); + } } if (! error_state) @@ -186,9 +270,67 @@ %! assert(d, single([-1, 0; 0, 3]), sqrt (eps('single'))); %! assert(v, [-x, x; x, x], sqrt (eps('single'))); +%!test +%! A = [1, 2; -1, 1]; B = [3, 3; 1, 2]; +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); + +%!test +%! A = single([1, 2; -1, 1]); B = single([3, 3; 1, 2]); +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single'))); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single'))); + +%!test +%! A = [1, 2; 2, 1]; B = [3, -2; -2, 3]; +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); + +%!test +%! A = single([1, 2; 2, 1]); B = single([3, -2; -2, 3]); +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single'))); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single'))); + +%!test +%! A = [1+3i, 2+i; 2-i, 1+3i]; B = [5+9i, 2+i; 2-i, 5+9i]; +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); + +%!test +%! A = single([1+3i, 2+i; 2-i, 1+3i]); B = single([5+9i, 2+i; 2-i, 5+9i]); +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single'))); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single'))); + +%!test +%! A = [1+3i, 2+3i; 3-8i, 8+3i]; B = [8+i, 3+i; 4-9i, 3+i]; +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); + +%!test +%! A = single([1+3i, 2+3i; 3-8i, 8+3i]); B = single([8+i, 3+i; 4-9i, 3+i]); +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single'))); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single'))); + +%!test +%! A = [1, 2; 3, 8]; B = [8, 3; 4, 3]; +%! [v, d] = eig (A, B); +%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); +%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); + %!error <Invalid call to eig.*> eig (); -%!error <Invalid call to eig.*> eig ([1, 2; 3, 4], 2); +%!error <Invalid call to eig.*> eig ([1, 2; 3, 4], [4, 3; 2, 1], 1); +%!error eig ([1, 2; 3, 4], 2); %!error eig ([1, 2; 3, 4; 5, 6]); +%!error eig ("abcd"); +%!error eig ([1 2 ; 2 3], "abcd"); +%!error eig (false, [1 2 ; 2 3]); */