Mercurial > hg > octave-terminal
changeset 538:8e134d3b21c9
[project @ 1994-07-21 22:40:04 by jwe]
Initial revision
author | jwe |
---|---|
date | Thu, 21 Jul 1994 22:40:04 +0000 |
parents | 4ecbfd3c3710 |
children | 5ec10a984241 |
files | liboctave/CmplxQRP.cc liboctave/CmplxQRP.h liboctave/dbleQRP.cc liboctave/dbleQRP.h |
diffstat | 4 files changed, 459 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/liboctave/CmplxQRP.cc @@ -0,0 +1,151 @@ +// -*- 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 <assert.h> + +#include "CmplxQRP.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (zgeqpf) (const int*, const int*, Complex*, const int*, + int*, Complex*, Complex*, double*, int*); +} + +// It would be best to share some of this code with ComplexQR class... + +ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type) +{ + assert (qr_type != QR::raw); + + 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 = 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 ()); + + + work = new Complex[n]; + double *rwork = new double[2*n]; + int *jpvt = new int[n]; + +// Clear Pivot vector (code to enforce a certain permutation would go +// here...) + + for (int i = 0; i < n; i++) + jpvt[i] = 0; + + F77_FCN (zgeqpf) (&m, &n, tmp_data, &m, jpvt, tau, work, rwork, &info); + +// Form Permutation matrix (if economy is requested, return the +// indices only!) + + if (qr_type == QR::economy && m > n) + { + p.resize (1, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (0, j) = jpvt[j]; + } + else + { + p.resize (n, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (jpvt[j]-1, j) = 1.0; + } + + delete [] jpvt; + delete [] work; + + int m2; + if (qr_type == QR::economy && m > n) + { + m2 = n; + r.resize (n, n, 0.0); + } + else + { + m2 = m; + 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: *** +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/CmplxQRP.h @@ -0,0 +1,84 @@ +// -*- 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_ComplexQRP_h) +#define octave_ComplexQRP_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "CmplxQR.h" + +extern "C++" { + +class ComplexQRP : public ComplexQR +{ +public: + + ComplexQRP (void) {} + + ComplexQRP (const ComplexMatrix& A, QR::type qr_type = QR::std); + + ComplexQRP (const ComplexQRP& a); + + ComplexQRP& operator = (const ComplexQRP& a); + + Matrix P (void) const; + + friend ostream& operator << (ostream& os, const ComplexQRP& a); + +private: + + Matrix p; +}; + +inline ComplexQRP::ComplexQRP (const ComplexQRP& a) : ComplexQR (a) +{ + p = a.p; +} + +inline ComplexQRP& ComplexQRP::operator = (const ComplexQRP& a) +{ + ComplexQR::operator = (a); + p = a.p; + return *this; +} + +inline Matrix ComplexQRP::P (void) const +{ + return p; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/dbleQRP.cc @@ -0,0 +1,140 @@ +// -*- 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 <assert.h> + +#include "dbleQRP.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +extern "C" +{ + int F77_FCN (dgeqpf) (const int*, const int*, double*, const int*, + int*, double*, double*, int*); +} + +// It would be best to share some of this code with QR class... + +QRP::QRP (const Matrix& a, QR::type qr_type) +{ + assert (qr_type != QR::raw); + + 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 = 3*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 ()); + + int *jpvt = new int[n]; + +// Clear Pivot vector (code to enforce a certain permutation would go +// here...) + + for (int i = 0; i < n; i++) + jpvt[i] = 0; + + F77_FCN (dgeqpf) (&m, &n, tmp_data, &m, jpvt, tau, work, &info); + +// Form Permutation matrix (if economy is requested, return the +// indices only!) + + if (qr_type == QR::economy && m > n) + { + p.resize (1, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (0, j) = jpvt[j]; + } + else + { + p.resize (n, n, 0.0); + for (int j = 0; j < n; j++) + p.elem (jpvt[j]-1, j) = 1.0; + } + + delete [] jpvt; + delete [] work; + + int m2; + if (qr_type == QR::economy && m > n) + { + m2 = n; + r.resize (n, n, 0.0); + } + else + { + m2 = m; + 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: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/dbleQRP.h @@ -0,0 +1,84 @@ +// -*- 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_QRP_h) +#define octave_QRP_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ostream; + +#include "dbleQR.h" + +extern "C++" { + +class QRP : public QR +{ +public: + + QRP (void) {} + + QRP (const Matrix& A, QR::type qr_type = QR::std); + + QRP (const QRP& a); + + QRP& operator = (const QRP& a); + + Matrix P (void) const; + + friend ostream& operator << (ostream& os, const QRP& a); + +private: + + Matrix p; +}; + +inline QRP::QRP (const QRP& a) : QR (a) +{ + p = a.p; +} + +inline QRP& QRP::operator = (const QRP& a) +{ + QR::operator = (a); + p = a.p; + return *this; +} + +inline Matrix QRP::P (void) const +{ + return p; +} + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/