2
0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-25 18:50:40 +08:00

if EIGEN_NICE_RANDOM is defined, the random functions will return numbers with

few bits left of the comma and for floating-point types will never return zero.
This replaces the custom functions in test/main.h, so one does not anymore need
to think about that when writing tests.
This commit is contained in:
Benoit Jacob 2008-09-01 17:31:21 +00:00
parent 49ff9b204c
commit 46fe7a3d9e
15 changed files with 78 additions and 83 deletions

@ -89,7 +89,14 @@ inline float ei_pow(float x, float y) { return std::pow(x, y); }
template<> inline float ei_random(float a, float b) template<> inline float ei_random(float a, float b)
{ {
#ifdef EIGEN_NICE_RANDOM
int i;
do { i = ei_random<int>(256*int(a),256*int(b));
} while(i==0);
return i/256.f;
#else
return a + (b-a) * std::rand() / RAND_MAX; return a + (b-a) * std::rand() / RAND_MAX;
#endif
} }
template<> inline float ei_random() template<> inline float ei_random()
{ {
@ -123,7 +130,14 @@ inline double ei_pow(double x, double y) { return std::pow(x, y); }
template<> inline double ei_random(double a, double b) template<> inline double ei_random(double a, double b)
{ {
#ifdef EIGEN_NICE_RANDOM
int i;
do { i= ei_random<int>(256*int(a),256*int(b));
} while(i==0);
return i/256.;
#else
return a + (b-a) * std::rand() / RAND_MAX; return a + (b-a) * std::rand() / RAND_MAX;
#endif
} }
template<> inline double ei_random() template<> inline double ei_random()
{ {

@ -41,19 +41,19 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
if (ei_is_same_type<RealScalar,float>::ret) if (ei_is_same_type<RealScalar,float>::ret)
largerEps = 1e-3f; largerEps = 1e-3f;
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols), m3(rows, cols),
mzero = MatrixType::Zero(rows, cols), mzero = MatrixType::Zero(rows, cols),
identity = SquareMatrixType::Identity(rows, rows), identity = SquareMatrixType::Identity(rows, rows),
square = test_random_matrix<SquareMatrixType>(rows, rows); square = SquareMatrixType::Random(rows, rows);
VectorType v1 = test_random_matrix<VectorType>(rows), VectorType v1 = VectorType::Random(rows),
v2 = test_random_matrix<VectorType>(rows), v2 = VectorType::Random(rows),
v3 = test_random_matrix<VectorType>(rows), v3 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
Scalar s1 = test_random<Scalar>(), Scalar s1 = ei_random<Scalar>(),
s2 = test_random<Scalar>(); s2 = ei_random<Scalar>();
// check basic compatibility of adjoint, transpose, conjugate // check basic compatibility of adjoint, transpose, conjugate
VERIFY_IS_APPROX(m1.transpose().conjugate().adjoint(), m1); VERIFY_IS_APPROX(m1.transpose().conjugate().adjoint(), m1);

@ -38,12 +38,12 @@ template<typename MatrixType> void scalarAdd(const MatrixType& m)
int rows = m.rows(); int rows = m.rows();
int cols = m.cols(); int cols = m.cols();
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols); m3(rows, cols);
Scalar s1 = test_random<Scalar>(), Scalar s1 = ei_random<Scalar>(),
s2 = test_random<Scalar>(); s2 = ei_random<Scalar>();
VERIFY_IS_APPROX(m1.cwise() + s1, s1 + m1.cwise()); VERIFY_IS_APPROX(m1.cwise() + s1, s1 + m1.cwise());
VERIFY_IS_APPROX(m1.cwise() + s1, MatrixType::Constant(rows,cols,s1) + m1); VERIFY_IS_APPROX(m1.cwise() + s1, MatrixType::Constant(rows,cols,s1) + m1);

@ -34,18 +34,18 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
// this test relies a lot on Random.h, and there's not much more that we can do // this test relies a lot on Random.h, and there's not much more that we can do
// to test it, hence I consider that we will have tested Random.h // to test it, hence I consider that we will have tested Random.h
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols), m3(rows, cols),
mzero = MatrixType::Zero(rows, cols), mzero = MatrixType::Zero(rows, cols),
identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
::Identity(rows, rows), ::Identity(rows, rows),
square = test_random_matrix<Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> >(rows, rows); square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
VectorType v1 = test_random_matrix<VectorType>(rows), VectorType v1 = VectorType::Random(rows),
v2 = test_random_matrix<VectorType>(rows), v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
Scalar x = test_random<Scalar>(); Scalar x = ei_random<Scalar>();
int r = ei_random<int>(0, rows-1), int r = ei_random<int>(0, rows-1),
c = ei_random<int>(0, cols-1); c = ei_random<int>(0, cols-1);

@ -43,12 +43,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
MatrixType a0 = test_random_matrix<MatrixType>(rows,cols); MatrixType a0 = MatrixType::Random(rows,cols);
VectorType vecB = test_random_matrix<VectorType>(rows); VectorType vecB = VectorType::Random(rows);
MatrixType matB = test_random_matrix<MatrixType>(rows,cols); MatrixType matB = MatrixType::Random(rows,cols);
SquareMatrixType symm = a0 * a0.adjoint(); SquareMatrixType symm = a0 * a0.adjoint();
// let's make sure the matrix is not singular or near singular // let's make sure the matrix is not singular or near singular
MatrixType a1 = test_random_matrix<MatrixType>(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols);
symm += a1 * a1.adjoint(); symm += a1 * a1.adjoint();
#ifdef HAS_GSL #ifdef HAS_GSL

@ -43,16 +43,16 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
int rows = m.rows(); int rows = m.rows();
int cols = m.cols(); int cols = m.cols();
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols), m3(rows, cols),
mzero = MatrixType::Zero(rows, cols), mzero = MatrixType::Zero(rows, cols),
mones = MatrixType::Ones(rows, cols), mones = MatrixType::Ones(rows, cols),
identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
::Identity(rows, rows), ::Identity(rows, rows),
square = test_random_matrix<Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> >(rows, rows); square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
VectorType v1 = test_random_matrix<VectorType>(rows), VectorType v1 = VectorType::Random(rows),
v2 = test_random_matrix<VectorType>(rows), v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
int r = ei_random<int>(0, rows-1), int r = ei_random<int>(0, rows-1),

