Mercurial > octave-libtiff
view liboctave/numeric/lu.cc @ 30092:51040836259d
maint: rename "m_l" to "m_L" for clarity in class lu.
* lu.cc, lu.h: rename "m_l" to "m_L" for clarity in class lu.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 31 Aug 2021 12:54:59 -0700 |
parents | 43622407af81 |
children | 796f54d4ddbf |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1996-2021 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // 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 // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <algorithm> #include "CColVector.h" #include "CMatrix.h" #include "PermMatrix.h" #include "dColVector.h" #include "dMatrix.h" #include "fCColVector.h" #include "fCMatrix.h" #include "fColVector.h" #include "fMatrix.h" #include "lo-error.h" #include "lo-lapack-proto.h" #include "lo-qrupdate-proto.h" #include "lu.h" #include "oct-locbuf.h" namespace octave { namespace math { // FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type> // but m_ipvt is an Array<octave_f77_int_type>. This could cause // trouble for large arrays if octave_f77_int_type is 32-bits but // octave_idx_type is 64. Since this constructor is called from // Fluupdate, it could be given values that are out of range. We // should ensure that the values are within range here. template <typename T> lu<T>::lu (const T& l, const T& u, const PermMatrix& p) : m_a_fact (u), m_L (l), m_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 m_L.dims () == dim_vector (); } template <typename T> void lu<T>::unpack (void) { if (packed ()) { m_L = L (); m_a_fact = U (); // FIXME: sub-optimal // FIXME: getp returns Array<octave_idx_type> but m_ipvt is // Array<octave_f77_int_type>. However, getp produces its // result from a valid m_ipvt array so validation should not be // necessary. OTOH, it might be better to have a version of // getp that doesn't cause us to convert from // Array<octave_f77_int_type> to Array<octave_idx_type> and // back again. m_ipvt = getp (); } } template <typename T> T lu<T>::L (void) const { if (packed ()) { octave_idx_type a_nr = m_a_fact.rows (); octave_idx_type a_nc = m_a_fact.columns (); 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) = m_a_fact.xelem (i, j); } return l; } else return m_L; } template <typename T> T lu<T>::U (void) const { if (packed ()) { octave_idx_type a_nr = m_a_fact.rows (); octave_idx_type a_nc = m_a_fact.columns (); 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) = m_a_fact.xelem (i, j); } return u; } else return m_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 m_a_fact; } template <typename T> Array<octave_idx_type> lu<T>::getp (void) const { if (packed ()) { octave_idx_type a_nr = m_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 < m_ipvt.numel (); i++) { octave_idx_type k = m_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 m_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 = m_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 (m_a_fact.rows (), m_a_fact.columns ()); for (octave_idx_type i = 0; i < k; i++) { if (m_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 <> OCTAVE_API lu<Matrix>::lu (const Matrix& a) { F77_INT a_nr = to_f77_int (a.rows ()); F77_INT a_nc = to_f77_int (a.columns ()); F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); m_ipvt.resize (dim_vector (mn, 1)); F77_INT *pipvt = m_ipvt.fortran_vec (); m_a_fact = a; double *tmp_data = m_a_fact.fortran_vec (); F77_INT info = 0; F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); for (F77_INT i = 0; i < mn; i++) pipvt[i] -= 1; } #if defined (HAVE_QRUPDATE_LUU) template <> OCTAVE_API void lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v) { if (packed ()) unpack (); Matrix& l = m_L; Matrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != 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 <> OCTAVE_API void lu<Matrix>::update (const Matrix& u, const Matrix& v) { if (packed ()) unpack (); Matrix& l = m_L; Matrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); for (volatile F77_INT i = 0; i < u_nc; 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 <> OCTAVE_API void lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v) { if (packed ()) unpack (); Matrix& l = m_L; Matrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); ColumnVector utmp = u; ColumnVector vtmp = v; OCTAVE_LOCAL_BUFFER (double, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, m_ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } template <> OCTAVE_API void lu<Matrix>::update_piv (const Matrix& u, const Matrix& v) { if (packed ()) unpack (); Matrix& l = m_L; Matrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); OCTAVE_LOCAL_BUFFER (double, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment for (volatile F77_INT i = 0; i < u_nc; 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, m_ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); } for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } #endif template <> OCTAVE_API lu<FloatMatrix>::lu (const FloatMatrix& a) { F77_INT a_nr = to_f77_int (a.rows ()); F77_INT a_nc = to_f77_int (a.columns ()); F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); m_ipvt.resize (dim_vector (mn, 1)); F77_INT *pipvt = m_ipvt.fortran_vec (); m_a_fact = a; float *tmp_data = m_a_fact.fortran_vec (); F77_INT info = 0; F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); for (F77_INT i = 0; i < mn; i++) pipvt[i] -= 1; } #if defined (HAVE_QRUPDATE_LUU) template <> OCTAVE_API void lu<FloatMatrix>::update (const FloatColumnVector& u, const FloatColumnVector& v) { if (packed ()) unpack (); FloatMatrix& l = m_L; FloatMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != 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 <> OCTAVE_API void lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v) { if (packed ()) unpack (); FloatMatrix& l = m_L; FloatMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); for (volatile F77_INT i = 0; i < u_nc; 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 <> OCTAVE_API void lu<FloatMatrix>::update_piv (const FloatColumnVector& u, const FloatColumnVector& v) { if (packed ()) unpack (); FloatMatrix& l = m_L; FloatMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); FloatColumnVector utmp = u; FloatColumnVector vtmp = v; OCTAVE_LOCAL_BUFFER (float, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, m_ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } template <> OCTAVE_API void lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v) { if (packed ()) unpack (); FloatMatrix& l = m_L; FloatMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); OCTAVE_LOCAL_BUFFER (float, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment for (volatile F77_INT i = 0; i < u_nc; 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, m_ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); } for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } #endif template <> OCTAVE_API lu<ComplexMatrix>::lu (const ComplexMatrix& a) { F77_INT a_nr = to_f77_int (a.rows ()); F77_INT a_nc = to_f77_int (a.columns ()); F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); m_ipvt.resize (dim_vector (mn, 1)); F77_INT *pipvt = m_ipvt.fortran_vec (); m_a_fact = a; Complex *tmp_data = m_a_fact.fortran_vec (); F77_INT info = 0; F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data), a_nr, pipvt, info)); for (F77_INT i = 0; i < mn; i++) pipvt[i] -= 1; } #if defined (HAVE_QRUPDATE_LUU) template <> OCTAVE_API void lu<ComplexMatrix>::update (const ComplexColumnVector& u, const ComplexColumnVector& v) { if (packed ()) unpack (); ComplexMatrix& l = m_L; ComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); ComplexColumnVector utmp = u; ComplexColumnVector vtmp = v; F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()))); } template <> OCTAVE_API void lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v) { if (packed ()) unpack (); ComplexMatrix& l = m_L; ComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); for (volatile F77_INT i = 0; i < u_nc; i++) { ComplexColumnVector utmp = u.column (i); ComplexColumnVector vtmp = v.column (i); F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()))); } } template <> OCTAVE_API void lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u, const ComplexColumnVector& v) { if (packed ()) unpack (); ComplexMatrix& l = m_L; ComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); ComplexColumnVector utmp = u; ComplexColumnVector vtmp = v; OCTAVE_LOCAL_BUFFER (Complex, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, m_ipvt.fortran_vec (), F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()), F77_DBLE_CMPLX_ARG (w))); for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } template <> OCTAVE_API void lu<ComplexMatrix>::update_piv (const ComplexMatrix& u, const ComplexMatrix& v) { if (packed ()) unpack (); ComplexMatrix& l = m_L; ComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); OCTAVE_LOCAL_BUFFER (Complex, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment for (volatile F77_INT i = 0; i < u_nc; i++) { ComplexColumnVector utmp = u.column (i); ComplexColumnVector vtmp = v.column (i); F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, m_ipvt.fortran_vec (), F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()), F77_DBLE_CMPLX_ARG (w))); } for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } #endif template <> OCTAVE_API lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a) { F77_INT a_nr = to_f77_int (a.rows ()); F77_INT a_nc = to_f77_int (a.columns ()); F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); m_ipvt.resize (dim_vector (mn, 1)); F77_INT *pipvt = m_ipvt.fortran_vec (); m_a_fact = a; FloatComplex *tmp_data = m_a_fact.fortran_vec (); F77_INT info = 0; F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr, pipvt, info)); for (F77_INT i = 0; i < mn; i++) pipvt[i] -= 1; } #if defined (HAVE_QRUPDATE_LUU) template <> OCTAVE_API void lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v) { if (packed ()) unpack (); FloatComplexMatrix& l = m_L; FloatComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); FloatComplexColumnVector utmp = u; FloatComplexColumnVector vtmp = v; F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m, F77_CMPLX_ARG (r.fortran_vec ()), k, F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ()))); } template <> OCTAVE_API void lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v) { if (packed ()) unpack (); FloatComplexMatrix& l = m_L; FloatComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); for (volatile F77_INT i = 0; i < u_nc; i++) { FloatComplexColumnVector utmp = u.column (i); FloatComplexColumnVector vtmp = v.column (i); F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m, F77_CMPLX_ARG (r.fortran_vec ()), k, F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ()))); } } template <> OCTAVE_API void lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v) { if (packed ()) unpack (); FloatComplexMatrix& l = m_L; FloatComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nel = to_f77_int (u.numel ()); F77_INT v_nel = to_f77_int (v.numel ()); if (u_nel != m || v_nel != n) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); FloatComplexColumnVector utmp = u; FloatComplexColumnVector vtmp = v; OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m, F77_CMPLX_ARG (r.fortran_vec ()), k, m_ipvt.fortran_vec (), F77_CONST_CMPLX_ARG (utmp.data ()), F77_CONST_CMPLX_ARG (vtmp.data ()), F77_CMPLX_ARG (w))); for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } template <> OCTAVE_API void lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u, const FloatComplexMatrix& v) { if (packed ()) unpack (); FloatComplexMatrix& l = m_L; FloatComplexMatrix& r = m_a_fact; F77_INT m = to_f77_int (l.rows ()); F77_INT n = to_f77_int (r.columns ()); F77_INT k = to_f77_int (l.columns ()); F77_INT u_nr = to_f77_int (u.rows ()); F77_INT u_nc = to_f77_int (u.columns ()); F77_INT v_nr = to_f77_int (v.rows ()); F77_INT v_nc = to_f77_int (v.columns ()); if (u_nr != m || v_nr != n || u_nc != v_nc) (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment for (volatile F77_INT i = 0; i < u_nc; i++) { FloatComplexColumnVector utmp = u.column (i); FloatComplexColumnVector vtmp = v.column (i); F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m, F77_CMPLX_ARG (r.fortran_vec ()), k, m_ipvt.fortran_vec (), F77_CONST_CMPLX_ARG (utmp.data ()), F77_CONST_CMPLX_ARG (vtmp.data ()), F77_CMPLX_ARG (w))); } for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement } #endif // Instantiations we need. template class lu<Matrix>; template class lu<FloatMatrix>; template class lu<ComplexMatrix>; template class lu<FloatComplexMatrix>; } }