From 273f80384626c16b716abf55a96ebe05f856ea8f Mon Sep 17 00:00:00 2001 From: Melven Roehrig-Zoellner Date: Fri, 9 Dec 2022 18:50:12 +0000 Subject: [PATCH] Add BDCSVD_LAPACKE binding --- Eigen/SVD | 5 +- Eigen/src/SVD/BDCSVD.h | 2 +- Eigen/src/SVD/BDCSVD_LAPACKE.h | 163 ++++++++++++++++++++++++++++++ Eigen/src/SVD/JacobiSVD_LAPACKE.h | 12 ++- Eigen/src/misc/lapacke_helpers.h | 1 + doc/UsingBlasLapackBackends.dox | 6 ++ 6 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 Eigen/src/SVD/BDCSVD_LAPACKE.h diff --git a/Eigen/SVD b/Eigen/SVD index 345179496..8241c7386 100644 --- a/Eigen/SVD +++ b/Eigen/SVD @@ -36,14 +36,17 @@ #include "src/SVD/SVDBase.h" #include "src/SVD/JacobiSVD.h" #include "src/SVD/BDCSVD.h" -#if defined(EIGEN_USE_LAPACKE) && !defined(EIGEN_USE_LAPACKE_STRICT) +#ifdef EIGEN_USE_LAPACKE #ifdef EIGEN_USE_MKL #include "mkl_lapacke.h" #else #include "src/misc/lapacke.h" #endif +#ifndef EIGEN_USE_LAPACKE_STRICT #include "src/SVD/JacobiSVD_LAPACKE.h" #endif +#include "src/SVD/BDCSVD_LAPACKE.h" +#endif #include "src/Core/util/ReenableStupidWarnings.h" diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h index 1b10c0e42..a69fdcabb 100644 --- a/Eigen/src/SVD/BDCSVD.h +++ b/Eigen/src/SVD/BDCSVD.h @@ -236,7 +236,6 @@ public: } private: - void allocate(Index rows, Index cols, unsigned int computationOptions); BDCSVD& compute_impl(const MatrixType& matrix, unsigned int computationOptions); void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift); void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V); @@ -254,6 +253,7 @@ private: void computeBaseCase(SVDType& svd, Index n, Index firstCol, Index firstRowW, Index firstColW, Index shift); protected: + void allocate(Index rows, Index cols, unsigned int computationOptions); MatrixXr m_naiveU, m_naiveV; MatrixXr m_computed; Index m_nRec; diff --git a/Eigen/src/SVD/BDCSVD_LAPACKE.h b/Eigen/src/SVD/BDCSVD_LAPACKE.h new file mode 100644 index 000000000..d4cc173c4 --- /dev/null +++ b/Eigen/src/SVD/BDCSVD_LAPACKE.h @@ -0,0 +1,163 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2022 Melven Roehrig-Zoellner +// Copyright (c) 2011, Intel Corporation. All rights reserved. +// +// This file is based on the JacobiSVD_LAPACKE.h originally from Intel - +// see license notice below: +/* + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Singular Value Decomposition - SVD (divide and conquer variant) + ******************************************************************************** +*/ +#ifndef EIGEN_BDCSVD_LAPACKE_H +#define EIGEN_BDCSVD_LAPACKE_H + +namespace Eigen { + +namespace internal { + +namespace lapacke_helpers { + +/** \internal Specialization for the data types supported by LAPACKe */ + +// defining a derived class to allow access to protected members +template +class BDCSVD_LAPACKE : public BDCSVD { + typedef BDCSVD SVD; + typedef typename SVD::MatrixType MatrixType; + typedef typename SVD::Scalar Scalar; + typedef typename SVD::RealScalar RealScalar; + +public: + // construct this by moving from a parent object + BDCSVD_LAPACKE(SVD&& svd) : SVD(std::move(svd)) {} + + void compute_impl_lapacke(const MatrixType& matrix, unsigned int computationOptions) { + + SVD::allocate(matrix.rows(), matrix.cols(), computationOptions); + + SVD::m_nonzeroSingularValues = SVD::m_diagSize; + + // prepare arguments to ?gesdd + const lapack_int matrix_order = lapack_storage_of(matrix); + const char jobz = (SVD::m_computeFullU || SVD::m_computeFullV) ? 'A' : (SVD::m_computeThinU || SVD::m_computeThinV) ? 'S' : 'N'; + const lapack_int u_cols = (jobz == 'A') ? to_lapack(SVD::m_rows) : (jobz == 'S') ? to_lapack(SVD::m_diagSize) : 1; + const lapack_int vt_rows = (jobz == 'A') ? to_lapack(SVD::m_cols) : (jobz == 'S') ? to_lapack(SVD::m_diagSize) : 1; + lapack_int ldu, ldvt; + Scalar *u, *vt, dummy; + MatrixType localU; + if (SVD::computeU() && !(SVD::m_computeThinU && SVD::m_computeFullV) ) { + ldu = to_lapack(SVD::m_matrixU.outerStride()); + u = SVD::m_matrixU.data(); + } else if (SVD::computeV()) { + localU.resize(SVD::m_rows, u_cols); + ldu = to_lapack(localU.outerStride()); + u = localU.data(); + } else { ldu=1; u=&dummy; } + MatrixType localV; + if (SVD::computeU() || SVD::computeV()) { + localV.resize(vt_rows, SVD::m_cols); + ldvt = to_lapack(localV.outerStride()); + vt = localV.data(); + } else { ldvt=1; vt=&dummy; } + MatrixType temp; temp = matrix; + + // actual call to ?gesdd + lapack_int info = gesdd( matrix_order, jobz, to_lapack(SVD::m_rows), to_lapack(SVD::m_cols), + to_lapack(temp.data()), to_lapack(temp.outerStride()), (RealScalar*)SVD::m_singularValues.data(), + to_lapack(u), ldu, to_lapack(vt), ldvt); + + // Check the result of the LAPACK call + if (info < 0 || !SVD::m_singularValues.allFinite()) { + // this includes info == -4 => NaN entry in A + SVD::m_info = InvalidInput; + } else if (info > 0 ) { + SVD::m_info = NoConvergence; + } else { + SVD::m_info = Success; + if (SVD::m_computeThinU && SVD::m_computeFullV) { + SVD::m_matrixU = localU.leftCols(SVD::m_matrixU.cols()); + } + if (SVD::computeV()) { + SVD::m_matrixV = localV.adjoint().leftCols(SVD::m_matrixV.cols()); + } + } + SVD::m_isInitialized = true; + } +}; + +template +BDCSVD& BDCSVD_wrapper(BDCSVD& svd, const MatrixType_& matrix, int computationOptions) +{ + // we need to move to the wrapper type and back + BDCSVD_LAPACKE tmpSvd(std::move(svd)); + tmpSvd.compute_impl_lapacke(matrix, computationOptions); + svd = std::move(tmpSvd); + return svd; +} + +} // end namespace lapacke_helpers + +} // end namespace internal + +#define EIGEN_LAPACKE_SDD(EIGTYPE, EIGCOLROW, OPTIONS) \ +template<> inline \ +BDCSVD, OPTIONS>& \ +BDCSVD, OPTIONS>::compute_impl(const Matrix& matrix, unsigned int computationOptions) {\ + return internal::lapacke_helpers::BDCSVD_wrapper(*this, matrix, computationOptions); \ +} + +#define EIGEN_LAPACK_SDD_OPTIONS(OPTIONS) \ + EIGEN_LAPACKE_SDD(double, ColMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(float, ColMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(dcomplex, ColMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(scomplex, ColMajor, OPTIONS) \ +\ + EIGEN_LAPACKE_SDD(double, RowMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(float, RowMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(dcomplex, RowMajor, OPTIONS) \ + EIGEN_LAPACKE_SDD(scomplex, RowMajor, OPTIONS) + +EIGEN_LAPACK_SDD_OPTIONS(0) +EIGEN_LAPACK_SDD_OPTIONS(ComputeThinU) +EIGEN_LAPACK_SDD_OPTIONS(ComputeThinV) +EIGEN_LAPACK_SDD_OPTIONS(ComputeFullU) +EIGEN_LAPACK_SDD_OPTIONS(ComputeFullV) +EIGEN_LAPACK_SDD_OPTIONS(ComputeThinU | ComputeThinV) +EIGEN_LAPACK_SDD_OPTIONS(ComputeFullU | ComputeFullV) +EIGEN_LAPACK_SDD_OPTIONS(ComputeThinU | ComputeFullV) +EIGEN_LAPACK_SDD_OPTIONS(ComputeFullU | ComputeThinV) + +#undef EIGEN_LAPACK_SDD_OPTIONS + +#undef EIGEN_LAPACKE_SDD + +} // end namespace Eigen + +#endif // EIGEN_BDCSVD_LAPACKE_H diff --git a/Eigen/src/SVD/JacobiSVD_LAPACKE.h b/Eigen/src/SVD/JacobiSVD_LAPACKE.h index c5717e9fd..93244cd98 100644 --- a/Eigen/src/SVD/JacobiSVD_LAPACKE.h +++ b/Eigen/src/SVD/JacobiSVD_LAPACKE.h @@ -72,8 +72,16 @@ JacobiSVD, OPTION } else { ldvt=1; vt=&dummy; }\ Matrix superb; superb.resize(m_diagSize, 1); \ MatrixType m_temp; m_temp = matrix; \ - LAPACKE_##LAPACKE_PREFIX##gesvd( matrix_order, jobu, jobvt, internal::convert_index(m_rows), internal::convert_index(m_cols), (LAPACKE_TYPE*)m_temp.data(), lda, (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \ - if (computeV()) m_matrixV = localV.adjoint(); \ + lapack_int info = LAPACKE_##LAPACKE_PREFIX##gesvd( matrix_order, jobu, jobvt, internal::convert_index(m_rows), internal::convert_index(m_cols), (LAPACKE_TYPE*)m_temp.data(), lda, (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \ + /* Check the result of the LAPACK call */ \ + if (info < 0 || !m_singularValues.allFinite()) { \ + m_info = InvalidInput; \ + } else if (info > 0 ) { \ + m_info = NoConvergence; \ + } else { \ + m_info = Success; \ + if (computeV()) m_matrixV = localV.adjoint(); \ + } \ /* for(int i=0;i\code ?gesvd \endcode +Singular value decomposition \n \c EIGEN_USE_LAPACKE \n \c EIGEN_USE_LAPACKE_STRICT \code +BDCSVD svd; +svd.compute(m1); +\endcode\code +?gesdd +\endcode Eigen-value decompositions \n \c EIGEN_USE_LAPACKE \n \c EIGEN_USE_LAPACKE_STRICT \code EigenSolver es(m1); ComplexEigenSolver ces(m1);