@ -45,12 +45,12 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
RealScalar largerEps = 10*test_precision<RealScalar>(); RealScalar largerEps = 10*test_precision<RealScalar>();
MatrixType a = test_random_matrix<MatrixType>(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
MatrixType a1 = test_random_matrix<MatrixType>(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols);
MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
MatrixType b = test_random_matrix<MatrixType>(rows,cols); MatrixType b = MatrixType::Random(rows,cols);
MatrixType b1 = test_random_matrix<MatrixType>(rows,cols); MatrixType b1 = MatrixType::Random(rows,cols);
MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1; MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
SelfAdjointEigenSolver<MatrixType> eiSymm(symmA); SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);

@ -49,10 +49,10 @@ template<typename Scalar> void geometry(void)
typedef Translation<Scalar,3> Translation3; typedef Translation<Scalar,3> Translation3;
Quaternion q1, q2; Quaternion q1, q2;
Vector3 v0 = test_random_matrix<Vector3>(), Vector3 v0 = Vector3::Random(),
v1 = test_random_matrix<Vector3>(), v1 = Vector3::Random(),
v2 = test_random_matrix<Vector3>(); v2 = Vector3::Random();
Vector2 u0 = test_random_matrix<Vector2>(); Vector2 u0 = Vector2::Random();
Matrix3 matrot1; Matrix3 matrot1;
Scalar a = ei_random<Scalar>(-M_PI, M_PI); Scalar a = ei_random<Scalar>(-M_PI, M_PI);
@ -134,7 +134,7 @@ template<typename Scalar> void geometry(void)
t1.setIdentity(); t1.setIdentity();
t1.linear() = q1.toRotationMatrix(); t1.linear() = q1.toRotationMatrix();
v0 << 50, 2, 1;//= test_random_matrix<Vector3>().cwiseProduct(Vector3(10,2,0.5)); v0 << 50, 2, 1;//= ei_random_matrix<Vector3>().cwiseProduct(Vector3(10,2,0.5));
t0.scale(v0); t0.scale(v0);
t1.prescale(v0); t1.prescale(v0);
@ -169,8 +169,8 @@ template<typename Scalar> void geometry(void)
// 2D transformation // 2D transformation
Transform2 t20, t21; Transform2 t20, t21;
Vector2 v20 = test_random_matrix<Vector2>(); Vector2 v20 = Vector2::Random();
Vector2 v21 = test_random_matrix<Vector2>(); Vector2 v21 = Vector2::Random();
for (int k=0; k<2; ++k) for (int k=0; k<2; ++k)
if (ei_abs(v21[k])<1e-3) v21[k] = 1e-3; if (ei_abs(v21[k])<1e-3) v21[k] = 1e-3;
t21.setIdentity(); t21.setIdentity();

