Mercurial > octave
changeset 21271:7e67c7f82fc1
better use of templates for lu factorization classes
* liboctave/numeric/lu.h, liboctave/numeric/lu.cc:
New files generated from base-lu.h, base-lu.cc, CmplxLU.cc, CmplxLU.h,
dbleLU.cc, dbleLU.h, fCmplxLU.cc, fCmplxLU.h, floatLU.cc, and
floatLU.h and converted to templates.
* liboctave/numeric/module.mk: Update.
* lu.cc, mx-defs.h, mx-ext.h, eigs-base.cc: Use new classes.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 16 Feb 2016 12:58:32 -0500 |
| parents | 230e186e292d |
| children | 987c1a79d33f |
| files | libinterp/corefcn/lu.cc liboctave/numeric/CmplxLU.cc liboctave/numeric/CmplxLU.h liboctave/numeric/base-lu.cc liboctave/numeric/base-lu.h liboctave/numeric/dbleLU.cc liboctave/numeric/dbleLU.h liboctave/numeric/eigs-base.cc liboctave/numeric/fCmplxLU.cc liboctave/numeric/fCmplxLU.h liboctave/numeric/floatLU.cc liboctave/numeric/floatLU.h liboctave/numeric/lu.cc liboctave/numeric/lu.h liboctave/numeric/module.mk liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 17 files changed, 1000 insertions(+), 1480 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/lu.cc Tue Feb 16 11:13:55 2016 -0500 +++ b/libinterp/corefcn/lu.cc Tue Feb 16 12:58:32 2016 -0500 @@ -24,10 +24,7 @@ # include <config.h> #endif -#include "CmplxLU.h" -#include "dbleLU.h" -#include "fCmplxLU.h" -#include "floatLU.h" +#include "lu.h" #include "sparse-lu.h" #include "defun.h" @@ -40,7 +37,7 @@ template <typename MT> static octave_value -get_lu_l (const base_lu<MT>& fact) +get_lu_l (const lu<MT>& fact) { MT L = fact.L (); if (L.is_square ()) @@ -51,7 +48,7 @@ template <typename MT> static octave_value -get_lu_u (const base_lu<MT>& fact) +get_lu_u (const lu<MT>& fact) { MT U = fact.U (); if (U.is_square () && fact.regular ()) @@ -344,7 +341,7 @@ { FloatMatrix m = arg.float_matrix_value (); - FloatLU fact (m); + lu<FloatMatrix> fact (m); switch (nargout) { @@ -378,7 +375,7 @@ { Matrix m = arg.matrix_value (); - LU fact (m); + lu<Matrix> fact (m); switch (nargout) { @@ -415,7 +412,7 @@ { FloatComplexMatrix m = arg.float_complex_matrix_value (); - FloatComplexLU fact (m); + lu<FloatComplexMatrix> fact (m); switch (nargout) { @@ -449,7 +446,7 @@ { ComplexMatrix m = arg.complex_matrix_value (); - ComplexLU fact (m); + lu<ComplexMatrix> fact (m); switch (nargout) { @@ -650,7 +647,7 @@ FloatMatrix x = argx.float_matrix_value (); FloatMatrix y = argy.float_matrix_value (); - FloatLU fact (L, U, P); + lu<FloatMatrix> fact (L, U, P); if (pivoted) fact.update_piv (x, y); else @@ -668,7 +665,7 @@ Matrix x = argx.matrix_value (); Matrix y = argy.matrix_value (); - LU fact (L, U, P); + lu<Matrix> fact (L, U, P); if (pivoted) fact.update_piv (x, y); else @@ -691,7 +688,7 @@ FloatComplexMatrix x = argx.float_complex_matrix_value (); FloatComplexMatrix y = argy.float_complex_matrix_value (); - FloatComplexLU fact (L, U, P); + lu<FloatComplexMatrix> fact (L, U, P); if (pivoted) fact.update_piv (x, y); else @@ -709,7 +706,7 @@ ComplexMatrix x = argx.complex_matrix_value (); ComplexMatrix y = argy.complex_matrix_value (); - ComplexLU fact (L, U, P); + lu<ComplexMatrix> fact (L, U, P); if (pivoted) fact.update_piv (x, y); else
--- a/liboctave/numeric/CmplxLU.cc Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,222 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague, a.s. - -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 "CmplxLU.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" -#include "CColVector.h" - -// Instantiate the base LU class for the types we need. - -#include "base-lu.h" -#include "base-lu.cc" - -template class base_lu <ComplexMatrix>; - -// Define the constructor for this particular derivation. - -extern "C" -{ - F77_RET_T - F77_FUNC (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&, - Complex*, const octave_idx_type&, - octave_idx_type*, octave_idx_type&); - -#ifdef HAVE_QRUPDATE_LUU - F77_RET_T - F77_FUNC (zlu1up, ZLU1UP) (const octave_idx_type&, const octave_idx_type&, - Complex *, const octave_idx_type&, - Complex *, const octave_idx_type&, - Complex *, Complex *); - - F77_RET_T - F77_FUNC (zlup1up, ZLUP1UP) (const octave_idx_type&, const octave_idx_type&, - Complex *, const octave_idx_type&, - Complex *, const octave_idx_type&, - octave_idx_type *, const Complex *, - const Complex *, Complex *); -#endif -} - -ComplexLU::ComplexLU (const ComplexMatrix& a) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - ipvt.resize (dim_vector (mn, 1)); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - a_fact = a; - Complex *tmp_data = a_fact.fortran_vec (); - - octave_idx_type info = 0; - - F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - - for (octave_idx_type i = 0; i < mn; i++) - pipvt[i] -= 1; -} - -#ifdef HAVE_QRUPDATE_LUU - -void ComplexLU::update (const ComplexColumnVector& u, - const ComplexColumnVector& v) -{ - if (packed ()) - unpack (); - - ComplexMatrix& l = l_fact; - ComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - ComplexColumnVector utmp = u; - ComplexColumnVector vtmp = v; - F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); -} - -void ComplexLU::update (const ComplexMatrix& u, const ComplexMatrix& v) -{ - if (packed ()) - unpack (); - - ComplexMatrix& l = l_fact; - ComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - ComplexColumnVector utmp = u.column (i); - ComplexColumnVector vtmp = v.column (i); - F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); - } -} - -void ComplexLU::update_piv (const ComplexColumnVector& u, - const ComplexColumnVector& v) -{ - if (packed ()) - unpack (); - - ComplexMatrix& l = l_fact; - ComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - ComplexColumnVector utmp = u; - ComplexColumnVector vtmp = v; - OCTAVE_LOCAL_BUFFER (Complex, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -void ComplexLU::update_piv (const ComplexMatrix& u, const ComplexMatrix& v) -{ - if (packed ()) - unpack (); - - ComplexMatrix& l = l_fact; - ComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - OCTAVE_LOCAL_BUFFER (Complex, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - ComplexColumnVector utmp = u.column (i); - ComplexColumnVector vtmp = v.column (i); - F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - } - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -#else - -void ComplexLU::update (const ComplexColumnVector&, const ComplexColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void ComplexLU::update (const ComplexMatrix&, const ComplexMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void ComplexLU::update_piv (const ComplexColumnVector&, - const ComplexColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void ComplexLU::update_piv (const ComplexMatrix&, const ComplexMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -#endif
--- a/liboctave/numeric/CmplxLU.h Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +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_CmplxLU_h) -#define octave_CmplxLU_h 1 - -#include "octave-config.h" - -#include "base-lu.h" -#include "dMatrix.h" -#include "CMatrix.h" - -class -OCTAVE_API -ComplexLU : public base_lu <ComplexMatrix> -{ -public: - - ComplexLU (void) - : base_lu <ComplexMatrix> () { } - - ComplexLU (const ComplexMatrix& a); - - ComplexLU (const ComplexLU& a) - : base_lu <ComplexMatrix> (a) { } - - ComplexLU (const ComplexMatrix& l, const ComplexMatrix& u, - const PermMatrix& p) - : base_lu <ComplexMatrix> (l, u, p) { } - - ComplexLU& operator = (const ComplexLU& a) - { - if (this != &a) - base_lu <ComplexMatrix> :: operator = (a); - - return *this; - } - - ~ComplexLU (void) { } - - void update (const ComplexColumnVector& u, const ComplexColumnVector& v); - - void update (const ComplexMatrix& u, const ComplexMatrix& v); - - void update_piv (const ComplexColumnVector& u, const ComplexColumnVector& v); - - void update_piv (const ComplexMatrix& u, const ComplexMatrix& v); -}; - -#endif
--- a/liboctave/numeric/base-lu.cc Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,192 +0,0 @@ -/* - -Copyright (C) 1996-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 "base-lu.h" - -template <typename lu_type> -base_lu<lu_type>::base_lu (const lu_type& l, const lu_type& u, - const PermMatrix& p) - : a_fact (u), l_fact (l), ipvt (p.transpose ().col_perm_vec ()) -{ - if (l.columns () != u.rows ()) - (*current_liboctave_error_handler) ("lu: dimension mismatch"); -} - -template <typename lu_type> -bool -base_lu <lu_type> :: packed (void) const -{ - return l_fact.dims () == dim_vector (); -} - -template <typename lu_type> -void -base_lu <lu_type> :: unpack (void) -{ - if (packed ()) - { - l_fact = L (); - a_fact = U (); // FIXME: sub-optimal - ipvt = getp (); - } -} - -template <typename lu_type> -lu_type -base_lu <lu_type> :: L (void) const -{ - if (packed ()) - { - octave_idx_type a_nr = a_fact.rows (); - octave_idx_type a_nc = a_fact.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - lu_type l (a_nr, mn, lu_elt_type (0.0)); - - for (octave_idx_type i = 0; i < a_nr; i++) - { - if (i < a_nc) - l.xelem (i, i) = 1.0; - - for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++) - l.xelem (i, j) = a_fact.xelem (i, j); - } - - return l; - } - else - return l_fact; -} - -template <typename lu_type> -lu_type -base_lu <lu_type> :: U (void) const -{ - if (packed ()) - { - octave_idx_type a_nr = a_fact.rows (); - octave_idx_type a_nc = a_fact.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - lu_type u (mn, a_nc, lu_elt_type (0.0)); - - for (octave_idx_type i = 0; i < mn; i++) - { - for (octave_idx_type j = i; j < a_nc; j++) - u.xelem (i, j) = a_fact.xelem (i, j); - } - - return u; - } - else - return a_fact; -} - -template <typename lu_type> -lu_type -base_lu <lu_type> :: Y (void) const -{ - if (! packed ()) - (*current_liboctave_error_handler) - ("lu: Y () not implemented for unpacked form"); - - return a_fact; -} - -template <typename lu_type> -Array<octave_idx_type> -base_lu <lu_type> :: getp (void) const -{ - if (packed ()) - { - octave_idx_type a_nr = a_fact.rows (); - - Array<octave_idx_type> pvt (dim_vector (a_nr, 1)); - - for (octave_idx_type i = 0; i < a_nr; i++) - pvt.xelem (i) = i; - - for (octave_idx_type i = 0; i < ipvt.numel (); i++) - { - octave_idx_type k = ipvt.xelem (i); - - if (k != i) - { - octave_idx_type tmp = pvt.xelem (k); - pvt.xelem (k) = pvt.xelem (i); - pvt.xelem (i) = tmp; - } - } - - return pvt; - } - else - return ipvt; -} - -template <typename lu_type> -PermMatrix -base_lu <lu_type> :: P (void) const -{ - return PermMatrix (getp (), false); -} - -template <typename lu_type> -ColumnVector -base_lu <lu_type> :: P_vec (void) const -{ - octave_idx_type a_nr = a_fact.rows (); - - ColumnVector p (a_nr); - - Array<octave_idx_type> pvt = getp (); - - for (octave_idx_type i = 0; i < a_nr; i++) - p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1); - - return p; -} - -template <typename lu_type> -bool -base_lu<lu_type>::regular (void) const -{ - bool retval = true; - - octave_idx_type k = std::min (a_fact.rows (), a_fact.columns ()); - - for (octave_idx_type i = 0; i < k; i++) - { - if (a_fact(i, i) == lu_elt_type ()) - { - retval = false; - break; - } - } - - return retval; -}
--- a/liboctave/numeric/base-lu.h Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -/* - -Copyright (C) 1996-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_base_lu_h) -#define octave_base_lu_h 1 - -#include "octave-config.h" - -#include "MArray.h" -#include "dColVector.h" -#include "PermMatrix.h" - -template <typename lu_type> -class -base_lu -{ -public: - - typedef typename lu_type::element_type lu_elt_type; - - base_lu (void) - : a_fact (), l_fact (), ipvt () { } - - base_lu (const base_lu& a) - : a_fact (a.a_fact), l_fact (a.l_fact), ipvt (a.ipvt) { } - - base_lu (const lu_type& l, const lu_type& u, - const PermMatrix& p); - - base_lu& operator = (const base_lu& a) - { - if (this != &a) - { - a_fact = a.a_fact; - l_fact = a.l_fact; - ipvt = a.ipvt; - } - return *this; - } - - virtual ~base_lu (void) { } - - bool packed (void) const; - - void unpack (void); - - lu_type L (void) const; - - lu_type U (void) const; - - lu_type Y (void) const; - - PermMatrix P (void) const; - - ColumnVector P_vec (void) const; - - bool regular (void) const; - -protected: - - Array<octave_idx_type> getp (void) const; - - lu_type a_fact; - lu_type l_fact; - - Array<octave_idx_type> ipvt; -}; - -#endif
--- a/liboctave/numeric/dbleLU.cc Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,219 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague, a.s. - -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 "dbleLU.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" -#include "dColVector.h" - -// Instantiate the base LU class for the types we need. - -#include "base-lu.h" -#include "base-lu.cc" - -template class base_lu <Matrix>; - -// Define the constructor for this particular derivation. - -extern "C" -{ - F77_RET_T - F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&, - double*, const octave_idx_type&, - octave_idx_type*, octave_idx_type&); - -#ifdef HAVE_QRUPDATE_LUU - F77_RET_T - F77_FUNC (dlu1up, DLU1UP) (const octave_idx_type&, const octave_idx_type&, - double *, const octave_idx_type&, - double *, const octave_idx_type&, - double *, double *); - - F77_RET_T - F77_FUNC (dlup1up, DLUP1UP) (const octave_idx_type&, const octave_idx_type&, - double *, const octave_idx_type&, - double *, const octave_idx_type&, - octave_idx_type *, const double *, - const double *, double *); -#endif -} - -LU::LU (const Matrix& a) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - ipvt.resize (dim_vector (mn, 1)); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - a_fact = a; - double *tmp_data = a_fact.fortran_vec (); - - octave_idx_type info = 0; - - F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - - for (octave_idx_type i = 0; i < mn; i++) - pipvt[i] -= 1; -} - -#ifdef HAVE_QRUPDATE_LUU - -void LU::update (const ColumnVector& u, const ColumnVector& v) -{ - if (packed ()) - unpack (); - - Matrix& l = l_fact; - Matrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - ColumnVector utmp = u; - ColumnVector vtmp = v; - F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); -} - -void LU::update (const Matrix& u, const Matrix& v) -{ - if (packed ()) - unpack (); - - Matrix& l = l_fact; - Matrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - ColumnVector utmp = u.column (i); - ColumnVector vtmp = v.column (i); - F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); - } -} - -void LU::update_piv (const ColumnVector& u, const ColumnVector& v) -{ - if (packed ()) - unpack (); - - Matrix& l = l_fact; - Matrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - ColumnVector utmp = u; - ColumnVector vtmp = v; - OCTAVE_LOCAL_BUFFER (double, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -void LU::update_piv (const Matrix& u, const Matrix& v) -{ - if (packed ()) - unpack (); - - Matrix& l = l_fact; - Matrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - OCTAVE_LOCAL_BUFFER (double, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - ColumnVector utmp = u.column (i); - ColumnVector vtmp = v.column (i); - F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - } - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -#else - -void LU::update (const ColumnVector&, const ColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void LU::update (const Matrix&, const Matrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void LU::update_piv (const ColumnVector&, const ColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void LU::update_piv (const Matrix&, const Matrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -#endif
--- a/liboctave/numeric/dbleLU.h Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,65 +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_dbleLU_h) -#define octave_dbleLU_h 1 - -#include "octave-config.h" - -#include "base-lu.h" -#include "dMatrix.h" - -class -OCTAVE_API -LU : public base_lu <Matrix> -{ -public: - - LU (void) : base_lu <Matrix> () { } - - LU (const Matrix& a); - - LU (const LU& a) : base_lu <Matrix> (a) { } - - LU (const Matrix& l, const Matrix& u, const PermMatrix& p) - : base_lu <Matrix> (l, u, p) { } - - LU& operator = (const LU& a) - { - if (this != &a) - base_lu <Matrix> :: operator = (a); - - return *this; - } - - ~LU (void) { } - - void update (const ColumnVector& u, const ColumnVector& v); - - void update (const Matrix& u, const Matrix& v); - - void update_piv (const ColumnVector& u, const ColumnVector& v); - - void update_piv (const Matrix& u, const Matrix& v); -}; - -#endif
--- a/liboctave/numeric/eigs-base.cc Tue Feb 16 11:13:55 2016 -0500 +++ b/liboctave/numeric/eigs-base.cc Tue Feb 16 12:58:32 2016 -0500 @@ -30,11 +30,10 @@ #include <iostream> #include "CSparse.h" -#include "CmplxLU.h" +#include "lu.h" #include "MatrixType.h" #include "chol.h" #include "dSparse.h" -#include "dbleLU.h" #include "eigs-base.h" #include "f77-fcn.h" #include "mx-ops.h" @@ -510,7 +509,7 @@ p[i*(n+1)] -= sigma; } - LU fact (AminusSigmaB); + lu<Matrix> fact (AminusSigmaB); L = fact.P ().transpose () * fact.L (); U = fact.U (); @@ -674,7 +673,7 @@ p[i*(n+1)] -= sigma; } - ComplexLU fact (AminusSigmaB); + lu<ComplexMatrix> fact (AminusSigmaB); L = fact.P ().transpose () * fact.L (); U = fact.U ();
--- a/liboctave/numeric/fCmplxLU.cc Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,229 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague, a.s. - -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 "fCmplxLU.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" -#include "fCColVector.h" - -// Instantiate the base LU class for the types we need. - -#include "base-lu.h" -#include "base-lu.cc" - -template class base_lu <FloatComplexMatrix>; - -// Define the constructor for this particular derivation. - -extern "C" -{ - F77_RET_T - F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - octave_idx_type*, octave_idx_type&); - -#ifdef HAVE_QRUPDATE_LUU - F77_RET_T - F77_FUNC (clu1up, CLU1UP) (const octave_idx_type&, const octave_idx_type&, - FloatComplex *, const octave_idx_type&, - FloatComplex *, const octave_idx_type&, - FloatComplex *, FloatComplex *); - - F77_RET_T - F77_FUNC (clup1up, CLUP1UP) (const octave_idx_type&, const octave_idx_type&, - FloatComplex *, const octave_idx_type&, - FloatComplex *, const octave_idx_type&, - octave_idx_type *, const FloatComplex *, - const FloatComplex *, FloatComplex *); -#endif -} - -FloatComplexLU::FloatComplexLU (const FloatComplexMatrix& a) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - ipvt.resize (dim_vector (mn, 1)); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - a_fact = a; - FloatComplex *tmp_data = a_fact.fortran_vec (); - - octave_idx_type info = 0; - - F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - - for (octave_idx_type i = 0; i < mn; i++) - pipvt[i] -= 1; -} - -#ifdef HAVE_QRUPDATE_LUU - -void FloatComplexLU::update (const FloatComplexColumnVector& u, - const FloatComplexColumnVector& v) -{ - if (packed ()) - unpack (); - - FloatComplexMatrix& l = l_fact; - FloatComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () == m && v.numel () == n) - { - FloatComplexColumnVector utmp = u; - FloatComplexColumnVector vtmp = v; - F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); - } - else - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); -} - -void FloatComplexLU::update (const FloatComplexMatrix& u, - const FloatComplexMatrix& v) -{ - if (packed ()) - unpack (); - - FloatComplexMatrix& l = l_fact; - FloatComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - FloatComplexColumnVector utmp = u.column (i); - FloatComplexColumnVector vtmp = v.column (i); - F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); - } -} - -void FloatComplexLU::update_piv (const FloatComplexColumnVector& u, - const FloatComplexColumnVector& v) -{ - if (packed ()) - unpack (); - - FloatComplexMatrix& l = l_fact; - FloatComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - FloatComplexColumnVector utmp = u; - FloatComplexColumnVector vtmp = v; - OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -void FloatComplexLU::update_piv (const FloatComplexMatrix& u, - const FloatComplexMatrix& v) -{ - if (packed ()) - unpack (); - - FloatComplexMatrix& l = l_fact; - FloatComplexMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - FloatComplexColumnVector utmp = u.column (i); - FloatComplexColumnVector vtmp = v.column (i); - F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - } - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -#else - -void FloatComplexLU::update (const FloatComplexColumnVector&, - const FloatComplexColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatComplexLU::update (const FloatComplexMatrix&, - const FloatComplexMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatComplexLU::update_piv (const FloatComplexColumnVector&, - const FloatComplexColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatComplexLU::update_piv (const FloatComplexMatrix&, - const FloatComplexMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -#endif
--- a/liboctave/numeric/fCmplxLU.h Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,71 +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_fCmplxLU_h) -#define octave_fCmplxLU_h 1 - -#include "octave-config.h" - -#include "base-lu.h" -#include "dMatrix.h" -#include "fCMatrix.h" - -class -OCTAVE_API -FloatComplexLU : public base_lu <FloatComplexMatrix> -{ -public: - - FloatComplexLU (void) - : base_lu <FloatComplexMatrix> () { } - - FloatComplexLU (const FloatComplexMatrix& a); - - FloatComplexLU (const FloatComplexLU& a) - : base_lu <FloatComplexMatrix> (a) { } - - FloatComplexLU (const FloatComplexMatrix& l, const FloatComplexMatrix& u, - const PermMatrix& p) - : base_lu <FloatComplexMatrix> (l, u, p) { } - - FloatComplexLU& operator = (const FloatComplexLU& a) - { - if (this != &a) - base_lu <FloatComplexMatrix> :: operator = (a); - - return *this; - } - - ~FloatComplexLU (void) { } - - void update (const FloatComplexColumnVector& u, - const FloatComplexColumnVector& v); - - void update (const FloatComplexMatrix& u, const FloatComplexMatrix& v); - - void update_piv (const FloatComplexColumnVector& u, - const FloatComplexColumnVector& v); - - void update_piv (const FloatComplexMatrix& u, const FloatComplexMatrix& v); -}; - -#endif
--- a/liboctave/numeric/floatLU.cc Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,221 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2009 VZLU Prague, a.s. - -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 "floatLU.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" -#include "fColVector.h" - -// Instantiate the base LU class for the types we need. - -#include "base-lu.h" -#include "base-lu.cc" - -template class base_lu <FloatMatrix>; - -// Define the constructor for this particular derivation. - -extern "C" -{ - F77_RET_T - F77_FUNC (sgetrf, SGETRF) (const octave_idx_type&, const octave_idx_type&, - float*, const octave_idx_type&, octave_idx_type*, - octave_idx_type&); - -#ifdef HAVE_QRUPDATE_LUU - F77_RET_T - F77_FUNC (slu1up, SLU1UP) (const octave_idx_type&, const octave_idx_type&, - float *, const octave_idx_type&, - float *, const octave_idx_type&, - float *, float *); - - F77_RET_T - F77_FUNC (slup1up, SLUP1UP) (const octave_idx_type&, const octave_idx_type&, - float *, const octave_idx_type&, - float *, const octave_idx_type&, - octave_idx_type *, const float *, - const float *, float *); -#endif -} - -FloatLU::FloatLU (const FloatMatrix& a) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); - - ipvt.resize (dim_vector (mn, 1)); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - a_fact = a; - float *tmp_data = a_fact.fortran_vec (); - - octave_idx_type info = 0; - - F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - - for (octave_idx_type i = 0; i < mn; i++) - pipvt[i] -= 1; -} - -#ifdef HAVE_QRUPDATE_LUU - -void FloatLU::update (const FloatColumnVector& u, const FloatColumnVector& v) -{ - if (packed ()) - unpack (); - - FloatMatrix& l = l_fact; - FloatMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - FloatColumnVector utmp = u; - FloatColumnVector vtmp = v; - F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); -} - -void FloatLU::update (const FloatMatrix& u, const FloatMatrix& v) -{ - if (packed ()) - unpack (); - - FloatMatrix& l = l_fact; - FloatMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - FloatColumnVector utmp = u.column (i); - FloatColumnVector vtmp = v.column (i); - F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - utmp.fortran_vec (), vtmp.fortran_vec ())); - } -} - -void FloatLU::update_piv (const FloatColumnVector& u, - const FloatColumnVector& v) -{ - if (packed ()) - unpack (); - - FloatMatrix& l = l_fact; - FloatMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.numel () != m || v.numel () != n) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - FloatColumnVector utmp = u; - FloatColumnVector vtmp = v; - OCTAVE_LOCAL_BUFFER (float, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -void FloatLU::update_piv (const FloatMatrix& u, const FloatMatrix& v) -{ - if (packed ()) - unpack (); - - FloatMatrix& l = l_fact; - FloatMatrix& r = a_fact; - - octave_idx_type m = l.rows (); - octave_idx_type n = r.columns (); - octave_idx_type k = l.columns (); - - if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) - (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); - - OCTAVE_LOCAL_BUFFER (float, w, m); - for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment - for (volatile octave_idx_type i = 0; i < u.cols (); i++) - { - FloatColumnVector utmp = u.column (i); - FloatColumnVector vtmp = v.column (i); - F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), - m, r.fortran_vec (), k, - ipvt.fortran_vec (), - utmp.data (), vtmp.data (), w)); - } - for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement -} - -#else - -void FloatLU::update (const FloatColumnVector&, const FloatColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatLU::update (const FloatMatrix&, const FloatMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatLU::update_piv (const FloatColumnVector&, const FloatColumnVector&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -void FloatLU::update_piv (const FloatMatrix&, const FloatMatrix&) -{ - (*current_liboctave_error_handler) - ("luupdate: support for qrupdate with LU updates " - "was unavailable or disabled when liboctave was built"); -} - -#endif
--- a/liboctave/numeric/floatLU.h Tue Feb 16 11:13:55 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +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_floatLU_h) -#define octave_floatLU_h 1 - -#include "octave-config.h" - -#include "base-lu.h" -#include "dMatrix.h" -#include "fMatrix.h" - -class -OCTAVE_API -FloatLU : public base_lu <FloatMatrix> -{ -public: - - FloatLU (void) : base_lu <FloatMatrix> () { } - - FloatLU (const FloatMatrix& a); - - FloatLU (const FloatLU& a) : base_lu <FloatMatrix> (a) { } - - FloatLU (const FloatMatrix& l, const FloatMatrix& u, - const PermMatrix& p) - : base_lu <FloatMatrix> (l, u, p) { } - - FloatLU& operator = (const FloatLU& a) - { - if (this != &a) - base_lu <FloatMatrix> :: operator = (a); - - return *this; - } - - ~FloatLU (void) { } - - void update (const FloatColumnVector& u, const FloatColumnVector& v); - - void update (const FloatMatrix& u, const FloatMatrix& v); - - void update_piv (const FloatColumnVector& u, const FloatColumnVector& v); - - void update_piv (const FloatMatrix& u, const FloatMatrix& v); -}; - -#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/lu.cc Tue Feb 16 12:58:32 2016 -0500 @@ -0,0 +1,884 @@ +/* + +Copyright (C) 1996-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 "CColVector.h" +#include "CMatrix.h" +#include "dColVector.h" +#include "dMatrix.h" +#include "f77-fcn.h" +#include "fCColVector.h" +#include "fCMatrix.h" +#include "fColVector.h" +#include "fMatrix.h" +#include "lo-error.h" +#include "lu.h" +#include "oct-locbuf.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&, + double*, const octave_idx_type&, + octave_idx_type*, octave_idx_type&); + +#ifdef HAVE_QRUPDATE_LUU + F77_RET_T + F77_FUNC (dlu1up, DLU1UP) (const octave_idx_type&, const octave_idx_type&, + double *, const octave_idx_type&, + double *, const octave_idx_type&, + double *, double *); + + F77_RET_T + F77_FUNC (dlup1up, DLUP1UP) (const octave_idx_type&, const octave_idx_type&, + double *, const octave_idx_type&, + double *, const octave_idx_type&, + octave_idx_type *, const double *, + const double *, double *); +#endif + + F77_RET_T + F77_FUNC (sgetrf, SGETRF) (const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, octave_idx_type*, + octave_idx_type&); + +#ifdef HAVE_QRUPDATE_LUU + F77_RET_T + F77_FUNC (slu1up, SLU1UP) (const octave_idx_type&, const octave_idx_type&, + float *, const octave_idx_type&, + float *, const octave_idx_type&, + float *, float *); + + F77_RET_T + F77_FUNC (slup1up, SLUP1UP) (const octave_idx_type&, const octave_idx_type&, + float *, const octave_idx_type&, + float *, const octave_idx_type&, + octave_idx_type *, const float *, + const float *, float *); +#endif + + F77_RET_T + F77_FUNC (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&, + Complex*, const octave_idx_type&, + octave_idx_type*, octave_idx_type&); + +#ifdef HAVE_QRUPDATE_LUU + F77_RET_T + F77_FUNC (zlu1up, ZLU1UP) (const octave_idx_type&, const octave_idx_type&, + Complex *, const octave_idx_type&, + Complex *, const octave_idx_type&, + Complex *, Complex *); + + F77_RET_T + F77_FUNC (zlup1up, ZLUP1UP) (const octave_idx_type&, const octave_idx_type&, + Complex *, const octave_idx_type&, + Complex *, const octave_idx_type&, + octave_idx_type *, const Complex *, + const Complex *, Complex *); +#endif + + F77_RET_T + F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + octave_idx_type*, octave_idx_type&); + +#ifdef HAVE_QRUPDATE_LUU + F77_RET_T + F77_FUNC (clu1up, CLU1UP) (const octave_idx_type&, const octave_idx_type&, + FloatComplex *, const octave_idx_type&, + FloatComplex *, const octave_idx_type&, + FloatComplex *, FloatComplex *); + + F77_RET_T + F77_FUNC (clup1up, CLUP1UP) (const octave_idx_type&, const octave_idx_type&, + FloatComplex *, const octave_idx_type&, + FloatComplex *, const octave_idx_type&, + octave_idx_type *, const FloatComplex *, + const FloatComplex *, FloatComplex *); +#endif +} + +template <typename T> +lu<T>::lu (const T& l, const T& u, + const PermMatrix& p) + : a_fact (u), l_fact (l), ipvt (p.transpose ().col_perm_vec ()) +{ + if (l.columns () != u.rows ()) + (*current_liboctave_error_handler) ("lu: dimension mismatch"); +} + +template <typename T> +bool +lu<T>::packed (void) const +{ + return l_fact.dims () == dim_vector (); +} + +template <typename T> +void +lu<T>::unpack (void) +{ + if (packed ()) + { + l_fact = L (); + a_fact = U (); // FIXME: sub-optimal + ipvt = getp (); + } +} + +template <typename T> +T +lu<T>::L (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + octave_idx_type a_nc = a_fact.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + T l (a_nr, mn, ELT_T (0.0)); + + for (octave_idx_type i = 0; i < a_nr; i++) + { + if (i < a_nc) + l.xelem (i, i) = 1.0; + + for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++) + l.xelem (i, j) = a_fact.xelem (i, j); + } + + return l; + } + else + return l_fact; +} + +template <typename T> +T +lu<T>::U (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + octave_idx_type a_nc = a_fact.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + T u (mn, a_nc, ELT_T (0.0)); + + for (octave_idx_type i = 0; i < mn; i++) + { + for (octave_idx_type j = i; j < a_nc; j++) + u.xelem (i, j) = a_fact.xelem (i, j); + } + + return u; + } + else + return a_fact; +} + +template <typename T> +T +lu<T>::Y (void) const +{ + if (! packed ()) + (*current_liboctave_error_handler) + ("lu: Y () not implemented for unpacked form"); + + return a_fact; +} + +template <typename T> +Array<octave_idx_type> +lu<T>::getp (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + + Array<octave_idx_type> pvt (dim_vector (a_nr, 1)); + + for (octave_idx_type i = 0; i < a_nr; i++) + pvt.xelem (i) = i; + + for (octave_idx_type i = 0; i < ipvt.numel (); i++) + { + octave_idx_type k = ipvt.xelem (i); + + if (k != i) + { + octave_idx_type tmp = pvt.xelem (k); + pvt.xelem (k) = pvt.xelem (i); + pvt.xelem (i) = tmp; + } + } + + return pvt; + } + else + return ipvt; +} + +template <typename T> +PermMatrix +lu<T>::P (void) const +{ + return PermMatrix (getp (), false); +} + +template <typename T> +ColumnVector +lu<T>::P_vec (void) const +{ + octave_idx_type a_nr = a_fact.rows (); + + ColumnVector p (a_nr); + + Array<octave_idx_type> pvt = getp (); + + for (octave_idx_type i = 0; i < a_nr; i++) + p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1); + + return p; +} + +template <typename T> +bool +lu<T>::regular (void) const +{ + bool retval = true; + + octave_idx_type k = std::min (a_fact.rows (), a_fact.columns ()); + + for (octave_idx_type i = 0; i < k; i++) + { + if (a_fact(i, i) == ELT_T ()) + { + retval = false; + break; + } + } + + return retval; +} + +#if ! defined (HAVE_QRUPDATE_LUU) + +template <typename T> +void +lu<T>::update (const VT&, const VT&) +{ + (*current_liboctave_error_handler) + ("luupdate: support for qrupdate with LU updates " + "was unavailable or disabled when liboctave was built"); +} + +template <typename T> +void +lu<T>::update (const T&, const T&) +{ + (*current_liboctave_error_handler) + ("luupdate: support for qrupdate with LU updates " + "was unavailable or disabled when liboctave was built"); +} + +template <typename T> +void +lu<T>::update_piv (const VT&, const VT&) +{ + (*current_liboctave_error_handler) + ("luupdate: support for qrupdate with LU updates " + "was unavailable or disabled when liboctave was built"); +} + +template <typename T> +void +lu<T>::update_piv (const T&, const T&) +{ + (*current_liboctave_error_handler) + ("luupdate: support for qrupdate with LU updates " + "was unavailable or disabled when liboctave was built"); +} + +#endif + +// Specializations. + +template <> +lu<Matrix>::lu (const Matrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + ipvt.resize (dim_vector (mn, 1)); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + a_fact = a; + double *tmp_data = a_fact.fortran_vec (); + + octave_idx_type info = 0; + + F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); + + for (octave_idx_type i = 0; i < mn; i++) + pipvt[i] -= 1; +} + +#ifdef HAVE_QRUPDATE_LUU + +template <> +void +lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v) +{ + if (packed ()) + unpack (); + + Matrix& l = l_fact; + Matrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + ColumnVector utmp = u; + ColumnVector vtmp = v; + F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); +} + +template <> +void +lu<Matrix>::update (const Matrix& u, const Matrix& v) +{ + if (packed ()) + unpack (); + + Matrix& l = l_fact; + Matrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + ColumnVector utmp = u.column (i); + ColumnVector vtmp = v.column (i); + F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); + } +} + +template <> +void +lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v) +{ + if (packed ()) + unpack (); + + Matrix& l = l_fact; + Matrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + ColumnVector utmp = u; + ColumnVector vtmp = v; + OCTAVE_LOCAL_BUFFER (double, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +template <> +void +lu<Matrix>::update_piv (const Matrix& u, const Matrix& v) +{ + if (packed ()) + unpack (); + + Matrix& l = l_fact; + Matrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + OCTAVE_LOCAL_BUFFER (double, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + ColumnVector utmp = u.column (i); + ColumnVector vtmp = v.column (i); + F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + } + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +#endif + +template <> +lu<FloatMatrix>::lu (const FloatMatrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + ipvt.resize (dim_vector (mn, 1)); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + a_fact = a; + float *tmp_data = a_fact.fortran_vec (); + + octave_idx_type info = 0; + + F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); + + for (octave_idx_type i = 0; i < mn; i++) + pipvt[i] -= 1; +} + +#ifdef HAVE_QRUPDATE_LUU + +template <> +void +lu<FloatMatrix>::update (const FloatColumnVector& u, const FloatColumnVector& v) +{ + if (packed ()) + unpack (); + + FloatMatrix& l = l_fact; + FloatMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + FloatColumnVector utmp = u; + FloatColumnVector vtmp = v; + F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); +} + +template <> +void +lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v) +{ + if (packed ()) + unpack (); + + FloatMatrix& l = l_fact; + FloatMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + FloatColumnVector utmp = u.column (i); + FloatColumnVector vtmp = v.column (i); + F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); + } +} + +template <> +void +lu<FloatMatrix>::update_piv (const FloatColumnVector& u, + const FloatColumnVector& v) +{ + if (packed ()) + unpack (); + + FloatMatrix& l = l_fact; + FloatMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + FloatColumnVector utmp = u; + FloatColumnVector vtmp = v; + OCTAVE_LOCAL_BUFFER (float, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +template <> +void +lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v) +{ + if (packed ()) + unpack (); + + FloatMatrix& l = l_fact; + FloatMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + OCTAVE_LOCAL_BUFFER (float, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + FloatColumnVector utmp = u.column (i); + FloatColumnVector vtmp = v.column (i); + F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + } + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +#endif + +template <> +lu<ComplexMatrix>::lu (const ComplexMatrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + ipvt.resize (dim_vector (mn, 1)); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + a_fact = a; + Complex *tmp_data = a_fact.fortran_vec (); + + octave_idx_type info = 0; + + F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); + + for (octave_idx_type i = 0; i < mn; i++) + pipvt[i] -= 1; +} + +#ifdef HAVE_QRUPDATE_LUU + +template <> +void +lu<ComplexMatrix>::update (const ComplexColumnVector& u, + const ComplexColumnVector& v) +{ + if (packed ()) + unpack (); + + ComplexMatrix& l = l_fact; + ComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + ComplexColumnVector utmp = u; + ComplexColumnVector vtmp = v; + F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); +} + +template <> +void +lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v) +{ + if (packed ()) + unpack (); + + ComplexMatrix& l = l_fact; + ComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + ComplexColumnVector utmp = u.column (i); + ComplexColumnVector vtmp = v.column (i); + F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); + } +} + +template <> +void +lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u, + const ComplexColumnVector& v) +{ + if (packed ()) + unpack (); + + ComplexMatrix& l = l_fact; + ComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + ComplexColumnVector utmp = u; + ComplexColumnVector vtmp = v; + OCTAVE_LOCAL_BUFFER (Complex, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +template <> +void +lu<ComplexMatrix>::update_piv (const ComplexMatrix& u, const ComplexMatrix& v) +{ + if (packed ()) + unpack (); + + ComplexMatrix& l = l_fact; + ComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + OCTAVE_LOCAL_BUFFER (Complex, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + ComplexColumnVector utmp = u.column (i); + ComplexColumnVector vtmp = v.column (i); + F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + } + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +#endif + +template <> +lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + ipvt.resize (dim_vector (mn, 1)); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + a_fact = a; + FloatComplex *tmp_data = a_fact.fortran_vec (); + + octave_idx_type info = 0; + + F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); + + for (octave_idx_type i = 0; i < mn; i++) + pipvt[i] -= 1; +} + +#ifdef HAVE_QRUPDATE_LUU + +template <> +void +lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u, + const FloatComplexColumnVector& v) +{ + if (packed ()) + unpack (); + + FloatComplexMatrix& l = l_fact; + FloatComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () == m && v.numel () == n) + { + FloatComplexColumnVector utmp = u; + FloatComplexColumnVector vtmp = v; + F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); + } + else + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); +} + +template <> +void +lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u, + const FloatComplexMatrix& v) +{ + if (packed ()) + unpack (); + + FloatComplexMatrix& l = l_fact; + FloatComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + FloatComplexColumnVector utmp = u.column (i); + FloatComplexColumnVector vtmp = v.column (i); + F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec ())); + } +} + +template <> +void +lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u, + const FloatComplexColumnVector& v) +{ + if (packed ()) + unpack (); + + FloatComplexMatrix& l = l_fact; + FloatComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.numel () != m || v.numel () != n) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + FloatComplexColumnVector utmp = u; + FloatComplexColumnVector vtmp = v; + OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +template <> +void +lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u, + const FloatComplexMatrix& v) +{ + if (packed ()) + unpack (); + + FloatComplexMatrix& l = l_fact; + FloatComplexMatrix& r = a_fact; + + octave_idx_type m = l.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = l.columns (); + + if (u.rows () != m || v.rows () != n || u.cols () != v.cols ()) + (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); + + OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); + for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment + for (volatile octave_idx_type i = 0; i < u.cols (); i++) + { + FloatComplexColumnVector utmp = u.column (i); + FloatComplexColumnVector vtmp = v.column (i); + F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), + m, r.fortran_vec (), k, + ipvt.fortran_vec (), + utmp.data (), vtmp.data (), w)); + } + for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement +} + +#endif + +// Instantiations we need. + +template class lu<Matrix>; + +template class lu<FloatMatrix>; + +template class lu<ComplexMatrix>; + +template class lu<FloatComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/lu.h Tue Feb 16 12:58:32 2016 -0500 @@ -0,0 +1,98 @@ +/* + +Copyright (C) 1996-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_lu_h) +#define octave_lu_h 1 + +#include "octave-config.h" + +#include "PermMatrix.h" + +template <typename T> +class +lu +{ +public: + + typedef typename T::column_vector_type VT; + typedef typename T::element_type ELT_T; + + lu (void) + : a_fact (), l_fact (), ipvt () { } + + lu (const T& a); + + lu (const lu& a) + : a_fact (a.a_fact), l_fact (a.l_fact), ipvt (a.ipvt) { } + + lu (const T& l, const T& u, const PermMatrix& p); + + lu& operator = (const lu& a) + { + if (this != &a) + { + a_fact = a.a_fact; + l_fact = a.l_fact; + ipvt = a.ipvt; + } + + return *this; + } + + virtual ~lu (void) { } + + bool packed (void) const; + + void unpack (void); + + T L (void) const; + + T U (void) const; + + T Y (void) const; + + PermMatrix P (void) const; + + ColumnVector P_vec (void) const; + + bool regular (void) const; + + void update (const VT& u, const VT& v); + + void update (const T& u, const T& v); + + void update_piv (const VT& u, const VT& v); + + void update_piv (const T& u, const T& v); + +protected: + + Array<octave_idx_type> getp (void) const; + + T a_fact; + T l_fact; + + Array<octave_idx_type> ipvt; +}; + +#endif
--- a/liboctave/numeric/module.mk Tue Feb 16 11:13:55 2016 -0500 +++ b/liboctave/numeric/module.mk Tue Feb 16 12:58:32 2016 -0500 @@ -8,7 +8,6 @@ LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in) NUMERIC_INC = \ - liboctave/numeric/CmplxLU.h \ liboctave/numeric/CmplxQR.h \ liboctave/numeric/CmplxQRP.h \ liboctave/numeric/CmplxSVD.h \ @@ -31,23 +30,19 @@ liboctave/numeric/aepbalance.h \ liboctave/numeric/base-dae.h \ liboctave/numeric/base-de.h \ - liboctave/numeric/base-lu.h \ liboctave/numeric/base-min.h \ liboctave/numeric/base-qr.h \ liboctave/numeric/bsxfun-decl.h \ liboctave/numeric/bsxfun.h \ liboctave/numeric/chol.h \ - liboctave/numeric/dbleLU.h \ liboctave/numeric/dbleQR.h \ liboctave/numeric/dbleQRP.h \ liboctave/numeric/dbleSVD.h \ liboctave/numeric/eigs-base.h \ - liboctave/numeric/fCmplxLU.h \ liboctave/numeric/fCmplxQR.h \ liboctave/numeric/fCmplxQRP.h \ liboctave/numeric/fCmplxSVD.h \ liboctave/numeric/fEIG.h \ - liboctave/numeric/floatLU.h \ liboctave/numeric/floatQR.h \ liboctave/numeric/floatQRP.h \ liboctave/numeric/floatSVD.h \ @@ -55,6 +50,7 @@ liboctave/numeric/hess.h \ liboctave/numeric/lo-mappers.h \ liboctave/numeric/lo-specfun.h \ + liboctave/numeric/lu.h \ liboctave/numeric/oct-convn.h \ liboctave/numeric/oct-fftw.h \ liboctave/numeric/oct-norm.h \ @@ -75,7 +71,6 @@ liboctave/numeric/randpoisson.c NUMERIC_SRC = \ - liboctave/numeric/CmplxLU.cc \ liboctave/numeric/CmplxQR.cc \ liboctave/numeric/CmplxQRP.cc \ liboctave/numeric/CmplxSVD.cc \ @@ -89,17 +84,14 @@ liboctave/numeric/Quad.cc \ liboctave/numeric/aepbalance.cc \ liboctave/numeric/chol.cc \ - liboctave/numeric/dbleLU.cc \ liboctave/numeric/dbleQR.cc \ liboctave/numeric/dbleQRP.cc \ liboctave/numeric/dbleSVD.cc \ liboctave/numeric/eigs-base.cc \ - liboctave/numeric/fCmplxLU.cc \ liboctave/numeric/fCmplxQR.cc \ liboctave/numeric/fCmplxQRP.cc \ liboctave/numeric/fCmplxSVD.cc \ liboctave/numeric/fEIG.cc \ - liboctave/numeric/floatLU.cc \ liboctave/numeric/floatQR.cc \ liboctave/numeric/floatQRP.cc \ liboctave/numeric/floatSVD.cc \ @@ -107,6 +99,7 @@ liboctave/numeric/hess.cc \ liboctave/numeric/lo-mappers.cc \ liboctave/numeric/lo-specfun.cc \ + liboctave/numeric/lu.cc \ liboctave/numeric/oct-convn.cc \ liboctave/numeric/oct-fftw.cc \ liboctave/numeric/oct-norm.cc \ @@ -120,7 +113,6 @@ $(NUMERIC_C_SRC) LIBOCTAVE_TEMPLATE_SRC += \ - liboctave/numeric/base-lu.cc \ liboctave/numeric/base-qr.cc \ liboctave/numeric/bsxfun-defs.cc
--- a/liboctave/operators/mx-defs.h Tue Feb 16 11:13:55 2016 -0500 +++ b/liboctave/operators/mx-defs.h Tue Feb 16 12:58:32 2016 -0500 @@ -75,10 +75,7 @@ class FloatSVD; class FloatComplexSVD; -class LU; -class ComplexLU; -class FloatLU; -class FloatComplexLU; +template <typename T> class lu; class QR; class ComplexQR;
--- a/liboctave/operators/mx-ext.h Tue Feb 16 11:13:55 2016 -0500 +++ b/liboctave/operators/mx-ext.h Tue Feb 16 12:58:32 2016 -0500 @@ -62,10 +62,7 @@ // Result of an LU decomposition. -#include "dbleLU.h" -#include "CmplxLU.h" -#include "floatLU.h" -#include "fCmplxLU.h" +#include "lu.h" // Result of a QR decomposition.