Mercurial > jwe > octave
changeset 457:3d4b4f0fa5ba
[project @ 1994-06-06 00:33:33 by jwe]
Initial revision
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxAEPBAL.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,85 @@ +// -*- 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 "CmplxAEPBAL.h" +#include "dMatrix.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgebal) (const char*, const int*, Complex*, const int*, + int*, int*, double*, int*, long, long); + + int F77_FCN (zgebak) (const char*, const char*, const int*, const int*, + const int*, double*, const int*, Complex*, + const int*, int*, long, long); +} + +int +ComplexAEPBALANCE::init (const ComplexMatrix& a, const char *balance_job) +{ + + int n = a.cols (); + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *scale = new double [n]; + +// Copy matrix into local structure. + + balanced_mat = a; + + F77_FCN (zgebal) (balance_job, &n, balanced_mat.fortran_vec (), + &n, &ilo, &ihi, scale, &info, 1L, 1L); + +// Initialize balancing matrix to identity. + + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat (i, i) = 1.0; + + F77_FCN (zgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, + balancing_mat.fortran_vec (), &n, &info, 1L, 1L); + + delete [] scale; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxAEPBAL.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,100 @@ +// -*- 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_ComplexAEPBALANCE_h) +#define octave_ComplexAEPBALANCE_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CMatrix.h" + +extern "C++" { + +class ComplexAEPBALANCE +{ +friend class ComplexMatrix; + +public: + + ComplexAEPBALANCE (void) {} + ComplexAEPBALANCE (const ComplexMatrix& a, const char *balance_job); + ComplexAEPBALANCE (const ComplexAEPBALANCE& a); + ComplexAEPBALANCE& operator = (const ComplexAEPBALANCE& a); + ComplexMatrix balanced_matrix (void) const; + ComplexMatrix balancing_matrix (void) const; + + friend ostream& operator << (ostream& os, const ComplexAEPBALANCE& a); + +private: + + int init (const ComplexMatrix& a, const char * balance_job); + + ComplexMatrix balanced_mat; + ComplexMatrix balancing_mat; +}; + +inline ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexMatrix& a, + const char * balance_job) +{ + init (a, balance_job); +} + +inline ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexAEPBALANCE& a) +{ + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; +} + +inline ComplexAEPBALANCE& +ComplexAEPBALANCE::operator = (const ComplexAEPBALANCE& a) +{ + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; + + return *this; +} + +inline ComplexMatrix ComplexAEPBALANCE::balanced_matrix (void) const +{ + return balanced_mat; +} + +inline ComplexMatrix ComplexAEPBALANCE::balancing_matrix (void) const +{ + return balancing_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxCHOL.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,82 @@ +// -*- 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 "CmplxCHOL.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zpotrf) (const char*, const int*, Complex*, const int*, + int*, long); +} + +int +ComplexCHOL::init (const ComplexMatrix& a) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) + ("ComplexCHOL requires square matrix"); + return -1; + } + + char uplo = 'U'; + + int n = a_nc; + int info; + + Complex *h = dup (a.data (), a.length ()); + + F77_FCN (zpotrf) (&uplo, &n, h, &n, &info, 1L); + + chol_mat = ComplexMatrix (h, n, n); + +// If someone thinks of a more graceful way of doing this (or faster for +// that matter :-)), please let me know! + + if (n > 1) + for (int j = 0; j < a_nc; j++) + for (int i = j+1; i < a_nr; i++) + chol_mat.elem (i, j) = 0.0; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxCHOL.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,96 @@ +// -*- 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_ComplexCHOL_h) +#define octave_ComplexCHOL_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CMatrix.h" + +extern "C++" { + +class ComplexCHOL +{ +friend class ComplexMatrix; + +public: + + ComplexCHOL (void) {} + ComplexCHOL (const ComplexMatrix& a); + ComplexCHOL (const ComplexMatrix& a, int& info); + ComplexCHOL (const ComplexCHOL& a); + ComplexCHOL& operator = (const ComplexCHOL& a); + ComplexMatrix chol_matrix (void) const; + + friend ostream& operator << (ostream& os, const ComplexCHOL& a); + +private: + + int init (const ComplexMatrix& a); + + ComplexMatrix chol_mat; +}; + +inline ComplexCHOL::ComplexCHOL (const ComplexMatrix& a) +{ + init (a); +} + +inline ComplexCHOL::ComplexCHOL (const ComplexMatrix& a, int& info) +{ + info = init (a); +} + +inline ComplexCHOL::ComplexCHOL (const ComplexCHOL& a) +{ + chol_mat = a.chol_mat; +} + +inline ComplexCHOL& +ComplexCHOL::operator = (const ComplexCHOL& a) +{ + chol_mat = a.chol_mat; + + return *this; +} + +inline ComplexMatrix ComplexCHOL::chol_matrix (void) const +{ + return chol_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxDET.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,74 @@ +// -*- 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 <math.h> +#include <float.h> + +#include <Complex.h> + +#include "CmplxDET.h" + +int +ComplexDET::value_will_overflow (void) const +{ + return det[2].real () + 1 > log10 (DBL_MAX) ? 1 : 0; +} + +int +ComplexDET::value_will_underflow (void) const +{ + return det[2].real () - 1 < log10 (DBL_MIN) ? 1 : 0; +} + +Complex +ComplexDET::coefficient (void) const +{ + return det[0]; +} + +int +ComplexDET::exponent (void) const +{ + return (int) (det[1].real ()); +} + +Complex +ComplexDET::value (void) const +{ + return det[0] * pow (10.0, det[1].real ()); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxDET.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,94 @@ +// -*- 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_ComplexDET_h) +#define octave_ComplexDET_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include <Complex.h> + +extern "C++" { + +class ComplexDET +{ +public: + + ComplexDET (void); + + ComplexDET (const ComplexDET& a); + + ComplexDET& operator = (const ComplexDET& a); + + int value_will_overflow (void) const; + int value_will_underflow (void) const; + Complex coefficient (void) const; + int exponent (void) const; + Complex value (void) const; + + friend ostream& operator << (ostream& os, const ComplexDET& a); + +private: + + ComplexDET (const Complex *d); + + Complex det [2]; +}; + +inline ComplexDET::ComplexDET (void) +{ +} + +inline ComplexDET::ComplexDET (const ComplexDET& a) +{ + det[0] = a.det[0]; + det[1] = a.det[1]; +} + +inline ComplexDET& ComplexDET::operator = (const ComplexDET& a) +{ + det[0] = a.det[0]; + det[1] = a.det[1]; + return *this; +} + +inline ComplexDET::ComplexDET (const Complex *d) +{ + det[0] = d[0]; + det[1] = d[1]; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxHESS.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,119 @@ +// -*- 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 "CmplxHESS.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgebal) (const char*, const int*, Complex*, const int*, + int*, int*, double*, int*, long, long); + + int F77_FCN (zgebak) (const char*, const char*, const int*, const int*, + const int*, double*, const int*, Complex*, + const int*, int*, long, long); + + int F77_FCN (zgehrd) (const int*, const int*, const int*, Complex*, + const int*, Complex*, Complex*, const int*, + int*, long, long); + + int F77_FCN (zunghr) (const int*, const int*, const int*, Complex*, + const int*, Complex*, Complex*, const int*, + int*, long, long); +} + +int +ComplexHESS::init (const ComplexMatrix& a) +{ + 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; + } + + char job = 'N'; + char side = 'R'; + + int n = a_nc; + int lwork = 32 * n; + int info; + int ilo; + int ihi; + + Complex *h = dup (a.data (), a.length ()); + + double *scale = new double [n]; + Complex *tau = new Complex [n-1]; + Complex *work = new Complex [lwork]; + Complex *z = new Complex [n*n]; + + F77_FCN (zgebal) (&job, &n, h, &n, &ilo, &ihi, scale, &info, 1L, 1L); + + F77_FCN (zgehrd) (&n, &ilo, &ihi, h, &n, tau, work, &lwork, &info, 1L, + 1L); + + copy (z, h, n*n); + + F77_FCN (zunghr) (&n, &ilo, &ihi, z, &n, tau, work, &lwork, &info, 1L, + 1L); + + F77_FCN (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); + +// 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 [] work; + delete [] tau; + delete [] scale; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxHESS.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,105 @@ +// -*- 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_ComplexHESS_h) +#define octave_ComplexHESS_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CMatrix.h" + +extern "C++" { + +class ComplexHESS +{ +friend class ComplexMatrix; + +public: + + ComplexHESS (void) {} + ComplexHESS (const ComplexMatrix& a); + ComplexHESS (const ComplexMatrix& a, int& info); + ComplexHESS (const ComplexHESS& a); + ComplexHESS& operator = (const ComplexHESS& a); + ComplexMatrix hess_matrix (void) const; + ComplexMatrix unitary_hess_matrix (void) const; + + friend ostream& operator << (ostream& os, const ComplexHESS& a); + +private: + + int init (const ComplexMatrix& a); + + ComplexMatrix hess_mat; + ComplexMatrix unitary_hess_mat; +}; + +inline ComplexHESS::ComplexHESS (const ComplexMatrix& a) +{ + init (a); +} + +inline ComplexHESS::ComplexHESS (const ComplexMatrix& a, int& info) +{ + info = init (a); +} + +inline ComplexHESS::ComplexHESS (const ComplexHESS& a) +{ + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; +} + +inline ComplexHESS& +ComplexHESS::operator = (const ComplexHESS& a) +{ + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; + + return *this; +} + +inline ComplexMatrix ComplexHESS::hess_matrix (void) const +{ + return hess_mat; +} + +inline ComplexMatrix ComplexHESS::unitary_hess_matrix (void) const +{ + return unitary_hess_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxLU.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,112 @@ +// -*- 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 "CmplxLU.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgesv) (const int*, const int*, Complex*, const int*, + int*, Complex*, const int*, int*); +} + +ComplexLU::ComplexLU (const ComplexMatrix& a) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) + { + (*current_liboctave_error_handler) ("ComplexLU requires square matrix"); + return; + } + + int n = a_nr; + + int *ipvt = new int [n]; + int *pvt = new int [n]; + Complex *tmp_data = dup (a.data (), a.length ()); + int info = 0; + int zero = 0; + Complex b; + + F77_FCN (zgesv) (&n, &zero, tmp_data, &n, ipvt, &b, &n, &info); + + ComplexMatrix A_fact (tmp_data, n, n); + + int i; + + for (i = 0; i < n; i++) + { + ipvt[i] -= 1; + pvt[i] = i; + } + + for (i = 0; i < n - 1; i++) + { + int k = ipvt[i]; + if (k != i) + { + int tmp = pvt[k]; + pvt[k] = pvt[i]; + pvt[i] = tmp; + } + } + + l.resize (n, n, 0.0); + u.resize (n, n, 0.0); + p.resize (n, n, 0.0); + + for (i = 0; i < n; i++) + { + p.elem (i, pvt[i]) = 1.0; + + int j; + + l.elem (i, i) = 1.0; + for (j = 0; j < i; j++) + l.elem (i, j) = A_fact.elem (i, j); + + for (j = i; j < n; j++) + u.elem (i, j) = A_fact.elem (i, j); + } + + delete [] ipvt; + delete [] pvt; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxLU.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,104 @@ +// -*- 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_ComplexLU_h) +#define octave_Complex_LU_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" +#include "CMatrix.h" + +extern "C++" { + +class ComplexLU +{ +friend class ComplexMatrix; + +public: + + ComplexLU (void) {} + + ComplexLU (const ComplexMatrix& a); + + ComplexLU (const ComplexLU& a); + + ComplexLU& operator = (const ComplexLU& a); + + ComplexMatrix L (void) const; + ComplexMatrix U (void) const; + Matrix P (void) const; + + friend ostream& operator << (ostream& os, const ComplexLU& a); + +private: + + ComplexMatrix l; + ComplexMatrix u; + Matrix p; +}; + +inline ComplexLU::ComplexLU (const ComplexLU& a) +{ + l = a.l; + u = a.u; + p = a.p; +} + +inline ComplexLU& ComplexLU::operator = (const ComplexLU& a) +{ + l = a.l; + u = a.u; + p = a.p; + return *this; +} + +inline ComplexMatrix ComplexLU::L (void) const +{ + return l; +} + +inline ComplexMatrix ComplexLU::U (void) const +{ + return u; +} + +inline Matrix ComplexLU::P (void) const +{ + return p; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxQR.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,101 @@ +// -*- 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 "CmplxQR.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgeqrf) (const int*, const int*, Complex*, const int*, + Complex*, Complex*, const int*, int*); + + int F77_FCN (zungqr) (const int*, const int*, const int*, Complex*, + const int*, Complex*, Complex*, const int*, int*); +} + +ComplexQR::ComplexQR (const ComplexMatrix& a) +{ + int m = a.rows (); + int n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) + ("ComplexQR must have non-empty matrix"); + return; + } + + Complex *tmp_data; + int min_mn = m < n ? m : n; + Complex *tau = new Complex[min_mn]; + int lwork = 32*n; + Complex *work = new Complex[lwork]; + int info = 0; + + if (m > n) + { + tmp_data = new Complex [m*m]; + copy (tmp_data, a.data (), a.length ()); + } + else + tmp_data = dup (a.data (), a.length ()); + + F77_FCN (zgeqrf) (&m, &n, tmp_data, &m, tau, work, &lwork, &info); + + delete [] work; + + r.resize (m, n, 0.0); + for (int j = 0; j < n; j++) + { + int limit = j < min_mn-1 ? j : min_mn-1; + for (int i = 0; i <= limit; i++) + r.elem (i, j) = tmp_data[m*j+i]; + } + + lwork = 32*m; + work = new Complex[lwork]; + + F77_FCN (zungqr) (&m, &m, &min_mn, tmp_data, &m, tau, work, &lwork, &info); + + q = ComplexMatrix (tmp_data, m, m); + + delete [] tau; + delete [] work; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxQR.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,92 @@ +// -*- 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_ComplexQR_h) +#define octave_ComplexQR_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CMatrix.h" + +extern "C++" { + +class ComplexQR +{ +public: + + ComplexQR (void) {} + + ComplexQR (const ComplexMatrix& A); + + ComplexQR (const ComplexQR& a); + + ComplexQR& operator = (const ComplexQR& a); + + ComplexMatrix Q (void) const; + ComplexMatrix R (void) const; + + friend ostream& operator << (ostream& os, const ComplexQR& a); + +private: + + ComplexMatrix q; + ComplexMatrix r; +}; + +inline ComplexQR::ComplexQR (const ComplexQR& a) +{ + q = a.q; + r = a.r; +} + +inline ComplexQR& ComplexQR::operator = (const ComplexQR& a) +{ + q = a.q; + r = a.r; + return *this; +} + +inline ComplexMatrix ComplexQR::Q (void) const +{ + return q; +} + +inline ComplexMatrix ComplexQR::R (void) const +{ + return r; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxSCHUR.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,135 @@ +// -*- 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 "CmplxSCHUR.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgeesx) (const char*, const char*, int (*)(Complex*), + const char*, const int*, Complex*, const int*, + int*, Complex*, Complex*, const int*, double*, + double*, Complex*, const int*, double*, int*, + int*, long, long); +} + +static int +complex_select_ana (Complex *a) +{ + return a->real () < 0.0; +} + +static int +complex_select_dig (Complex *a) +{ + return (abs (*a) < 1.0); +} + +int +ComplexSCHUR::init (const ComplexMatrix& a, const char *ord) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) + ("ComplexSCHUR requires square matrix"); + return -1; + } + + char jobvs = 'V'; + char sort; + if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd') + sort = 'S'; + else + sort = 'N'; + + char sense = 'N'; + + int n = a_nc; + int lwork = 8 * n; + int info; + int sdim; + double rconde; + double rcondv; + + double *rwork = new double [n]; + +// bwork is not referenced for non-ordered Schur. + + int *bwork = (int *) NULL; + if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd') + bwork = new int [n]; + + Complex *s = dup (a.data (), a.length ()); + + Complex *work = new Complex [lwork]; + Complex *q = new Complex [n*n]; + Complex *w = new Complex [n]; + + if (*ord == 'A' || *ord == 'a') + { + F77_FCN (zgeesx) (&jobvs, &sort, complex_select_ana, &sense, + &n, s, &n, &sdim, w, q, &n, &rconde, &rcondv, + work, &lwork, rwork, bwork, &info, 1L, 1L); + } + else if (*ord == 'D' || *ord == 'd') + { + F77_FCN (zgeesx) (&jobvs, &sort, complex_select_dig, &sense, + &n, s, &n, &sdim, w, q, &n, &rconde, &rcondv, + work, &lwork, rwork, bwork, &info, 1L, 1L); + } + else + { + F77_FCN (zgeesx) (&jobvs, &sort, (void *) 0, &sense, &n, s, + &n, &sdim, w, q, &n, &rconde, &rcondv, work, + &lwork, rwork, bwork, &info, 1L, 1L); + } + + schur_mat = ComplexMatrix (s,n,n); + unitary_mat = ComplexMatrix (q,n,n); + + delete [] w; + delete [] work; + delete [] rwork; + delete [] bwork; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxSCHUR.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,110 @@ +// -*- 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_ComplexSCHUR_h) +#define octave_ComplexSCHUR_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CMatrix.h" + +extern "C++" { + +class ComplexSCHUR +{ +friend class ComplexMatrix; + +public: + + ComplexSCHUR (void) {} + + ComplexSCHUR (const ComplexMatrix& a, const char *ord); + ComplexSCHUR (const ComplexMatrix& a, const char *ord, int& info); + + ComplexSCHUR (const ComplexSCHUR& a, const char *ord); + + ComplexSCHUR& operator = (const ComplexSCHUR& a); + + ComplexMatrix schur_matrix (void) const; + ComplexMatrix unitary_matrix (void) const; + + friend ostream& operator << (ostream& os, const ComplexSCHUR& a); + +private: + + int init (const ComplexMatrix& a, const char *ord); + + ComplexMatrix schur_mat; + ComplexMatrix unitary_mat; +}; + +inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord) +{ + init (a,ord); +} + +inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord, + int& info) +{ + info = init (a,ord); +} + +inline ComplexSCHUR::ComplexSCHUR (const ComplexSCHUR& a, const char *ord) +{ + schur_mat = a.schur_mat; + unitary_mat = a.unitary_mat; +} + +inline ComplexSCHUR& +ComplexSCHUR::operator = (const ComplexSCHUR& a) +{ + schur_mat = a.schur_mat; + unitary_mat = a.unitary_mat; + + return *this; +} + +inline ComplexMatrix ComplexSCHUR::schur_matrix (void) const +{ + return schur_mat; +} + +inline ComplexMatrix ComplexSCHUR::unitary_matrix (void) const +{ + return unitary_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxSVD.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,89 @@ +// -*- 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 "CmplxSVD.h" +#include "mx-inlines.cc" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgesvd) (const char*, const char*, const int*, + const int*, Complex*, const int*, double*, + Complex*, const int*, Complex*, const int*, + Complex*, const int*, double*, int*, long, long); +} + +int +ComplexSVD::init (const ComplexMatrix& a) +{ + int info; + + int m = a.rows (); + int n = a.cols (); + + char jobu = 'A'; + char jobv = 'A'; + + Complex *tmp_data = dup (a.data (), a.length ()); + + int min_mn = m < n ? m : n; + int max_mn = m > n ? m : n; + + Complex *u = new Complex[m*m]; + double *s_vec = new double[min_mn]; + Complex *vt = new Complex[n*n]; + + int lwork = 2*min_mn + max_mn; + Complex *work = new Complex[lwork]; + + int lrwork = 5*max_mn; + double *rwork = new double[lrwork]; + + F77_FCN (zgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m, + vt, &n, work, &lwork, rwork, &info, 1L, 1L); + + left_sm = ComplexMatrix (u, m, m); + sigma = DiagMatrix (s_vec, m, n); + ComplexMatrix vt_m (vt, n, n); + right_sm = ComplexMatrix (vt_m.hermitian ()); + + delete [] tmp_data; + delete [] work; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxSVD.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,119 @@ +// -*- 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_ComplexSVD_h) +#define octave_ComplexSVD_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dDiagMatrix.h" +#include "CMatrix.h" + +extern "C++" { + +class ComplexSVD +{ +friend class ComplexMatrix; + +public: + + ComplexSVD (void) {} + + ComplexSVD (const ComplexMatrix& a); + ComplexSVD (const ComplexMatrix& a, int& info); + + ComplexSVD (const ComplexSVD& a); + + ComplexSVD& operator = (const ComplexSVD& a); + + DiagMatrix singular_values (void) const; + ComplexMatrix left_singular_matrix (void) const; + ComplexMatrix right_singular_matrix (void) const; + + friend ostream& operator << (ostream& os, const ComplexSVD& a); + +private: + + int init (const ComplexMatrix& a); + + DiagMatrix sigma; + ComplexMatrix left_sm; + ComplexMatrix right_sm; +}; + +inline ComplexSVD::ComplexSVD (const ComplexMatrix& a) +{ + init (a); +} + +inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, int& info) +{ + info = init (a); +} + +inline ComplexSVD::ComplexSVD (const ComplexSVD& a) +{ + sigma = a.sigma; + left_sm = a.left_sm; + right_sm = a.right_sm; +} + +inline ComplexSVD& +ComplexSVD::operator = (const ComplexSVD& a) +{ + sigma = a.sigma; + left_sm = a.left_sm; + right_sm = a.right_sm; + + return *this; +} + +inline DiagMatrix ComplexSVD::singular_values (void) const +{ + return sigma; +} + +inline ComplexMatrix ComplexSVD::left_singular_matrix (void) const +{ + return left_sm; +} + +inline ComplexMatrix ComplexSVD::right_singular_matrix (void) const +{ + return right_sm; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleAEPBAL.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,91 @@ +// -*- 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 "dbleAEPBAL.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgebal) (const char*, const int*, double*, + const int*, int*, int*, double*, + int*, long, long); + + int F77_FCN (dgebak) (const char*, const char*, const int*, const int*, + const int*, double*, const int*, double*, const int*, + int*, long, long); +} + +int +AEPBALANCE::init (const Matrix& a, const char *balance_job) +{ + int a_nc = a.cols (); + if (a.rows () != a_nc) + { + (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); + return -1; + } + + int n = a_nc; + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *scale = new double [n]; + +// Copy matrix into local structure. + + balanced_mat = a; + + F77_FCN (dgebal) (balance_job, &n, balanced_mat.fortran_vec (), + &n, &ilo, &ihi, scale, &info, 1L, 1L); + +// Initialize balancing matrix to identity. + + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; + + F77_FCN (dgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, + balancing_mat.fortran_vec (), &n, &info, 1L, 1L); + + delete [] scale; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleAEPBAL.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,101 @@ +// -*- 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_AEPBALANCE_h) +#define octave_AEPBALANCE_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class AEPBALANCE +{ +friend class Matrix; + +public: + + AEPBALANCE (void) {} + + AEPBALANCE (const Matrix& a, const char *balance_job); + + AEPBALANCE (const AEPBALANCE& a); + + AEPBALANCE& operator = (const AEPBALANCE& a); + Matrix balanced_matrix (void) const; + Matrix balancing_matrix (void) const; + friend ostream& operator << (ostream& os, const AEPBALANCE& a); + +private: + + int init (const Matrix& a, const char * balance_job); + + Matrix balanced_mat; + Matrix balancing_mat; +}; + +inline AEPBALANCE::AEPBALANCE (const Matrix& a,const char * balance_job) +{ + init (a, balance_job); +} + +inline AEPBALANCE::AEPBALANCE (const AEPBALANCE& a) +{ + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; +} + +inline AEPBALANCE& +AEPBALANCE::operator = (const AEPBALANCE& a) +{ + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; + + return *this; +} + +inline Matrix AEPBALANCE::balanced_matrix (void) const +{ + return balanced_mat; +} + +inline Matrix AEPBALANCE::balancing_matrix (void) const +{ + return balancing_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleCHOL.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,81 @@ +// -*- 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 "dbleCHOL.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dpotrf) (const char*, const int*, double*, const int*, + int*, long); +} + +int +CHOL::init (const Matrix& a) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) ("CHOL requires square matrix"); + return -1; + } + + char uplo = 'U'; + + int n = a_nc; + int info; + + double *h = dup (a.data (), a.length ()); + + F77_FCN (dpotrf) (&uplo, &n, h, &n, &info, 1L); + + chol_mat = Matrix (h, n, n); + +// If someone thinks of a more graceful way of doing this (or faster for +// that matter :-)), please let me know! + + if (n > 1) + for (int j = 0; j < a_nc; j++) + for (int i = j+1; i < a_nr; i++) + chol_mat.elem (i, j) = 0.0; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleCHOL.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,98 @@ +// -*- 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_CHOL_h) +#define octave_CHOL_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class CHOL +{ +friend class Matrix; + +public: + + CHOL (void) {} + + CHOL (const Matrix& a); + CHOL (const Matrix& a, int& info); + + CHOL (const CHOL& a); + + CHOL& operator = (const CHOL& a); + Matrix chol_matrix (void) const; + friend ostream& operator << (ostream& os, const CHOL& a); + +private: + + int init (const Matrix& a); + + Matrix chol_mat; +}; + +inline CHOL::CHOL (const Matrix& a) +{ + init (a); +} + +inline CHOL::CHOL (const Matrix& a, int& info) +{ + info = init (a); +} + +inline CHOL::CHOL (const CHOL& a) +{ + chol_mat = a.chol_mat; +} + +inline CHOL& +CHOL::operator = (const CHOL& a) +{ + chol_mat = a.chol_mat; + + return *this; +} + +inline Matrix CHOL::chol_matrix (void) const +{ + return chol_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleDET.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,72 @@ +// -*- 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 <math.h> +#include <float.h> + +#include "dbleDET.h" + +int +DET::value_will_overflow (void) const +{ + return det[2] + 1 > log10 (DBL_MAX) ? 1 : 0; +} + +int +DET::value_will_underflow (void) const +{ + return det[2] - 1 < log10 (DBL_MIN) ? 1 : 0; +} + +double +DET::coefficient (void) const +{ + return det[0]; +} + +int +DET::exponent (void) const +{ + return (int) det[1]; +} + +double +DET::value (void) const +{ + return det[0] * pow (10.0, det[1]); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleDET.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,92 @@ +// -*- 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_DET_h) +#define octave_DET_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +extern "C++" { + +class DET +{ +public: + + DET (void); + + DET (const DET& a); + + DET& operator = (const DET& a); + + int value_will_overflow (void) const; + int value_will_underflow (void) const; + double coefficient (void) const; + int exponent (void) const; + double value (void) const; + + friend ostream& operator << (ostream& os, const DET& a); + +private: + + DET (const double *d); + + double det [2]; +}; + +inline DET::DET (void) +{ +} + +inline DET::DET (const DET& a) +{ + det[0] = a.det[0]; + det[1] = a.det[1]; +} + +inline DET& DET::operator = (const DET& a) +{ + det[0] = a.det[0]; + det[1] = a.det[1]; + return *this; +} + +inline DET::DET (const double *d) +{ + det[0] = d[0]; + det[1] = d[1]; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleGEPBAL.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,190 @@ +// -*- 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 <math.h> + +#include "dbleGEPBAL.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgebak) (const char*, const char*, const int*, const int*, + const int*, double*, const int*, double*, const int*, + int*, long, long); + + int F77_FCN (reduce) (const int*, const int*, double*, + const int*, double*, + int*, int*, double*, double*); + + int F77_FCN (scaleg) (const int*, const int*, double*, + const int*, double*, + const int*, const int*, double*, double*, double*); + + int F77_FCN (gradeq) (const int*, const int*, double*, + const int*, double*, + int*, int*, double*, double*); +} + +int +GEPBALANCE::init (const Matrix& a, const Matrix& b, const char *balance_job) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + int b_nr = b.rows (); + if (a_nr != a_nc || a_nr != b_nr || b_nr != b.cols ()) + { + (*current_liboctave_error_handler) + ("GEPBALANCE requires square matrices of the same size"); + return -1; + } + + int n = a_nc; + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *cscale = new double [n]; + double *cperm = new double [n]; + Matrix wk (n, 6, 0.0); + +// Back out the permutations: +// +// cscale contains the exponents of the column scaling factors in its +// ilo through ihi locations and the reducing column permutations in +// its first ilo-1 and its ihi+1 through n locations. +// +// cperm contains the column permutations applied in grading the a and b +// submatrices in its ilo through ihi locations. +// +// wk contains the exponents of the row scaling factors in its ilo +// through ihi locations, the reducing row permutations in its first +// ilo-1 and its ihi+1 through n locations, and the row permutations +// applied in grading the a and b submatrices in its n+ilo through +// n+ihi locations. + +// Copy matrices into local structure. + + balanced_a_mat = a; + balanced_b_mat = b; + +// Initialize balancing matrices to identity. + + left_balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + left_balancing_mat (i, i) = 1.0; + + right_balancing_mat = left_balancing_mat; + +// Check for permutation option. + + if (*balance_job == 'P' || *balance_job == 'B') + { + F77_FCN (reduce) (&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cscale, wk.fortran_vec ()); + } + else + { + +// Set up for scaling later. + + ilo = 1; + ihi = n; + } + +// Check for scaling option. + + if ((*balance_job == 'S' || *balance_job == 'B') && ilo != ihi) + { + F77_FCN (scaleg) (&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cscale, cperm, wk.fortran_vec ()); + } + else + { + +// Set scaling data to 0's. + + for (int tmp = ilo-1; tmp < ihi; tmp++) + { + cscale[tmp] = 0.0; + wk.elem (tmp, 0) = 0.0; + } + } + +// Scaleg returns exponents, not values, so... + + for (int tmp = ilo-1; tmp < ihi; tmp++) + { + cscale[tmp] = pow (2.0, cscale[tmp]); + wk.elem (tmp, 0) = pow (2.0, -wk.elem (tmp, 0)); + } + +// Column permutations/scaling. + + F77_FCN (dgebak) (balance_job, "R", &n, &ilo, &ihi, cscale, &n, + right_balancing_mat.fortran_vec (), &n, &info, 1L, + 1L); + +// Row permutations/scaling. + + F77_FCN (dgebak) (balance_job, "L", &n, &ilo, &ihi, &wk.elem (0, 0), &n, + left_balancing_mat.fortran_vec (), &n, &info, 1L, 1L); + +// XXX FIXME XXX --- these four lines need to be added and debugged. +// GEPBALANCE::init will work without them, though, so here they are. + +#if 0 + if ((*balance_job == 'P' || *balance_job == 'B') && ilo != ihi) + { + F77_FCN (gradeq) (&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cperm, &wk.elem (0, 1)); + } +#endif + +// Transpose for aa = cc*a*dd convention... + left_balancing_mat = left_balancing_mat.transpose (); + + delete [] cscale; + delete [] cperm; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleGEPBAL.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,120 @@ +// -*- 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_GEPBALANCE_h) +#define octave_GEPBALANCE_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class GEPBALANCE +{ +friend class Matrix; + +public: + + GEPBALANCE (void) {} + + GEPBALANCE (const Matrix& a, const Matrix &, const char *balance_job); + + GEPBALANCE (const GEPBALANCE& a); + + GEPBALANCE& operator = (const GEPBALANCE& a); + Matrix balanced_a_matrix (void) const; + Matrix balanced_b_matrix (void) const; + Matrix left_balancing_matrix (void) const; + Matrix right_balancing_matrix (void) const; + friend ostream& operator << (ostream& os, const GEPBALANCE& a); + +private: + + int init (const Matrix& a, const Matrix& b, const char * balance_job); + + Matrix balanced_a_mat; + Matrix balanced_b_mat; + Matrix left_balancing_mat; + Matrix right_balancing_mat; +}; + +inline GEPBALANCE::GEPBALANCE (const Matrix& a, const Matrix& b, + const char * balance_job) +{ + init (a, b, balance_job); +} + +inline GEPBALANCE::GEPBALANCE (const GEPBALANCE& a) +{ + balanced_a_mat = a.balanced_a_mat; + balanced_b_mat = a.balanced_b_mat; + left_balancing_mat = a.left_balancing_mat; + right_balancing_mat = a.right_balancing_mat; +} + +inline GEPBALANCE& +GEPBALANCE::operator = (const GEPBALANCE& a) +{ + balanced_a_mat = a.balanced_a_mat; + balanced_b_mat = a.balanced_b_mat; + left_balancing_mat = a.left_balancing_mat; + right_balancing_mat = a.right_balancing_mat; + + return *this; +} + +inline Matrix GEPBALANCE::balanced_a_matrix (void) const +{ + return balanced_a_mat; +} + +inline Matrix GEPBALANCE::balanced_b_matrix (void) const +{ + return balanced_b_mat; +} + +inline Matrix GEPBALANCE::left_balancing_matrix (void) const +{ + return left_balancing_mat; +} + +inline Matrix GEPBALANCE::right_balancing_matrix (void) const +{ + return right_balancing_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleHESS.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,123 @@ +// -*- 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 "dbleHESS.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgebal) (const char*, const int*, double*, + const int*, int*, int*, double*, + int*, long, long); + + int F77_FCN (dgebak) (const char*, const char*, const int*, const int*, + const int*, double*, const int*, double*, const int*, + int*, long, long); + + int F77_FCN (dgehrd) (const int*, const int*, const int*, + double*, const int*, double*, double*, + const int*, int*, long, long); + + int F77_FCN (dorghr) (const int*, const int*, const int*, + double*, const int*, double*, double*, + const int*, int*, long, long); +} + +int +HESS::init (const Matrix& a) +{ + 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'; + + int n = a_nc; + int lwork = 32 * n; + int info; + int ilo; + int ihi; + + double *h = dup (a.data (), a.length ()); + + double *tau = new double [n+1]; + double *scale = new double [n]; + double *z = new double [n*n]; + double *work = new double [lwork]; + + F77_FCN (dgebal) (&jobbal, &n, h, &n, &ilo, &ihi, scale, &info, + 1L, 1L); + + F77_FCN (dgehrd) (&n, &ilo, &ihi, h, &n, tau, work, &lwork, &info, + 1L, 1L); + + copy (z, h, n*n); + + F77_FCN (dorghr) (&n, &ilo, &ihi, z, &n, tau, work, &lwork, &info, + 1L, 1L); + + F77_FCN (dgebak) (&jobbal, &side, &n, &ilo, &ihi, scale, &n, z, &n, + &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 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; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleHESS.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,107 @@ +// -*- 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_HESS_h) +#define octave_HESS_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class HESS +{ +friend class Matrix; + +public: + + HESS (void) {} + + HESS (const Matrix& a); + HESS (const Matrix&a, int& info); + + HESS (const HESS& a); + + HESS& operator = (const HESS& a); + Matrix hess_matrix (void) const; + Matrix unitary_hess_matrix (void) const; + friend ostream& operator << (ostream& os, const HESS& a); + +private: + + int init (const Matrix& a); + + Matrix hess_mat; + Matrix unitary_hess_mat; +}; + +inline HESS::HESS (const Matrix& a) +{ + init (a); +} + +inline HESS::HESS (const Matrix& a, int& info) +{ + info = init (a); +} + +inline HESS::HESS (const HESS& a) +{ + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; +} + +inline HESS& +HESS::operator = (const HESS& a) +{ + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; + + return *this; +} + +inline Matrix HESS::hess_matrix (void) const +{ + return hess_mat; +} + +inline Matrix HESS::unitary_hess_matrix (void) const +{ + return unitary_hess_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleLU.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,112 @@ +// -*- 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 "dbleLU.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgesv) (const int*, const int*, double*, const int*, + int*, double*, const int*, int*); +} + +LU::LU (const Matrix& a) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) + { + (*current_liboctave_error_handler) ("LU requires square matrix"); + return; + } + + int n = a_nr; + + int *ipvt = new int [n]; + int *pvt = new int [n]; + double *tmp_data = dup (a.data (), a.length ()); + int info = 0; + int zero = 0; + double b; + + F77_FCN (dgesv) (&n, &zero, tmp_data, &n, ipvt, &b, &n, &info); + + Matrix A_fact (tmp_data, n, n); + + int i; + + for (i = 0; i < n; i++) + { + ipvt[i] -= 1; + pvt[i] = i; + } + + for (i = 0; i < n - 1; i++) + { + int k = ipvt[i]; + if (k != i) + { + int tmp = pvt[k]; + pvt[k] = pvt[i]; + pvt[i] = tmp; + } + } + + l.resize (n, n, 0.0); + u.resize (n, n, 0.0); + p.resize (n, n, 0.0); + + for (i = 0; i < n; i++) + { + p.elem (i, pvt[i]) = 1.0; + + int j; + + l.elem (i, i) = 1.0; + for (j = 0; j < i; j++) + l.elem (i, j) = A_fact.elem (i, j); + + for (j = i; j < n; j++) + u.elem (i, j) = A_fact.elem (i, j); + } + + delete [] ipvt; + delete [] pvt; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleLU.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,103 @@ +// -*- 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_LU_h) +#define octave_LU_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class LU +{ +friend class Matrix; + +public: + + LU (void) {} + + LU (const Matrix& a); + + LU (const LU& a); + + LU& operator = (const LU& a); + + Matrix L (void) const; + Matrix U (void) const; + Matrix P (void) const; + + friend ostream& operator << (ostream& os, const LU& a); + +private: + + Matrix l; + Matrix u; + Matrix p; +}; + +inline LU::LU (const LU& a) +{ + l = a.l; + u = a.u; + p = a.p; +} + +inline LU& LU::operator = (const LU& a) +{ + l = a.l; + u = a.u; + p = a.p; + return *this; +} + +inline Matrix LU::L (void) const +{ + return l; +} + +inline Matrix LU::U (void) const +{ + return u; +} + +inline Matrix LU::P (void) const +{ + return p; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleQR.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,100 @@ +// -*- 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 "dbleQR.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgeqrf) (const int*, const int*, double*, const int*, + double*, double*, const int*, int*); + + int F77_FCN (dorgqr) (const int*, const int*, const int*, double*, + const int*, double*, double*, const int*, int*); +} + +QR::QR (const Matrix& a) +{ + int m = a.rows (); + int n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) ("QR must have non-empty matrix"); + return; + } + + double *tmp_data; + int min_mn = m < n ? m : n; + double *tau = new double[min_mn]; + int lwork = 32*n; + double *work = new double[lwork]; + int info = 0; + + if (m > n) + { + tmp_data = new double [m*m]; + copy (tmp_data, a.data (), a.length ()); + } + else + tmp_data = dup (a.data (), a.length ()); + + F77_FCN (dgeqrf) (&m, &n, tmp_data, &m, tau, work, &lwork, &info); + + delete [] work; + + r.resize (m, n, 0.0); + for (int j = 0; j < n; j++) + { + int limit = j < min_mn-1 ? j : min_mn-1; + for (int i = 0; i <= limit; i++) + r.elem (i, j) = tmp_data[m*j+i]; + } + + lwork = 32*m; + work = new double[lwork]; + + F77_FCN (dorgqr) (&m, &m, &min_mn, tmp_data, &m, tau, work, &lwork, &info); + + q = Matrix (tmp_data, m, m); + + delete [] tau; + delete [] work; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleQR.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,92 @@ +// -*- 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_QR_h) +#define octave_QR_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class QR +{ +public: + + QR (void) {} + + QR (const Matrix& A); + + QR (const QR& a); + + QR& operator = (const QR& a); + + Matrix Q (void) const; + Matrix R (void) const; + + friend ostream& operator << (ostream& os, const QR& a); + +private: + + Matrix q; + Matrix r; +}; + +inline QR::QR (const QR& a) +{ + q = a.q; + r = a.r; +} + +inline QR& QR::operator = (const QR& a) +{ + q = a.q; + r = a.r; + return *this; +} + +inline Matrix QR::Q (void) const +{ + return q; +} + +inline Matrix QR::R (void) const +{ + return r; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleSCHUR.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,152 @@ +// -*- 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 "dbleSCHUR.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgeesx) (const char*, const char*, + int (*)(double*, double*), const char*, + const int*, double*, const int*, int*, double*, + double*, double*, const int*, double*, double*, + double*, const int*, int*, const int*, int*, + int*, long, long); +} + +static int +select_ana (double *a, double *b) +{ + return (*a < 0.0); +} + +static int +select_dig (double *a, double *b) +{ + return (hypot (*a, *b) < 1.0); +} + +int +SCHUR::init (const Matrix& a, const char *ord) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) ("SCHUR requires square matrix"); + return -1; + } + + char jobvs = 'V'; + char sort; + + if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd') + sort = 'S'; + else + sort = 'N'; + + char sense = 'N'; + + int n = a_nc; + int lwork = 8 * n; + int liwork = 1; + int info; + int sdim; + double rconde; + double rcondv; + + double *s = dup (a.data (), a.length ()); + + double *wr = new double [n]; + double *wi = new double [n]; + double *q = new double [n*n]; + double *work = new double [lwork]; + +// These are not referenced for the non-ordered Schur routine. + + int *iwork = (int *) NULL; + int *bwork = (int *) NULL; + if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd') + { + iwork = new int [liwork]; + bwork = new int [n]; + } + + if (*ord == 'A' || *ord == 'a') + { + F77_FCN (dgeesx) (&jobvs, &sort, select_ana, &sense, &n, s, &n, + &sdim, wr, wi, q, &n, &rconde, &rcondv, work, + &lwork, iwork, &liwork, bwork, &info, 1L, 1L); + } + else if (*ord == 'D' || *ord == 'd') + { + F77_FCN (dgeesx) (&jobvs, &sort, select_dig, &sense, &n, s, &n, + &sdim, wr, wi, q, &n, &rconde, &rcondv, work, + &lwork, iwork, &liwork, bwork, &info, 1L, 1L); + + } + else + { + F77_FCN (dgeesx) (&jobvs, &sort, (void *) 0, &sense, &n, s, + &n, &sdim, wr, wi, q, &n, &rconde, &rcondv, + work, &lwork, iwork, &liwork, bwork, &info, + 1L, 1L); + } + + schur_mat = Matrix (s, n, n); + unitary_mat = Matrix (q, n, n); + + delete [] wr; + delete [] wi; + delete [] work; + delete [] iwork; + delete [] bwork; + + return info; +} + +ostream& +operator << (ostream& os, const SCHUR& a) +{ + os << a.schur_matrix () << "\n"; + os << a.unitary_matrix () << "\n"; + + return os; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleSCHUR.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,109 @@ +// -*- 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_SCHUR_h) +#define octave_SCHUR_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dMatrix.h" + +extern "C++" { + +class SCHUR +{ +friend class Matrix; + +public: + + SCHUR (void) {} + + SCHUR (const Matrix& a, const char *ord); + SCHUR (const Matrix& a, const char *ord, int& info); + + SCHUR (const SCHUR& a, const char *ord); + + SCHUR& operator = (const SCHUR& a); + + Matrix schur_matrix (void) const; + Matrix unitary_matrix (void) const; + + friend ostream& operator << (ostream& os, const SCHUR& a); + +private: + + int init (const Matrix& a, const char *ord); + + Matrix schur_mat; + Matrix unitary_mat; +}; + +inline SCHUR::SCHUR (const Matrix& a, const char *ord) +{ + init (a, ord); +} + +inline SCHUR::SCHUR (const Matrix& a, const char *ord, int& info) +{ + info = init (a, ord); +} + +inline SCHUR::SCHUR (const SCHUR& a, const char *ord) +{ + schur_mat = a.schur_mat; + unitary_mat = a.unitary_mat; +} + +inline SCHUR& +SCHUR::operator = (const SCHUR& a) +{ + schur_mat = a.schur_mat; + unitary_mat = a.unitary_mat; + + return *this; +} + +inline Matrix SCHUR::schur_matrix (void) const +{ + return schur_mat; +} + +inline Matrix SCHUR::unitary_matrix (void) const +{ + return unitary_mat; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleSVD.cc Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,98 @@ +// -*- 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 "dbleSVD.h" +#include "mx-inlines.cc" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgesvd) (const char*, const char*, const int*, + const int*, double*, const int*, double*, + double*, const int*, double*, const int*, + double*, const int*, int*, long, long); +} + +int +SVD::init (const Matrix& a) +{ + int info; + + int m = a.rows (); + int n = a.cols (); + + char jobu = 'A'; + char jobv = 'A'; + + double *tmp_data = dup (a.data (), a.length ()); + + int min_mn = m < n ? m : n; + int max_mn = m > n ? m : n; + + double *u = new double[m*m]; + double *s_vec = new double[min_mn]; + double *vt = new double[n*n]; + + int tmp1 = 3*min_mn + max_mn; + int tmp2 = 5*min_mn - 4; + int lwork = tmp1 > tmp2 ? tmp1 : tmp2; + double *work = new double[lwork]; + + F77_FCN (dgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m, + vt, &n, work, &lwork, &info, 1L, 1L); + + left_sm = Matrix (u, m, m); + sigma = DiagMatrix (s_vec, m, n); + Matrix vt_m (vt, n, n); + right_sm = Matrix (vt_m.transpose ()); + + delete [] tmp_data; + delete [] work; + + return info; +} + +ostream& +operator << (ostream& os, const SVD& a) +{ + os << a.left_singular_matrix () << "\n"; + os << a.singular_values () << "\n"; + os << a.right_singular_matrix () << "\n"; + + return os; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleSVD.h Mon Jun 06 00:33:51 1994 +0000 @@ -0,0 +1,119 @@ +// -*- 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_SVD_h) +#define octave_SVD_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dDiagMatrix.h" +#include "dMatrix.h" + +extern "C++" { + +class SVD +{ +friend class Matrix; + +public: + + SVD (void) {} + + SVD (const Matrix& a); + SVD (const Matrix& a, int& info); + + SVD (const SVD& a); + + SVD& operator = (const SVD& a); + + DiagMatrix singular_values (void) const; + Matrix left_singular_matrix (void) const; + Matrix right_singular_matrix (void) const; + + friend ostream& operator << (ostream& os, const SVD& a); + +private: + + int init (const Matrix& a); + + DiagMatrix sigma; + Matrix left_sm; + Matrix right_sm; +}; + +inline SVD::SVD (const Matrix& a) +{ + init (a); +} + +inline SVD::SVD (const Matrix& a, int& info) +{ + info = init (a); +} + +inline SVD::SVD (const SVD& a) +{ + sigma = a.sigma; + left_sm = a.left_sm; + right_sm = a.right_sm; +} + +inline SVD& +SVD::operator = (const SVD& a) +{ + sigma = a.sigma; + left_sm = a.left_sm; + right_sm = a.right_sm; + + return *this; +} + +inline DiagMatrix SVD::singular_values (void) const +{ + return sigma; +} + +inline Matrix SVD::left_singular_matrix (void) const +{ + return left_sm; +} + +inline Matrix SVD::right_singular_matrix (void) const +{ + return right_sm; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/