@ -38,7 +38,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2(rows, cols), m2(rows, cols),
mzero = MatrixType::Zero(rows, cols), mzero = MatrixType::Zero(rows, cols),
identity = MatrixType::Identity(rows, rows); identity = MatrixType::Identity(rows, rows);
@ -46,7 +46,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
if (ei_is_same_type<RealScalar,float>::ret) if (ei_is_same_type<RealScalar,float>::ret)
{ {
// let's build a more stable to inverse matrix // let's build a more stable to inverse matrix
MatrixType a = test_random_matrix<MatrixType>(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
m1 += m1 * m1.adjoint() + a * a.adjoint(); m1 += m1 * m1.adjoint() + a * a.adjoint();
} }

@ -38,13 +38,13 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
// this test relies a lot on Random.h, and there's not much more that we can do // this test relies a lot on Random.h, and there's not much more that we can do
// to test it, hence I consider that we will have tested Random.h // to test it, hence I consider that we will have tested Random.h
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols), m3(rows, cols),
mzero = MatrixType::Zero(rows, cols); mzero = MatrixType::Zero(rows, cols);
Scalar s1 = test_random<Scalar>(); Scalar s1 = ei_random<Scalar>();
while (ei_abs(s1)<1e-3) s1 = test_random<Scalar>(); while (ei_abs(s1)<1e-3) s1 = ei_random<Scalar>();
int r = ei_random<int>(0, rows-1), int r = ei_random<int>(0, rows-1),
c = ei_random<int>(0, cols-1); c = ei_random<int>(0, cols-1);

