1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 #include "main.h"
11 #include <Eigen/LU>
12
13 template<typename Derived>
doSomeRankPreservingOperations(Eigen::MatrixBase<Derived> & m)14 void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
15 {
16 typedef typename Derived::RealScalar RealScalar;
17 for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
18 {
19 RealScalar d = Eigen::ei_random<RealScalar>(-1,1);
20 int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
21 int j;
22 do {
23 j = Eigen::ei_random<int>(0,m.rows()-1);
24 } while (i==j); // j is another one (must be different)
25 m.row(i) += d * m.row(j);
26
27 i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number
28 do {
29 j = Eigen::ei_random<int>(0,m.cols()-1);
30 } while (i==j); // j is another one (must be different)
31 m.col(i) += d * m.col(j);
32 }
33 }
34
lu_non_invertible()35 template<typename MatrixType> void lu_non_invertible()
36 {
37 /* this test covers the following files:
38 LU.h
39 */
40 // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function
41 int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
42 int rank = ei_random<int>(1, std::min(rows, cols)-1);
43
44 MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
45 m1 = MatrixType::Random(rows,cols);
46 if(rows <= cols)
47 for(int i = rank; i < rows; i++) m1.row(i).setZero();
48 else
49 for(int i = rank; i < cols; i++) m1.col(i).setZero();
50 doSomeRankPreservingOperations(m1);
51
52 LU<MatrixType> lu(m1);
53 typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel();
54 typename LU<MatrixType>::ImageResultType m1image = lu.image();
55
56 VERIFY(rank == lu.rank());
57 VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
58 VERIFY(!lu.isInjective());
59 VERIFY(!lu.isInvertible());
60 VERIFY(lu.isSurjective() == (lu.rank() == rows));
61 VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
62 VERIFY(m1image.lu().rank() == rank);
63 MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
64 sidebyside << m1, m1image;
65 VERIFY(sidebyside.lu().rank() == rank);
66 m2 = MatrixType::Random(cols,cols2);
67 m3 = m1*m2;
68 m2 = MatrixType::Random(cols,cols2);
69 lu.solve(m3, &m2);
70 VERIFY_IS_APPROX(m3, m1*m2);
71 /* solve now always returns true
72 m3 = MatrixType::Random(rows,cols2);
73 VERIFY(!lu.solve(m3, &m2));
74 */
75 }
76
lu_invertible()77 template<typename MatrixType> void lu_invertible()
78 {
79 /* this test covers the following files:
80 LU.h
81 */
82 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
83 int size = ei_random<int>(10,200);
84
85 MatrixType m1(size, size), m2(size, size), m3(size, size);
86 m1 = MatrixType::Random(size,size);
87
88 if (ei_is_same_type<RealScalar,float>::ret)
89 {
90 // let's build a matrix more stable to inverse
91 MatrixType a = MatrixType::Random(size,size*2);
92 m1 += a * a.adjoint();
93 }
94
95 LU<MatrixType> lu(m1);
96 VERIFY(0 == lu.dimensionOfKernel());
97 VERIFY(size == lu.rank());
98 VERIFY(lu.isInjective());
99 VERIFY(lu.isSurjective());
100 VERIFY(lu.isInvertible());
101 VERIFY(lu.image().lu().isInvertible());
102 m3 = MatrixType::Random(size,size);
103 lu.solve(m3, &m2);
104 VERIFY_IS_APPROX(m3, m1*m2);
105 VERIFY_IS_APPROX(m2, lu.inverse()*m3);
106 m3 = MatrixType::Random(size,size);
107 VERIFY(lu.solve(m3, &m2));
108 }
109
test_eigen2_lu()110 void test_eigen2_lu()
111 {
112 for(int i = 0; i < g_repeat; i++) {
113 CALL_SUBTEST_1( lu_non_invertible<MatrixXf>() );
114 CALL_SUBTEST_2( lu_non_invertible<MatrixXd>() );
115 CALL_SUBTEST_3( lu_non_invertible<MatrixXcf>() );
116 CALL_SUBTEST_4( lu_non_invertible<MatrixXcd>() );
117 CALL_SUBTEST_1( lu_invertible<MatrixXf>() );
118 CALL_SUBTEST_2( lu_invertible<MatrixXd>() );
119 CALL_SUBTEST_3( lu_invertible<MatrixXcf>() );
120 CALL_SUBTEST_4( lu_invertible<MatrixXcd>() );
121 }
122 }
123