aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/SVD
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/SVD')
-rw-r--r--Eigen/src/SVD/BDCSVD.h61
-rw-r--r--Eigen/src/SVD/CMakeLists.txt6
-rw-r--r--Eigen/src/SVD/JacobiSVD.h63
-rw-r--r--Eigen/src/SVD/JacobiSVD_LAPACKE.h (renamed from Eigen/src/SVD/JacobiSVD_MKL.h)48
-rw-r--r--Eigen/src/SVD/SVDBase.h3
5 files changed, 86 insertions, 95 deletions
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h
index 3552c87bf..25fca6f4d 100644
--- a/Eigen/src/SVD/BDCSVD.h
+++ b/Eigen/src/SVD/BDCSVD.h
@@ -11,7 +11,7 @@
// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
// Copyright (C) 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
-// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2014-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -21,6 +21,7 @@
#define EIGEN_BDCSVD_H
// #define EIGEN_BDCSVD_DEBUG_VERBOSE
// #define EIGEN_BDCSVD_SANITY_CHECKS
+
namespace Eigen {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
@@ -49,6 +50,18 @@ struct traits<BDCSVD<_MatrixType> >
*
* \tparam _MatrixType the type of the matrix of which we are computing the SVD decomposition
*
+ * This class first reduces the input matrix to bi-diagonal form using class UpperBidiagonalization,
+ * and then performs a divide-and-conquer diagonalization. Small blocks are diagonalized using class JacobiSVD.
+ * You can control the switching size with the setSwitchSize() method, default is 16.
+ * For small matrice (<16), it is thus preferable to directly use JacobiSVD. For larger ones, BDCSVD is highly
+ * recommended and can several order of magnitude faster.
+ *
+ * \warning this algorithm is unlikely to provide accurate result when compiled with unsafe math optimizations.
+ * For instance, this concerns Intel's compiler (ICC), which perfroms such optimization by default unless
+ * you compile with the \c -fp-model \c precise option. Likewise, the \c -ffast-math option of GCC or clang will
+ * significantly degrade the accuracy.
+ *
+ * \sa class JacobiSVD
*/
template<typename _MatrixType>
class BDCSVD : public SVDBase<BDCSVD<_MatrixType> >
@@ -228,6 +241,8 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
#endif
allocate(matrix.rows(), matrix.cols(), computationOptions);
using std::abs;
+
+ const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
//**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
if(matrix.cols() < m_algoswap)
@@ -266,7 +281,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
{
RealScalar a = abs(m_computed.coeff(i, i));
m_singularValues.coeffRef(i) = a * scale;
- if (a == 0)
+ if (a<considerZero)
{
m_nonzeroSingularValues = i;
m_singularValues.tail(m_diagSize - i - 1).setZero();
@@ -380,6 +395,7 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
using std::abs;
const Index n = lastCol - firstCol + 1;
const Index k = n/2;
+ const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar alphaK;
RealScalar betaK;
RealScalar r0;
@@ -434,7 +450,7 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
}
if (m_compV) m_naiveV(firstRowW+k, firstColW) = 1;
- if (r0 == 0)
+ if (r0<considerZero)
{
c0 = 1;
s0 = 0;
@@ -553,6 +569,8 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
template <typename MatrixType>
void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
{
+ const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
+ using std::abs;
ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
ArrayRef diag = m_workspace.head(n);
@@ -575,7 +593,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
while(actual_n>1 && diag(actual_n-1)==0) --actual_n;
Index m = 0; // size of the deflated problem
for(Index k=0;k<actual_n;++k)
- if(col0(k)!=0)
+ if(abs(col0(k))>considerZero)
m_workspaceI(m++) = k;
Map<ArrayXi> perm(m_workspaceI.data(),m);
@@ -600,7 +618,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
{
Index actual_n = n;
- while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
+ while(actual_n>1 && abs(col0(actual_n-1))<considerZero) --actual_n;
std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\n";
std::cout << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
std::cout << " check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
@@ -680,6 +698,7 @@ typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar
res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu));
}
return res;
+
}
template <typename MatrixType>
@@ -746,14 +765,14 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar muPrev, muCur;
if (shift == left)
{
- muPrev = (right - left) * 0.1;
+ muPrev = (right - left) * RealScalar(0.1);
if (k == actual_n-1) muCur = right - left;
- else muCur = (right - left) * 0.5;
+ else muCur = (right - left) * RealScalar(0.5);
}
else
{
- muPrev = -(right - left) * 0.1;
- muCur = -(right - left) * 0.5;
+ muPrev = -(right - left) * RealScalar(0.1);
+ muCur = -(right - left) * RealScalar(0.5);
}
RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
@@ -798,15 +817,15 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar leftShifted, rightShifted;
if (shift == left)
{
- leftShifted = RealScalar(1)/NumTraits<RealScalar>::highest();
+ leftShifted = (std::numeric_limits<RealScalar>::min)();
// I don't understand why the case k==0 would be special there:
// if (k == 0) rightShifted = right - left; else
- rightShifted = (k==actual_n-1) ? right : ((right - left) * 0.6); // theoretically we can take 0.5, but let's be safe
+ rightShifted = (k==actual_n-1) ? right : ((right - left) * RealScalar(0.6)); // theoretically we can take 0.5, but let's be safe
}
else
{
- leftShifted = -(right - left) * 0.6;
- rightShifted = -RealScalar(1)/NumTraits<RealScalar>::highest();
+ leftShifted = -(right - left) * RealScalar(0.6);
+ rightShifted = -(std::numeric_limits<RealScalar>::min)();
}
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
@@ -817,7 +836,10 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
if(!(fLeft * fRight<0))
+ {
+ std::cout << "fLeft: " << leftShifted << " - " << diagShifted.head(10).transpose() << "\n ; " << bool(left==shift) << " " << (left-shift) << "\n";
std::cout << k << " : " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; " << left << " - " << right << " -> " << leftShifted << " " << rightShifted << " shift=" << shift << "\n";
+ }
#endif
eigen_internal_assert(fLeft * fRight < 0);
@@ -1028,8 +1050,9 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
Diagonal<MatrixXr> fulldiag(m_computed);
VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
+ const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
- RealScalar epsilon_strict = NumTraits<RealScalar>::epsilon() * maxDiag;
+ RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
@@ -1082,7 +1105,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
{
// Check for total deflation
// If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting
- bool total_deflation = (col0.tail(length-1).array()==RealScalar(0)).all();
+ bool total_deflation = (col0.tail(length-1).array()<considerZero).all();
// Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge.
// First, compute the respective permutation.
@@ -1093,7 +1116,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
// Move deflated diagonal entries at the end.
for(Index i=1; i<length; ++i)
- if(diag(i)==0)
+ if(abs(diag(i))<considerZero)
permutation[p++] = i;
Index i=1, j=k+1;
@@ -1112,7 +1135,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
for(Index i=1; i<length; ++i)
{
Index pi = permutation[i];
- if(diag(pi)==0 || diag(0)<diag(pi))
+ if(abs(diag(pi))<considerZero || diag(0)<diag(pi))
permutation[i-1] = permutation[i];
else
{
@@ -1163,7 +1186,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
//condition 4.4
{
Index i = length-1;
- while(i>0 && (diag(i)==0 || col0(i)==0)) --i;
+ while(i>0 && (abs(diag(i))<considerZero || abs(col0(i))<considerZero)) --i;
for(; i>1;--i)
if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
{
@@ -1177,7 +1200,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
for(Index j=2;j<length;++j)
- assert(diag(j-1)<=diag(j) || diag(j)==0);
+ assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
#endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
diff --git a/Eigen/src/SVD/CMakeLists.txt b/Eigen/src/SVD/CMakeLists.txt
deleted file mode 100644
index 55efc44b1..000000000
--- a/Eigen/src/SVD/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_SVD_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_SVD_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SVD COMPONENT Devel
- )
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index 1940c8294..ea2bd62eb 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -419,38 +419,6 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
}
};
-template<typename MatrixType, typename RealScalar, typename Index>
-void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
- JacobiRotation<RealScalar> *j_left,
- JacobiRotation<RealScalar> *j_right)
-{
- using std::sqrt;
- using std::abs;
- Matrix<RealScalar,2,2> m;
- m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
- numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
- JacobiRotation<RealScalar> rot1;
- RealScalar t = m.coeff(0,0) + m.coeff(1,1);
- RealScalar d = m.coeff(1,0) - m.coeff(0,1);
- if(d == RealScalar(0))
- {
- rot1.s() = RealScalar(0);
- rot1.c() = RealScalar(1);
- }
- else
- {
- // If d!=0, then t/d cannot overflow because the magnitude of the
- // entries forming d are not too small compared to the ones forming t.
- RealScalar u = t / d;
- RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u));
- rot1.s() = RealScalar(1) / tmp;
- rot1.c() = u / tmp;
- }
- m.applyOnTheLeft(0,1,rot1);
- j_right->makeJacobi(m,0,1);
- *j_left = rot1 * j_right->transpose();
-}
-
template<typename _MatrixType, int QRPreconditioner>
struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
{
@@ -697,10 +665,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
// only worsening the precision of U and V as we accumulate more rotations
const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
- // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
- // FIXME What about considerering any denormal numbers as zero, using:
- // const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
- const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
+ // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
+ const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
// Scaling factor to reduce over/under-flows
RealScalar scale = matrix.cwiseAbs().maxCoeff();
@@ -745,7 +711,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
{
finished = false;
// perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
- // the complex to real operation returns true is the updated 2x2 block is not already diagonal
+ // the complex to real operation returns true if the updated 2x2 block is not already diagonal
if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
{
JacobiRotation<RealScalar> j_left, j_right;
@@ -759,7 +725,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
// keep track of the largest diagonal coefficient
- maxDiagEntry = numext::maxi(maxDiagEntry,numext::maxi(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
+ maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
}
}
}
@@ -770,9 +736,22 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
for(Index i = 0; i < m_diagSize; ++i)
{
- RealScalar a = abs(m_workMatrix.coeff(i,i));
- m_singularValues.coeffRef(i) = a;
- if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
+ // For a complex matrix, some diagonal coefficients might note have been
+ // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
+ // of some diagonal entry might not be null.
+ if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
+ {
+ RealScalar a = abs(m_workMatrix.coeff(i,i));
+ m_singularValues.coeffRef(i) = abs(a);
+ if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
+ }
+ else
+ {
+ // m_workMatrix.coeff(i,i) is already real, no difficulty:
+ RealScalar a = numext::real(m_workMatrix.coeff(i,i));
+ m_singularValues.coeffRef(i) = abs(a);
+ if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
+ }
}
m_singularValues *= scale;
@@ -802,7 +781,6 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
return *this;
}
-#ifndef __CUDACC__
/** \svd_module
*
* \return the singular value decomposition of \c *this computed by two-sided
@@ -816,7 +794,6 @@ MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
{
return JacobiSVD<PlainObject>(*this, computationOptions);
}
-#endif // __CUDACC__
} // end namespace Eigen
diff --git a/Eigen/src/SVD/JacobiSVD_MKL.h b/Eigen/src/SVD/JacobiSVD_LAPACKE.h
index 14e461c4e..50272154f 100644
--- a/Eigen/src/SVD/JacobiSVD_MKL.h
+++ b/Eigen/src/SVD/JacobiSVD_LAPACKE.h
@@ -25,21 +25,19 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Singular Value Decomposition - SVD.
********************************************************************************
*/
-#ifndef EIGEN_JACOBISVD_MKL_H
-#define EIGEN_JACOBISVD_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
+#ifndef EIGEN_JACOBISVD_LAPACKE_H
+#define EIGEN_JACOBISVD_LAPACKE_H
namespace Eigen {
-/** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_SVD(EIGTYPE, MKLTYPE, MKLRTYPE, MKLPREFIX, EIGCOLROW, MKLCOLROW) \
+#define EIGEN_LAPACKE_SVD(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \
template<> inline \
JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, ColPivHouseholderQRPreconditioner>& \
JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, ColPivHouseholderQRPreconditioner>::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix, unsigned int computationOptions) \
@@ -52,41 +50,41 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, ColPiv
/*const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();*/ \
m_nonzeroSingularValues = m_diagSize; \
\
- lapack_int lda = matrix.outerStride(), ldu, ldvt; \
- lapack_int matrix_order = MKLCOLROW; \
+ lapack_int lda = internal::convert_index<lapack_int>(matrix.outerStride()), ldu, ldvt; \
+ lapack_int matrix_order = LAPACKE_COLROW; \
char jobu, jobvt; \
- MKLTYPE *u, *vt, dummy; \
+ LAPACKE_TYPE *u, *vt, dummy; \
jobu = (m_computeFullU) ? 'A' : (m_computeThinU) ? 'S' : 'N'; \
jobvt = (m_computeFullV) ? 'A' : (m_computeThinV) ? 'S' : 'N'; \
if (computeU()) { \
- ldu = m_matrixU.outerStride(); \
- u = (MKLTYPE*)m_matrixU.data(); \
+ ldu = internal::convert_index<lapack_int>(m_matrixU.outerStride()); \
+ u = (LAPACKE_TYPE*)m_matrixU.data(); \
} else { ldu=1; u=&dummy; }\
MatrixType localV; \
- ldvt = (m_computeFullV) ? m_cols : (m_computeThinV) ? m_diagSize : 1; \
+ ldvt = (m_computeFullV) ? internal::convert_index<lapack_int>(m_cols) : (m_computeThinV) ? internal::convert_index<lapack_int>(m_diagSize) : 1; \
if (computeV()) { \
localV.resize(ldvt, m_cols); \
- vt = (MKLTYPE*)localV.data(); \
+ vt = (LAPACKE_TYPE*)localV.data(); \
} else { ldvt=1; vt=&dummy; }\
- Matrix<MKLRTYPE, Dynamic, Dynamic> superb; superb.resize(m_diagSize, 1); \
+ Matrix<LAPACKE_RTYPE, Dynamic, Dynamic> superb; superb.resize(m_diagSize, 1); \
MatrixType m_temp; m_temp = matrix; \
- LAPACKE_##MKLPREFIX##gesvd( matrix_order, jobu, jobvt, m_rows, m_cols, (MKLTYPE*)m_temp.data(), lda, (MKLRTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \
+ LAPACKE_##LAPACKE_PREFIX##gesvd( matrix_order, jobu, jobvt, internal::convert_index<lapack_int>(m_rows), internal::convert_index<lapack_int>(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(); \
/* for(int i=0;i<m_diagSize;i++) if (m_singularValues.coeffRef(i) < precision) { m_nonzeroSingularValues--; m_singularValues.coeffRef(i)=RealScalar(0);}*/ \
m_isInitialized = true; \
return *this; \
}
-EIGEN_MKL_SVD(double, double, double, d, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SVD(float, float, float , s, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SVD(dcomplex, MKL_Complex16, double, z, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SVD(scomplex, MKL_Complex8, float , c, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SVD(double, double, double, d, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SVD(float, float, float , s, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SVD(double, double, double, d, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_SVD(float, float, float , s, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_SVD(dcomplex, MKL_Complex16, double, z, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_SVD(scomplex, MKL_Complex8, float , c, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SVD(double, double, double, d, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SVD(float, float, float , s, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, RowMajor, LAPACK_ROW_MAJOR)
} // end namespace Eigen
-#endif // EIGEN_JACOBISVD_MKL_H
+#endif // EIGEN_JACOBISVD_LAPACKE_H
diff --git a/Eigen/src/SVD/SVDBase.h b/Eigen/src/SVD/SVDBase.h
index e2d77a761..cc90a3b75 100644
--- a/Eigen/src/SVD/SVDBase.h
+++ b/Eigen/src/SVD/SVDBase.h
@@ -130,10 +130,9 @@ public:
inline Index rank() const
{
using std::abs;
- using std::max;
eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
if(m_singularValues.size()==0) return 0;
- RealScalar premultiplied_threshold = (max)(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
+ RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
Index i = m_nonzeroSingularValues-1;
while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
return i+1;
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-19 14:55:35 +0000