Mercurial > hg > octave-nkf
changeset 462:07fabd96ac6a
[project @ 1994-06-06 00:58:18 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 06 Jun 1994 00:58:18 +0000 |
parents | 00f8b2242a18 |
children | 0bfa99f634be |
files | liboctave/EIG.cc liboctave/EIG.h |
diffstat | 2 files changed, 293 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/liboctave/EIG.cc @@ -0,0 +1,169 @@ +// -*- C++ -*- +/* + +Copyright (C) 1992, 1993, 1994 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#if defined (__GNUG__) +#pragma implementation +#endif + +#include "EIG.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgeev) (const char*, const char*, const int*, double*, + const int*, double*, double*, double*, + const int*, double*, const int*, double*, + const int*, int*, long, long); + + int F77_FCN (zgeev) (const char*, const char*, const int*, Complex*, + const int*, Complex*, Complex*, const int*, + Complex*, const int*, Complex*, const int*, + double*, int*, long, long); +} + +int +EIG::init (const Matrix& a) +{ + int a_nr = a.rows (); + if (a_nr != a.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + int n = a_nr; + + int info; + + char jobvl = 'N'; + char jobvr = 'V'; + + double *tmp_data = dup (a.data (), a.length ()); + double *wr = new double[n]; + double *wi = new double[n]; + Matrix vr (n, n); + double *pvr = vr.fortran_vec (); + int lwork = 8*n; + double *work = new double[lwork]; + + double dummy; + int idummy = 1; + + F77_FCN (dgeev) (&jobvl, &jobvr, &n, tmp_data, &n, wr, wi, &dummy, + &idummy, pvr, &n, work, &lwork, &info, 1L, 1L); + + lambda.resize (n); + v.resize (n, n); + + for (int j = 0; j < n; j++) + { + 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 + { + 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[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); + } + 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 ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + int n = a_nr; + + int info; + + char jobvl = 'N'; + char jobvr = 'V'; + + 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 ()); + + int lwork = 8*n; + Complex *work = new Complex[lwork]; + double *rwork = new double[4*n]; + + Complex dummy; + int idummy = 1; + + F77_FCN (zgeev) (&jobvl, &jobvr, &n, tmp_data, &n, pw, &dummy, + &idummy, pvr, &n, work, &lwork, rwork, &info, 1L, + 1L); + + delete [] tmp_data; + delete [] work; + delete [] rwork; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/EIG.h @@ -0,0 +1,124 @@ +// -*- C++ -*- +/* + +Copyright (C) 1992, 1993, 1994 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#if !defined (octave_EIG_h) +#define octave_EIG_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" +#include "CMatrix.h" +#include "CColVector.h" + +extern "C++" { + +class EIG +{ +friend class Matrix; +friend class ComplexMatrix; + +public: + + EIG (void) {} + + EIG (const Matrix& a); + EIG (const Matrix& a, int& info); + + EIG (const ComplexMatrix& a); + EIG (const ComplexMatrix& a, int& info); + + EIG (const EIG& a); + + EIG& operator = (const EIG& a); + + ComplexColumnVector eigenvalues (void) const; + ComplexMatrix eigenvectors (void) const; + + friend ostream& operator << (ostream& os, const EIG& a); + +private: + + int init (const Matrix& a); + int init (const ComplexMatrix& a); + + ComplexColumnVector lambda; + ComplexMatrix v; +}; + +inline EIG::EIG (const Matrix& a) +{ + init (a); +} + +inline EIG::EIG (const Matrix& a, int& info) +{ + info = init (a); +} + +inline EIG::EIG (const ComplexMatrix& a) +{ + init (a); +} + +inline EIG::EIG (const ComplexMatrix& a, int& info) +{ + info = init (a); +} + +inline EIG::EIG (const EIG& a) +{ + lambda = a.lambda; + v = a.v; +} + +inline EIG& EIG::operator = (const EIG& a) +{ + lambda = a.lambda; + v = a.v; + return *this; +} + +inline ComplexColumnVector EIG::eigenvalues (void) const +{ + return lambda; +} + +inline ComplexMatrix EIG::eigenvectors (void) const +{ + return v; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/