@ -56,7 +56,7 @@ template<typename MatrixType> void lu_non_invertible()
int rank = ei_random<int>(1, std::min(rows, cols)-1); int rank = ei_random<int>(1, std::min(rows, cols)-1);
MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
m1 = test_random_matrix<MatrixType>(rows,cols); m1 = MatrixType::Random(rows,cols);
if(rows <= cols) if(rows <= cols)
for(int i = rank; i < rows; i++) m1.row(i).setZero(); for(int i = rank; i < rows; i++) m1.row(i).setZero();
else else
@ -72,12 +72,12 @@ template<typename MatrixType> void lu_non_invertible()
VERIFY((m1 * lu.kernel()).isMuchSmallerThan(m1)); VERIFY((m1 * lu.kernel()).isMuchSmallerThan(m1));
lu.computeKernel(&k); lu.computeKernel(&k);
VERIFY((m1 * k).isMuchSmallerThan(m1)); VERIFY((m1 * k).isMuchSmallerThan(m1));
m2 = test_random_matrix<MatrixType>(cols,cols2); m2 = MatrixType::Random(cols,cols2);
m3 = m1*m2; m3 = m1*m2;
m2 = test_random_matrix<MatrixType>(cols,cols2); m2 = MatrixType::Random(cols,cols2);
lu.solve(m3, &m2); lu.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2); VERIFY_IS_APPROX(m3, m1*m2);
m3 = test_random_matrix<MatrixType>(rows,cols2); m3 = MatrixType::Random(rows,cols2);
VERIFY(!lu.solve(m3, &m2)); VERIFY(!lu.solve(m3, &m2));
} }
@ -90,12 +90,12 @@ template<typename MatrixType> void lu_invertible()
int size = ei_random<int>(10,200); int size = ei_random<int>(10,200);
MatrixType m1(size, size), m2(size, size), m3(size, size); MatrixType m1(size, size), m2(size, size), m3(size, size);
m1 = test_random_matrix<MatrixType>(size,size); m1 = MatrixType::Random(size,size);
if (ei_is_same_type<RealScalar,float>::ret) if (ei_is_same_type<RealScalar,float>::ret)
{ {
// let's build a matrix more stable to inverse // let's build a matrix more stable to inverse
MatrixType a = test_random_matrix<MatrixType>(size,size*2); MatrixType a = MatrixType::Random(size,size*2);
m1 += a * a.adjoint(); m1 += a * a.adjoint();
} }
@ -105,11 +105,11 @@ template<typename MatrixType> void lu_invertible()
VERIFY(lu.isInjective()); VERIFY(lu.isInjective());
VERIFY(lu.isSurjective()); VERIFY(lu.isSurjective());
VERIFY(lu.isInvertible()); VERIFY(lu.isInvertible());
m3 = test_random_matrix<MatrixType>(size,size); m3 = MatrixType::Random(size,size);
lu.solve(m3, &m2); lu.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2); VERIFY_IS_APPROX(m3, m1*m2);
VERIFY_IS_APPROX(m2, lu.inverse()*m3); VERIFY_IS_APPROX(m2, lu.inverse()*m3);
m3 = test_random_matrix<MatrixType>(size,size); m3 = MatrixType::Random(size,size);
VERIFY(lu.solve(m3, &m2)); VERIFY(lu.solve(m3, &m2));
} }

@ -135,6 +135,7 @@ namespace Eigen
#define EIGEN_INTERNAL_DEBUGGING #define EIGEN_INTERNAL_DEBUGGING
#define EIGEN_NICE_RANDOM
#include <Eigen/Core> #include <Eigen/Core>
namespace Eigen { namespace Eigen {
@ -221,26 +222,6 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m,
return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>()); return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>());
} }
template<typename T> T test_random();
template<> int test_random() { return ei_random<int>(-100,100); }
template<> float test_random() { return float(ei_random<int>(-1000,1000)) / 256.f; }
template<> double test_random() { return double(ei_random<int>(-1000,1000)) / 256.; }
template<> std::complex<float> test_random()
{ return std::complex<float>(test_random<float>(),test_random<float>()); }
template<> std::complex<double> test_random()
{ return std::complex<double>(test_random<double>(),test_random<double>()); }
template<typename MatrixType>
MatrixType test_random_matrix(int rows = MatrixType::RowsAtCompileTime, int cols = MatrixType::ColsAtCompileTime)
{
MatrixType res(rows, cols);
for (int j=0; j<cols; ++j)
for (int i=0; i<rows; ++i)
res.coeffRef(i,j) = test_random<typename MatrixType::Scalar>();
return res;
}
} // end namespace Eigen } // end namespace Eigen

