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
13template<typename Derived>
14void 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
35template<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
77template<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
110void 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