devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix imag_ref for real scalar types and clean up svd_fill.h

Browse Source 
libeigen/eigen!2303

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-03-15 19:56:01 -07:00
parent 929785924c
commit ea13a98dec
2 changed files with 45 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Core/MathFunctions.h
View File

@@ -170,18 +170,24 @@ struct imag_ref_default_impl {
template <typename Scalar>
struct imag_ref_default_impl<Scalar, false> {
EIGEN_DEVICE_FUNC constexpr static Scalar run(Scalar&) { return Scalar(0); }
EIGEN_DEVICE_FUNC constexpr static const Scalar run(const Scalar&) { return Scalar(0); }
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(Scalar&) { return RealScalar(0); }
EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(const Scalar&) { return RealScalar(0); }
};
template <typename Scalar>
struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
template <typename Scalar>
template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct imag_ref_retval {
typedef typename NumTraits<Scalar>::Real& type;
};
template <typename Scalar>
struct imag_ref_retval<Scalar, false> {
typedef typename NumTraits<Scalar>::Real type;
};
} // namespace internal
namespace numext {

47
test/svd_fill.h
View File

@@ -23,6 +23,15 @@ Array<T, 4, 1> four_denorms() {
return four_denorms<double>().cast<T>();
}
template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct maybe_set_imag_part {
static void run(Scalar& x, typename NumTraits<Scalar>::Real y) { numext::imag_ref(x) = y; }
};
template <typename Scalar>
struct maybe_set_imag_part<Scalar, false> {
static void run(Scalar&, typename NumTraits<Scalar>::Real) {}
};
template <typename MatrixType>
void svd_fill_random(MatrixType &m, int Option = 0) {
using std::pow;
@@ -35,9 +44,12 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
for (Index k = 0; k < diagSize; ++k) d(k) = d(k) * pow(RealScalar(10), internal::random<RealScalar>(-s, s));
bool dup = internal::random<int>(0, 10) < 3;
bool unit_uv =
internal::random<int>(0, 10) < (dup ? 7 : 3); // if we duplicate some diagonal entries, then increase the chance
// to preserve them using unitary U and V factors
// For SelfAdjoint, always use unitary U so the eigenvalues are exactly d
// and the conditioning is controlled. For Symmetric/general SVD, randomly
// choose unitary or arbitrary U/V.
bool unit_uv = (Option == SelfAdjoint) ||
internal::random<int>(0, 10) < (dup ? 7 : 3); // if we duplicate some diagonal entries, then increase
// the chance to preserve them using unitary U and V
// duplicate some singular values
if (dup) {
@@ -67,8 +79,11 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
RealScalar(1) / NumTraits<RealScalar>::highest(), (std::numeric_limits<RealScalar>::min)(),
pow((std::numeric_limits<RealScalar>::min)(), RealScalar(0.8));
if (Option == Symmetric) {
m = U * d.asDiagonal() * U.transpose();
if (Option == Symmetric || Option == SelfAdjoint) {
if (Option == SelfAdjoint)
m = U * d.asDiagonal() * U.adjoint();
else
m = U * d.asDiagonal() * U.transpose();
// randomly nullify some rows/columns
{
@@ -85,10 +100,21 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
for (Index k = 0; k < n; ++k) {
Index i = internal::random<Index>(0, m.rows() - 1);
Index j = internal::random<Index>(0, m.cols() - 1);
m(j, i) = m(i, j) = samples(internal::random<Index>(0, samples.size() - 1));
if (NumTraits<Scalar>::IsComplex)
*(&numext::real_ref(m(j, i)) + 1) = *(&numext::real_ref(m(i, j)) + 1) =
samples.real()(internal::random<Index>(0, samples.size() - 1));
Scalar val = samples(internal::random<Index>(0, samples.size() - 1));
maybe_set_imag_part<Scalar>::run(val, samples.real()(internal::random<Index>(0, samples.size() - 1)));
if (Option == SelfAdjoint) {
if (i == j) {
// Diagonal entries of a Hermitian matrix must be real.
m(i, i) = numext::real(val);
} else {
// Off-diagonal: m(j,i) = conj(m(i,j)) for Hermitian.
m(i, j) = val;
m(j, i) = numext::conj(val);
}
} else {
m(i, j) = val;
m(j, i) = val;
}
}
}
}
@@ -101,8 +127,7 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
Index i = internal::random<Index>(0, m.rows() - 1);
Index j = internal::random<Index>(0, m.cols() - 1);
m(i, j) = samples(internal::random<Index>(0, samples.size() - 1));
if (NumTraits<Scalar>::IsComplex)
*(&numext::real_ref(m(i, j)) + 1) = samples.real()(internal::random<Index>(0, samples.size() - 1));
maybe_set_imag_part<Scalar>::run(m(i, j), samples.real()(internal::random<Index>(0, samples.size() - 1)));
}
}
}