Mercurial > octave
changeset 21280:ebdf74c15722
better use of templates for qrp classes
* liboctave/numeric/qrp.h, liboctave/numeric/qrp.cc: New files for qrp
classes generated from CmplxQRP.cc, CmplxQRP.h, dbleQRP.cc, dbleQRP.h,
fCmplxQRP.cc, fCmplxQRP.h, floatQRP.cc, and floatQRP.h with classes
converted to templates.
* liboctave/numeric/module.mk: Update.
* qr.cc, mx-defs.h: Use new classes.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Wed, 17 Feb 2016 02:54:00 -0500 |
| parents | eb1524b07fe3 |
| children | b76955e83fe4 |
| files | libinterp/dldfcn/qr.cc liboctave/numeric/CmplxQRP.cc liboctave/numeric/CmplxQRP.h liboctave/numeric/dbleQRP.cc liboctave/numeric/dbleQRP.h liboctave/numeric/fCmplxQRP.cc liboctave/numeric/fCmplxQRP.h liboctave/numeric/floatQRP.cc liboctave/numeric/floatQRP.h liboctave/numeric/module.mk liboctave/numeric/qrp.cc liboctave/numeric/qrp.h liboctave/operators/mx-defs.h |
| diffstat | 13 files changed, 411 insertions(+), 751 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/qr.cc Wed Feb 17 02:57:21 2016 -0500 +++ b/libinterp/dldfcn/qr.cc Wed Feb 17 02:54:00 2016 -0500 @@ -26,11 +26,8 @@ # include <config.h> #endif -#include "CmplxQRP.h" -#include "dbleQRP.h" -#include "fCmplxQRP.h" -#include "floatQRP.h" #include "qr.h" +#include "qrp.h" #include "sparse-qr.h" @@ -304,7 +301,7 @@ default: { - FloatQRP fact (m, type); + qrp<FloatMatrix> fact (m, type); if (type == qr<FloatMatrix>::economy) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); @@ -340,7 +337,7 @@ default: { - FloatComplexQRP fact (m, type); + qrp<FloatComplexMatrix> fact (m, type); if (type == qr<FloatComplexMatrix>::economy) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else @@ -377,7 +374,7 @@ default: { - QRP fact (m, type); + qrp<Matrix> fact (m, type); if (type == qr<Matrix>::economy) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else @@ -412,7 +409,7 @@ default: { - ComplexQRP fact (m, type); + qrp<ComplexMatrix> fact (m, type); if (type == qr<ComplexMatrix>::economy) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else
--- a/liboctave/numeric/CmplxQRP.cc Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,109 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague - -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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <cassert> - -#include "CmplxQRP.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (zgeqp3, ZGEQP3) (const octave_idx_type&, const octave_idx_type&, - Complex*, const octave_idx_type&, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, double*, - octave_idx_type&); -} - -// It would be best to share some of this code with qr<ComplexMatrix> class... - -ComplexQRP::ComplexQRP (const ComplexMatrix& a, type qr_type) - : qr<ComplexMatrix> (), p () -{ - init (a, qr_type); -} - -void -ComplexQRP::init (const ComplexMatrix& a, type qr_type) -{ - assert (qr_type != qr<ComplexMatrix>::raw); - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - octave_idx_type min_mn = m < n ? m : n; - OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn); - - octave_idx_type info = 0; - - ComplexMatrix afact = a; - if (m > n && qr_type == qr<ComplexMatrix>::std) - afact.resize (m, m); - - MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); - - if (m > 0) - { - OCTAVE_LOCAL_BUFFER (double, rwork, 2*n); - - // workspace query. - Complex clwork; - F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - &clwork, -1, rwork, info)); - - // allocate buffer and do the job. - octave_idx_type lwork = clwork.real (); - lwork = std::max (lwork, static_cast<octave_idx_type> (1)); - OCTAVE_LOCAL_BUFFER (Complex, work, lwork); - F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - work, lwork, rwork, info)); - } - else - for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; - - // Form Permutation matrix (if economy is requested, return the - // indices only!) - - jpvt -= static_cast<octave_idx_type> (1); - p = PermMatrix (jpvt, true); - - - form (n, afact, tau, qr_type); -} - -RowVector -ComplexQRP::Pvec (void) const -{ - Array<double> pa (p.col_perm_vec ()); - RowVector pv (MArray<double> (pa) + 1.0); - return pv; -}
--- a/liboctave/numeric/CmplxQRP.h Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -/* - -Copyright (C) 1994-2015 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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_CmplxQRP_h) -#define octave_CmplxQRP_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "PermMatrix.h" -#include "CColVector.h" -#include "CMatrix.h" -#include "CRowVector.h" -#include "qr.h" - -class -OCTAVE_API -ComplexQRP : public qr<ComplexMatrix> -{ -public: - - typedef qr<ComplexMatrix>::type type; - - ComplexQRP (void) : qr<ComplexMatrix> (), p () { } - - ComplexQRP (const ComplexMatrix&, type = qr<ComplexMatrix>::std); - - ComplexQRP (const ComplexQRP& a) : qr<ComplexMatrix> (a), p (a.p) { } - - ComplexQRP& operator = (const ComplexQRP& a) - { - if (this != &a) - { - qr<ComplexMatrix>::operator = (a); - p = a.p; - } - return *this; - } - - ~ComplexQRP (void) { } - - void init (const ComplexMatrix&, type = qr<ComplexMatrix>::std); - - PermMatrix P (void) const { return p; } - - RowVector Pvec (void) const; - - friend std::ostream& operator << (std::ostream&, const ComplexQRP&); - -private: - - PermMatrix p; -}; - -#endif
--- a/liboctave/numeric/dbleQRP.cc Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,106 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague - -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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <cassert> - -#include "dbleQRP.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, - double*, const octave_idx_type&, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type&); -} - -// It would be best to share some of this code with qr<Matrix> class... - -QRP::QRP (const Matrix& a, type qr_type) - : qr<Matrix> (), p () -{ - init (a, qr_type); -} - -void -QRP::init (const Matrix& a, type qr_type) -{ - assert (qr_type != qr<Matrix>::raw); - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - octave_idx_type min_mn = m < n ? m : n; - OCTAVE_LOCAL_BUFFER (double, tau, min_mn); - - octave_idx_type info = 0; - - Matrix afact = a; - if (m > n && qr_type == qr<Matrix>::std) - afact.resize (m, m); - - MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); - - if (m > 0) - { - // workspace query. - double rlwork; - F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - &rlwork, -1, info)); - - // allocate buffer and do the job. - octave_idx_type lwork = rlwork; - lwork = std::max (lwork, static_cast<octave_idx_type> (1)); - OCTAVE_LOCAL_BUFFER (double, work, lwork); - F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - work, lwork, info)); - } - else - for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; - - // Form Permutation matrix (if economy is requested, return the - // indices only!) - - jpvt -= static_cast<octave_idx_type> (1); - p = PermMatrix (jpvt, true); - - - form (n, afact, tau, qr_type); -} - -RowVector -QRP::Pvec (void) const -{ - Array<double> pa (p.col_perm_vec ()); - RowVector pv (MArray<double> (pa) + 1.0); - return pv; -}
--- a/liboctave/numeric/dbleQRP.h Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -/* - -Copyright (C) 1994-2015 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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_dbleQRP_h) -#define octave_dbleQRP_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "PermMatrix.h" -#include "dColVector.h" -#include "dMatrix.h" -#include "dRowVector.h" -#include "qr.h" - -class -OCTAVE_API -QRP : public qr<Matrix> -{ -public: - - typedef qr<Matrix>::type type; - - QRP (void) : qr<Matrix> (), p () { } - - QRP (const Matrix&, type = qr<Matrix>::std); - - QRP (const QRP& a) : qr<Matrix> (a), p (a.p) { } - - QRP& operator = (const QRP& a) - { - if (this != &a) - { - qr<Matrix>::operator = (a); - p = a.p; - } - - return *this; - } - - ~QRP (void) { } - - void init (const Matrix&, type = qr<Matrix>::std); - - PermMatrix P (void) const { return p; } - - RowVector Pvec (void) const; - - friend std::ostream& operator << (std::ostream&, const QRP&); - -protected: - - PermMatrix p; -}; - -#endif
--- a/liboctave/numeric/fCmplxQRP.cc Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,108 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague - -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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <cassert> - -#include "fCmplxQRP.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (cgeqp3, CGEQP3) (const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - octave_idx_type*, FloatComplex*, FloatComplex*, - const octave_idx_type&, float*, octave_idx_type&); -} - -// It would be best to share some of this code with qr<FloatComplexMatrix> class... - -FloatComplexQRP::FloatComplexQRP (const FloatComplexMatrix& a, type qr_type) - : qr<FloatComplexMatrix> (), p () -{ - init (a, qr_type); -} - -void -FloatComplexQRP::init (const FloatComplexMatrix& a, type qr_type) -{ - assert (qr_type != qr<FloatComplexMatrix>::raw); - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - octave_idx_type min_mn = m < n ? m : n; - OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn); - - octave_idx_type info = 0; - - FloatComplexMatrix afact = a; - if (m > n && qr_type == qr<FloatComplexMatrix>::std) - afact.resize (m, m); - - MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); - - if (m > 0) - { - OCTAVE_LOCAL_BUFFER (float, rwork, 2*n); - - // workspace query. - FloatComplex clwork; - F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - &clwork, -1, rwork, info)); - - // allocate buffer and do the job. - octave_idx_type lwork = clwork.real (); - lwork = std::max (lwork, static_cast<octave_idx_type> (1)); - OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork); - F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - work, lwork, rwork, info)); - } - else - for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; - - // Form Permutation matrix (if economy is requested, return the - // indices only!) - - jpvt -= static_cast<octave_idx_type> (1); - p = PermMatrix (jpvt, true); - - - form (n, afact, tau, qr_type); -} - -FloatRowVector -FloatComplexQRP::Pvec (void) const -{ - Array<float> pa (p.col_perm_vec ()); - FloatRowVector pv (MArray<float> (pa) + 1.0f); - return pv; -}
--- a/liboctave/numeric/fCmplxQRP.h Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -/* - -Copyright (C) 1994-2015 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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_fCmplxQRP_h) -#define octave_fCmplxQRP_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "PermMatrix.h" -#include "fCColVector.h" -#include "fCMatrix.h" -#include "fCRowVector.h" -#include "qr.h" - -class -OCTAVE_API -FloatComplexQRP : public qr<FloatComplexMatrix> -{ -public: - - typedef qr<FloatComplexMatrix>::type type; - - FloatComplexQRP (void) : qr<FloatComplexMatrix> (), p () { } - - FloatComplexQRP (const FloatComplexMatrix&, type = qr<FloatComplexMatrix>::std); - - FloatComplexQRP (const FloatComplexQRP& a) : qr<FloatComplexMatrix> (a), p (a.p) { } - - FloatComplexQRP& operator = (const FloatComplexQRP& a) - { - if (this != &a) - { - qr<FloatComplexMatrix>::operator = (a); - p = a.p; - } - return *this; - } - - ~FloatComplexQRP (void) { } - - void init (const FloatComplexMatrix&, type = qr<FloatComplexMatrix>::std); - - PermMatrix P (void) const { return p; } - - FloatRowVector Pvec (void) const; - - friend std::ostream& operator << (std::ostream&, const FloatComplexQRP&); - -private: - - PermMatrix p; -}; - -#endif
--- a/liboctave/numeric/floatQRP.cc Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,106 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague - -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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <cassert> - -#include "floatQRP.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (sgeqp3, SGEQP3) (const octave_idx_type&, const octave_idx_type&, - float*, const octave_idx_type&, octave_idx_type*, - float*, float*, const octave_idx_type&, - octave_idx_type&); -} - -// It would be best to share some of this code with QR class... - -FloatQRP::FloatQRP (const FloatMatrix& a, type qr_type) - : qr<FloatMatrix> (), p () -{ - init (a, qr_type); -} - -void -FloatQRP::init (const FloatMatrix& a, type qr_type) -{ - assert (qr_type != qr<FloatMatrix>::raw); - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - octave_idx_type min_mn = m < n ? m : n; - OCTAVE_LOCAL_BUFFER (float, tau, min_mn); - - octave_idx_type info = 0; - - FloatMatrix afact = a; - if (m > n && qr_type == qr<FloatMatrix>::std) - afact.resize (m, m); - - MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); - - if (m > 0) - { - // workspace query. - float rlwork; - F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - &rlwork, -1, info)); - - // allocate buffer and do the job. - octave_idx_type lwork = rlwork; - lwork = std::max (lwork, static_cast<octave_idx_type> (1)); - OCTAVE_LOCAL_BUFFER (float, work, lwork); - F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), - m, jpvt.fortran_vec (), tau, - work, lwork, info)); - } - else - for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; - - // Form Permutation matrix (if economy is requested, return the - // indices only!) - - jpvt -= static_cast<octave_idx_type> (1); - p = PermMatrix (jpvt, true); - - - form (n, afact, tau, qr_type); -} - -FloatRowVector -FloatQRP::Pvec (void) const -{ - Array<float> pa (p.col_perm_vec ()); - FloatRowVector pv (MArray<float> (pa) + 1.0f); - return pv; -}
--- a/liboctave/numeric/floatQRP.h Wed Feb 17 02:57:21 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -/* - -Copyright (C) 1994-2015 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 3 of the License, 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, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_floatQRP_h) -#define octave_floatQRP_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "PermMatrix.h" -#include "fColVector.h" -#include "fMatrix.h" -#include "fRowVector.h" -#include "qr.h" - -class -OCTAVE_API -FloatQRP : public qr<FloatMatrix> -{ -public: - - typedef qr<FloatMatrix>::type type; - - FloatQRP (void) : qr<FloatMatrix> (), p () { } - - FloatQRP (const FloatMatrix&, type = qr<FloatMatrix>::std); - - FloatQRP (const FloatQRP& a) : qr<FloatMatrix> (a), p (a.p) { } - - FloatQRP& operator = (const FloatQRP& a) - { - if (this != &a) - { - qr<FloatMatrix>::operator = (a); - p = a.p; - } - - return *this; - } - - ~FloatQRP (void) { } - - void init (const FloatMatrix&, type = qr<FloatMatrix>::std); - - PermMatrix P (void) const { return p; } - - FloatRowVector Pvec (void) const; - - friend std::ostream& operator << (std::ostream&, const FloatQRP&); - -protected: - - PermMatrix p; -}; - -#endif
--- a/liboctave/numeric/module.mk Wed Feb 17 02:57:21 2016 -0500 +++ b/liboctave/numeric/module.mk Wed Feb 17 02:54:00 2016 -0500 @@ -8,7 +8,6 @@ LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in) NUMERIC_INC = \ - liboctave/numeric/CmplxQRP.h \ liboctave/numeric/CollocWt.h \ liboctave/numeric/DAE.h \ liboctave/numeric/DAEFunc.h \ @@ -32,11 +31,8 @@ liboctave/numeric/bsxfun-decl.h \ liboctave/numeric/bsxfun.h \ liboctave/numeric/chol.h \ - liboctave/numeric/dbleQRP.h \ liboctave/numeric/eigs-base.h \ - liboctave/numeric/fCmplxQRP.h \ liboctave/numeric/fEIG.h \ - liboctave/numeric/floatQRP.h \ liboctave/numeric/gepbalance.h \ liboctave/numeric/hess.h \ liboctave/numeric/lo-mappers.h \ @@ -48,6 +44,7 @@ liboctave/numeric/oct-rand.h \ liboctave/numeric/oct-spparms.h \ liboctave/numeric/qr.h \ + liboctave/numeric/qrp.h \ liboctave/numeric/randgamma.h \ liboctave/numeric/randmtzig.h \ liboctave/numeric/randpoisson.h \ @@ -64,7 +61,6 @@ liboctave/numeric/randpoisson.c NUMERIC_SRC = \ - liboctave/numeric/CmplxQRP.cc \ liboctave/numeric/CollocWt.cc \ liboctave/numeric/DASPK.cc \ liboctave/numeric/DASRT.cc \ @@ -75,11 +71,8 @@ liboctave/numeric/Quad.cc \ liboctave/numeric/aepbalance.cc \ liboctave/numeric/chol.cc \ - liboctave/numeric/dbleQRP.cc \ liboctave/numeric/eigs-base.cc \ - liboctave/numeric/fCmplxQRP.cc \ liboctave/numeric/fEIG.cc \ - liboctave/numeric/floatQRP.cc \ liboctave/numeric/gepbalance.cc \ liboctave/numeric/hess.cc \ liboctave/numeric/lo-mappers.cc \ @@ -91,6 +84,7 @@ liboctave/numeric/oct-rand.cc \ liboctave/numeric/oct-spparms.cc \ liboctave/numeric/qr.cc \ + liboctave/numeric/qrp.cc \ liboctave/numeric/schur.cc \ liboctave/numeric/sparse-chol.cc \ liboctave/numeric/sparse-dmsolve.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/qrp.cc Wed Feb 17 02:54:00 2016 -0500 @@ -0,0 +1,331 @@ +/* + +Copyright (C) 1994-2015 John W. Eatonn +Copyright (C) 2009 VZLU Prague + +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 3 of the License, 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, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include <cassert> + +#include "CMatrix.h" +#include "dMatrix.h" +#include "dRowVector.h" +#include "f77-fcn.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "fRowVector.h" +#include "lo-error.h" +#include "oct-locbuf.h" +#include "qrp.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, + double*, const octave_idx_type&, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (sgeqp3, SGEQP3) (const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, octave_idx_type*, + float*, float*, const octave_idx_type&, + octave_idx_type&); + F77_RET_T + F77_FUNC (zgeqp3, ZGEQP3) (const octave_idx_type&, const octave_idx_type&, + Complex*, const octave_idx_type&, + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, double*, + octave_idx_type&); + F77_RET_T + F77_FUNC (cgeqp3, CGEQP3) (const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + octave_idx_type*, FloatComplex*, FloatComplex*, + const octave_idx_type&, float*, octave_idx_type&); +} + +// Specialization. + +template <> +void +qrp<Matrix>::init (const Matrix& a, type qr_type) +{ + assert (qr_type != qr<Matrix>::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + octave_idx_type min_mn = m < n ? m : n; + OCTAVE_LOCAL_BUFFER (double, tau, min_mn); + + octave_idx_type info = 0; + + Matrix afact = a; + if (m > n && qr_type == qr<Matrix>::std) + afact.resize (m, m); + + MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); + + if (m > 0) + { + // workspace query. + double rlwork; + F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + &rlwork, -1, info)); + + // allocate buffer and do the job. + octave_idx_type lwork = rlwork; + lwork = std::max (lwork, static_cast<octave_idx_type> (1)); + OCTAVE_LOCAL_BUFFER (double, work, lwork); + F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + work, lwork, info)); + } + else + for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + jpvt -= static_cast<octave_idx_type> (1); + p = PermMatrix (jpvt, true); + + + form (n, afact, tau, qr_type); +} + +template <> +qrp<Matrix>::qrp (const Matrix& a, type qr_type) + : qr<Matrix> (), p () +{ + init (a, qr_type); +} + +template <> +RowVector +qrp<Matrix>::Pvec (void) const +{ + Array<double> pa (p.col_perm_vec ()); + RowVector pv (MArray<double> (pa) + 1.0); + return pv; +} + +template <> +void +qrp<FloatMatrix>::init (const FloatMatrix& a, type qr_type) +{ + assert (qr_type != qr<FloatMatrix>::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + octave_idx_type min_mn = m < n ? m : n; + OCTAVE_LOCAL_BUFFER (float, tau, min_mn); + + octave_idx_type info = 0; + + FloatMatrix afact = a; + if (m > n && qr_type == qr<FloatMatrix>::std) + afact.resize (m, m); + + MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); + + if (m > 0) + { + // workspace query. + float rlwork; + F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + &rlwork, -1, info)); + + // allocate buffer and do the job. + octave_idx_type lwork = rlwork; + lwork = std::max (lwork, static_cast<octave_idx_type> (1)); + OCTAVE_LOCAL_BUFFER (float, work, lwork); + F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + work, lwork, info)); + } + else + for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + jpvt -= static_cast<octave_idx_type> (1); + p = PermMatrix (jpvt, true); + + + form (n, afact, tau, qr_type); +} + +template <> +qrp<FloatMatrix>::qrp (const FloatMatrix& a, type qr_type) + : qr<FloatMatrix> (), p () +{ + init (a, qr_type); +} + +template <> +FloatRowVector +qrp<FloatMatrix>::Pvec (void) const +{ + Array<float> pa (p.col_perm_vec ()); + FloatRowVector pv (MArray<float> (pa) + 1.0f); + return pv; +} + +template <> +void +qrp<ComplexMatrix>::init (const ComplexMatrix& a, type qr_type) +{ + assert (qr_type != qr<ComplexMatrix>::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + octave_idx_type min_mn = m < n ? m : n; + OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn); + + octave_idx_type info = 0; + + ComplexMatrix afact = a; + if (m > n && qr_type == qr<ComplexMatrix>::std) + afact.resize (m, m); + + MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); + + if (m > 0) + { + OCTAVE_LOCAL_BUFFER (double, rwork, 2*n); + + // workspace query. + Complex clwork; + F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + &clwork, -1, rwork, info)); + + // allocate buffer and do the job. + octave_idx_type lwork = clwork.real (); + lwork = std::max (lwork, static_cast<octave_idx_type> (1)); + OCTAVE_LOCAL_BUFFER (Complex, work, lwork); + F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + work, lwork, rwork, info)); + } + else + for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + jpvt -= static_cast<octave_idx_type> (1); + p = PermMatrix (jpvt, true); + + + form (n, afact, tau, qr_type); +} + +template <> +qrp<ComplexMatrix>::qrp (const ComplexMatrix& a, type qr_type) + : qr<ComplexMatrix> (), p () +{ + init (a, qr_type); +} + +template <> +RowVector +qrp<ComplexMatrix>::Pvec (void) const +{ + Array<double> pa (p.col_perm_vec ()); + RowVector pv (MArray<double> (pa) + 1.0); + return pv; +} + +template <> +void +qrp<FloatComplexMatrix>::init (const FloatComplexMatrix& a, type qr_type) +{ + assert (qr_type != qr<FloatComplexMatrix>::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + octave_idx_type min_mn = m < n ? m : n; + OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn); + + octave_idx_type info = 0; + + FloatComplexMatrix afact = a; + if (m > n && qr_type == qr<FloatComplexMatrix>::std) + afact.resize (m, m); + + MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); + + if (m > 0) + { + OCTAVE_LOCAL_BUFFER (float, rwork, 2*n); + + // workspace query. + FloatComplex clwork; + F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + &clwork, -1, rwork, info)); + + // allocate buffer and do the job. + octave_idx_type lwork = clwork.real (); + lwork = std::max (lwork, static_cast<octave_idx_type> (1)); + OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork); + F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), + m, jpvt.fortran_vec (), tau, + work, lwork, rwork, info)); + } + else + for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + jpvt -= static_cast<octave_idx_type> (1); + p = PermMatrix (jpvt, true); + + + form (n, afact, tau, qr_type); +} + +template <> +qrp<FloatComplexMatrix>::qrp (const FloatComplexMatrix& a, type qr_type) + : qr<FloatComplexMatrix> (), p () +{ + init (a, qr_type); +} + +template <> +FloatRowVector +qrp<FloatComplexMatrix>::Pvec (void) const +{ + Array<float> pa (p.col_perm_vec ()); + FloatRowVector pv (MArray<float> (pa) + 1.0f); + return pv; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/qrp.h Wed Feb 17 02:54:00 2016 -0500 @@ -0,0 +1,72 @@ +/* + +Copyright (C) 1994-2015 John W. Eaton +Copyright (C) 2009 VZLU Prague + +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 3 of the License, 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, see +<http://www.gnu.org/licenses/>. + +*/ + +#if ! defined (octave_qrp_h) +#define octave_qrp_h 1 + +#include "octave-config.h" + +#include "PermMatrix.h" +#include "qr.h" + +template <typename T> +class +qrp : public qr<T> +{ +public: + + typedef typename T::real_row_vector_type RV_T; + + typedef typename qr<T>::type type; + + qrp (void) : qr<T> (), p () { } + + qrp (const T&, type = qr<T>::std); + + qrp (const qrp& a) : qr<T> (a), p (a.p) { } + + qrp& operator = (const qrp& a) + { + if (this != &a) + { + qr<T>::operator = (a); + p = a.p; + } + + return *this; + } + + ~qrp (void) { } + + void init (const T&, type = qr<T>::std); + + PermMatrix P (void) const { return p; } + + RV_T Pvec (void) const; + +private: + + PermMatrix p; +}; + +#endif
--- a/liboctave/operators/mx-defs.h Wed Feb 17 02:57:21 2016 -0500 +++ b/liboctave/operators/mx-defs.h Wed Feb 17 02:54:00 2016 -0500 @@ -76,10 +76,7 @@ template <typename T> class qr; -class QRP; -class ComplexQRP; -class FloatQRP; -class FloatComplexQRP; +template <typename T> class qrp; // Other data types we use but that don't always need to have full // declarations.