aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--test/eigensolver_selfadjoint.cpp16
-rw-r--r--test/svd_common.h61
-rw-r--r--test/svd_fill.h85
3 files changed, 92 insertions, 70 deletions
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 7b0077a6d..4748bdd0b 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -9,9 +9,12 @@
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "main.h"
+#include "svd_fill.h"
#include <limits>
#include <Eigen/Eigenvalues>
+
+
template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
{
typedef typename MatrixType::Index Index;
@@ -31,17 +34,8 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
MatrixType symmC = symmA;
- // randomly nullify some rows/columns
- {
- Index count = 1;//internal::random<Index>(-cols,cols);
- for(Index k=0; k<count; ++k)
- {
- Index i = internal::random<Index>(0,cols-1);
- symmA.row(i).setZero();
- symmA.col(i).setZero();
- }
- }
-
+ svd_fill_random(symmA,Symmetric);
+
symmA.template triangularView<StrictlyUpper>().setZero();
symmC.template triangularView<StrictlyUpper>().setZero();
diff --git a/test/svd_common.h b/test/svd_common.h
index b44b79124..b06a8a0f2 100644
--- a/test/svd_common.h
+++ b/test/svd_common.h
@@ -16,6 +16,8 @@
#error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h
#endif
+#include "svd_fill.h"
+
// Check that the matrix m is properly reconstructed and that the U and V factors are unitary
// The SVD must have already been computed.
template<typename SvdType, typename MatrixType>
@@ -257,65 +259,6 @@ void svd_test_all_computation_options(const MatrixType& m, bool full_only)
}
}
-template<typename MatrixType>
-void svd_fill_random(MatrixType &m)
-{
- typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- Index diagSize = (std::min)(m.rows(), m.cols());
- RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
- s = internal::random<RealScalar>(1,s);
- Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize);
- for(Index k=0; k<diagSize; ++k)
- d(k) = d(k)*std::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
-
- // duplicate some singular values
- if(dup)
- {
- Index n = internal::random<Index>(0,d.size()-1);
- for(Index i=0; i<n; ++i)
- d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1));
- }
-
- Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize);
- Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols());
- if(unit_uv)
- {
- // in very rare cases let's try with a pure diagonal matrix
- if(internal::random<int>(0,10) < 1)
- {
- U.setIdentity();
- VT.setIdentity();
- }
- else
- {
- createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U);
- createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT);
- }
- }
- else
- {
- U.setRandom();
- VT.setRandom();
- }
-
- m = U * d.asDiagonal() * VT;
-
- // (partly) cancel some coeffs
- if(!(dup && unit_uv))
- {
- Matrix<Scalar,Dynamic,1> samples(7);
- samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits<RealScalar>::highest(), 1./NumTraits<RealScalar>::highest();
- Index n = internal::random<Index>(0,m.size()-1);
- for(Index i=0; i<n; ++i)
- m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = samples(internal::random<Index>(0,6));
- }
-}
-
// work around stupid msvc error when constructing at compile time an expression that involves
// a division by zero, even if the numeric type has floating point
diff --git a/test/svd_fill.h b/test/svd_fill.h
new file mode 100644
index 000000000..aa75326d9
--- /dev/null
+++ b/test/svd_fill.h
@@ -0,0 +1,85 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2014-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This 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
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+template<typename MatrixType>
+void svd_fill_random(MatrixType &m, int Option = 0)
+{
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef typename MatrixType::Index Index;
+ Index diagSize = (std::min)(m.rows(), m.cols());
+ RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
+ s = internal::random<RealScalar>(1,s);
+ Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize);
+ for(Index k=0; k<diagSize; ++k)
+ d(k) = d(k)*std::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
+
+ // duplicate some singular values
+ if(dup)
+ {
+ Index n = internal::random<Index>(0,d.size()-1);
+ for(Index i=0; i<n; ++i)
+ d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1));
+ }
+
+ Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize);
+ Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols());
+ if(unit_uv)
+ {
+ // in very rare cases let's try with a pure diagonal matrix
+ if(internal::random<int>(0,10) < 1)
+ {
+ U.setIdentity();
+ VT.setIdentity();
+ }
+ else
+ {
+ createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U);
+ createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT);
+ }
+ }
+ else
+ {
+ U.setRandom();
+ VT.setRandom();
+ }
+
+ if(Option==Symmetric)
+ {
+ m = U * d.asDiagonal() * U.transpose();
+
+ // randomly nullify some rows/columns
+ {
+ Index count = internal::random<Index>(-1,1);
+ for(Index k=0; k<count; ++k)
+ {
+ Index i = internal::random<Index>(0,diagSize-1);
+ m.row(i).setZero();
+ m.col(i).setZero();
+ }
+ }
+ }
+ else
+ {
+ m = U * d.asDiagonal() * VT;
+ // (partly) cancel some coeffs
+ if(!(dup && unit_uv))
+ {
+ Matrix<Scalar,Dynamic,1> samples(7);
+ samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits<RealScalar>::highest(), 1./NumTraits<RealScalar>::highest();
+ Index n = internal::random<Index>(0,m.size()-1);
+ for(Index i=0; i<n; ++i)
+ m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = samples(internal::random<Index>(0,6));
+ }
+ }
+}
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 14:07:02 +0000