@ -31,7 +31,7 @@ template<typename MatrixType> void matrixSum(const MatrixType& m)
int rows = m.rows(); int rows = m.rows();
int cols = m.cols(); int cols = m.cols();
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols); MatrixType m1 = MatrixType::Random(rows, cols);
VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows, cols).sum(), Scalar(1)); VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows, cols).sum(), Scalar(1));
VERIFY_IS_APPROX(MatrixType::Ones(rows, cols).sum(), Scalar(rows*cols)); VERIFY_IS_APPROX(MatrixType::Ones(rows, cols).sum(), Scalar(rows*cols));
@ -45,7 +45,7 @@ template<typename VectorType> void vectorSum(const VectorType& w)
typedef typename VectorType::Scalar Scalar; typedef typename VectorType::Scalar Scalar;
int size = w.size(); int size = w.size();
VectorType v = test_random_matrix<VectorType>(size); VectorType v = VectorType::Random(size);
for(int i = 1; i < size; i++) for(int i = 1; i < size; i++)
{ {
Scalar s = Scalar(0); Scalar s = Scalar(0);

@ -35,9 +35,9 @@ template<typename MatrixType> void svd(const MatrixType& m)
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
MatrixType a = test_random_matrix<MatrixType>(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b =
test_random_matrix<Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> >(rows,1); Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1);
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1);
RealScalar largerEps = test_precision<RealScalar>(); RealScalar largerEps = test_precision<RealScalar>();
@ -56,7 +56,7 @@ template<typename MatrixType> void svd(const MatrixType& m)
{ {
if (ei_is_same_type<RealScalar,float>::ret) if (ei_is_same_type<RealScalar,float>::ret)
{ {
MatrixType a1 = test_random_matrix<MatrixType>(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols);
a += a * a.adjoint() + a1 * a1.adjoint(); a += a * a.adjoint() + a1 * a1.adjoint();
} }
SVD<MatrixType> svd(a); SVD<MatrixType> svd(a);

@ -35,8 +35,8 @@ template<typename MatrixType> void triangular(const MatrixType& m)
int rows = m.rows(); int rows = m.rows();
int cols = m.cols(); int cols = m.cols();
MatrixType m1 = test_random_matrix<MatrixType>(rows, cols), MatrixType m1 = MatrixType::Random(rows, cols),
m2 = test_random_matrix<MatrixType>(rows, cols), m2 = MatrixType::Random(rows, cols),
m3(rows, cols), m3(rows, cols),
m4(rows, cols), m4(rows, cols),
r1(rows, cols), r1(rows, cols),
@ -47,8 +47,8 @@ template<typename MatrixType> void triangular(const MatrixType& m)
::Identity(rows, rows), ::Identity(rows, rows),
square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
::Random(rows, rows); ::Random(rows, rows);
VectorType v1 = test_random_matrix<VectorType>(rows), VectorType v1 = VectorType::Random(rows),
v2 = test_random_matrix<VectorType>(rows), v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
MatrixType m1up = m1.template part<Eigen::Upper>(); MatrixType m1up = m1.template part<Eigen::Upper>();
@ -81,9 +81,9 @@ template<typename MatrixType> void triangular(const MatrixType& m)
m1.template part<Eigen::Lower>() = (m2.transpose() * m2).lazy(); m1.template part<Eigen::Lower>() = (m2.transpose() * m2).lazy();
VERIFY_IS_APPROX(m3.template part<Eigen::Lower>(), m1); VERIFY_IS_APPROX(m3.template part<Eigen::Lower>(), m1);
m1 = test_random_matrix<MatrixType>(rows, cols); m1 = MatrixType::Random(rows, cols);
for (int i=0; i<rows; ++i) for (int i=0; i<rows; ++i)
while (ei_abs2(m1(i,i))<1e-3) m1(i,i) = test_random<Scalar>(); while (ei_abs2(m1(i,i))<1e-3) m1(i,i) = ei_random<Scalar>();
Transpose<MatrixType> trm4(m4); Transpose<MatrixType> trm4(m4);
// test back and forward subsitution // test back and forward subsitution