xref: /aosp_15_r20/external/eigen/test/inverse.cpp (revision bf2c37156dfe67e5dfebd6d394bad8b2ab5804d4)
1*bf2c3715SXin Li // This file is part of Eigen, a lightweight C++ template library
2*bf2c3715SXin Li // for linear algebra.
3*bf2c3715SXin Li //
4*bf2c3715SXin Li // Copyright (C) 2008 Gael Guennebaud <[email protected]>
5*bf2c3715SXin Li // Copyright (C) 2008 Benoit Jacob <[email protected]>
6*bf2c3715SXin Li //
7*bf2c3715SXin Li // This Source Code Form is subject to the terms of the Mozilla
8*bf2c3715SXin Li // Public License v. 2.0. If a copy of the MPL was not distributed
9*bf2c3715SXin Li // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10*bf2c3715SXin Li 
11*bf2c3715SXin Li #include "main.h"
12*bf2c3715SXin Li #include <Eigen/LU>
13*bf2c3715SXin Li 
14*bf2c3715SXin Li template<typename MatrixType>
inverse_for_fixed_size(const MatrixType &,typename internal::enable_if<MatrixType::SizeAtCompileTime==Dynamic>::type * =0)15*bf2c3715SXin Li void inverse_for_fixed_size(const MatrixType&, typename internal::enable_if<MatrixType::SizeAtCompileTime==Dynamic>::type* = 0)
16*bf2c3715SXin Li {
17*bf2c3715SXin Li }
18*bf2c3715SXin Li 
19*bf2c3715SXin Li template<typename MatrixType>
inverse_for_fixed_size(const MatrixType & m1,typename internal::enable_if<MatrixType::SizeAtCompileTime!=Dynamic>::type * =0)20*bf2c3715SXin Li void inverse_for_fixed_size(const MatrixType& m1, typename internal::enable_if<MatrixType::SizeAtCompileTime!=Dynamic>::type* = 0)
21*bf2c3715SXin Li {
22*bf2c3715SXin Li   using std::abs;
23*bf2c3715SXin Li 
24*bf2c3715SXin Li   MatrixType m2, identity = MatrixType::Identity();
25*bf2c3715SXin Li 
26*bf2c3715SXin Li   typedef typename MatrixType::Scalar Scalar;
27*bf2c3715SXin Li   typedef typename NumTraits<Scalar>::Real RealScalar;
28*bf2c3715SXin Li   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
29*bf2c3715SXin Li 
30*bf2c3715SXin Li   //computeInverseAndDetWithCheck tests
31*bf2c3715SXin Li   //First: an invertible matrix
32*bf2c3715SXin Li   bool invertible;
33*bf2c3715SXin Li   Scalar det;
34*bf2c3715SXin Li 
35*bf2c3715SXin Li   m2.setZero();
36*bf2c3715SXin Li   m1.computeInverseAndDetWithCheck(m2, det, invertible);
37*bf2c3715SXin Li   VERIFY(invertible);
38*bf2c3715SXin Li   VERIFY_IS_APPROX(identity, m1*m2);
39*bf2c3715SXin Li   VERIFY_IS_APPROX(det, m1.determinant());
40*bf2c3715SXin Li 
41*bf2c3715SXin Li   m2.setZero();
42*bf2c3715SXin Li   m1.computeInverseWithCheck(m2, invertible);
43*bf2c3715SXin Li   VERIFY(invertible);
44*bf2c3715SXin Li   VERIFY_IS_APPROX(identity, m1*m2);
45*bf2c3715SXin Li 
46*bf2c3715SXin Li   //Second: a rank one matrix (not invertible, except for 1x1 matrices)
47*bf2c3715SXin Li   VectorType v3 = VectorType::Random();
48*bf2c3715SXin Li   MatrixType m3 = v3*v3.transpose(), m4;
49*bf2c3715SXin Li   m3.computeInverseAndDetWithCheck(m4, det, invertible);
50*bf2c3715SXin Li   VERIFY( m1.rows()==1 ? invertible : !invertible );
51*bf2c3715SXin Li   VERIFY_IS_MUCH_SMALLER_THAN(abs(det-m3.determinant()), RealScalar(1));
52*bf2c3715SXin Li   m3.computeInverseWithCheck(m4, invertible);
53*bf2c3715SXin Li   VERIFY( m1.rows()==1 ? invertible : !invertible );
54*bf2c3715SXin Li 
55*bf2c3715SXin Li   // check with submatrices
56*bf2c3715SXin Li   {
57*bf2c3715SXin Li     Matrix<Scalar, MatrixType::RowsAtCompileTime+1, MatrixType::RowsAtCompileTime+1, MatrixType::Options> m5;
58*bf2c3715SXin Li     m5.setRandom();
59*bf2c3715SXin Li     m5.topLeftCorner(m1.rows(),m1.rows()) = m1;
60*bf2c3715SXin Li     m2 = m5.template topLeftCorner<MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime>().inverse();
61*bf2c3715SXin Li     VERIFY_IS_APPROX( (m5.template topLeftCorner<MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime>()), m2.inverse() );
62*bf2c3715SXin Li   }
63*bf2c3715SXin Li }
64*bf2c3715SXin Li 
inverse(const MatrixType & m)65*bf2c3715SXin Li template<typename MatrixType> void inverse(const MatrixType& m)
66*bf2c3715SXin Li {
67*bf2c3715SXin Li   /* this test covers the following files:
68*bf2c3715SXin Li      Inverse.h
69*bf2c3715SXin Li   */
70*bf2c3715SXin Li   Index rows = m.rows();
71*bf2c3715SXin Li   Index cols = m.cols();
72*bf2c3715SXin Li 
73*bf2c3715SXin Li   typedef typename MatrixType::Scalar Scalar;
74*bf2c3715SXin Li 
75*bf2c3715SXin Li   MatrixType m1(rows, cols),
76*bf2c3715SXin Li              m2(rows, cols),
77*bf2c3715SXin Li              identity = MatrixType::Identity(rows, rows);
78*bf2c3715SXin Li   createRandomPIMatrixOfRank(rows,rows,rows,m1);
79*bf2c3715SXin Li   m2 = m1.inverse();
80*bf2c3715SXin Li   VERIFY_IS_APPROX(m1, m2.inverse() );
81*bf2c3715SXin Li 
82*bf2c3715SXin Li   VERIFY_IS_APPROX((Scalar(2)*m2).inverse(), m2.inverse()*Scalar(0.5));
83*bf2c3715SXin Li 
84*bf2c3715SXin Li   VERIFY_IS_APPROX(identity, m1.inverse() * m1 );
85*bf2c3715SXin Li   VERIFY_IS_APPROX(identity, m1 * m1.inverse() );
86*bf2c3715SXin Li 
87*bf2c3715SXin Li   VERIFY_IS_APPROX(m1, m1.inverse().inverse() );
88*bf2c3715SXin Li 
89*bf2c3715SXin Li   // since for the general case we implement separately row-major and col-major, test that
90*bf2c3715SXin Li   VERIFY_IS_APPROX(MatrixType(m1.transpose().inverse()), MatrixType(m1.inverse().transpose()));
91*bf2c3715SXin Li 
92*bf2c3715SXin Li   inverse_for_fixed_size(m1);
93*bf2c3715SXin Li 
94*bf2c3715SXin Li   // check in-place inversion
95*bf2c3715SXin Li   if(MatrixType::RowsAtCompileTime>=2 && MatrixType::RowsAtCompileTime<=4)
96*bf2c3715SXin Li   {
97*bf2c3715SXin Li     // in-place is forbidden
98*bf2c3715SXin Li     VERIFY_RAISES_ASSERT(m1 = m1.inverse());
99*bf2c3715SXin Li   }
100*bf2c3715SXin Li   else
101*bf2c3715SXin Li   {
102*bf2c3715SXin Li     m2 = m1.inverse();
103*bf2c3715SXin Li     m1 = m1.inverse();
104*bf2c3715SXin Li     VERIFY_IS_APPROX(m1,m2);
105*bf2c3715SXin Li   }
106*bf2c3715SXin Li }
107*bf2c3715SXin Li 
108*bf2c3715SXin Li template<typename Scalar>
inverse_zerosized()109*bf2c3715SXin Li void inverse_zerosized()
110*bf2c3715SXin Li {
111*bf2c3715SXin Li   Matrix<Scalar,Dynamic,Dynamic> A(0,0);
112*bf2c3715SXin Li   {
113*bf2c3715SXin Li     Matrix<Scalar,0,1> b, x;
114*bf2c3715SXin Li     x = A.inverse() * b;
115*bf2c3715SXin Li   }
116*bf2c3715SXin Li   {
117*bf2c3715SXin Li     Matrix<Scalar,Dynamic,Dynamic> b(0,1), x;
118*bf2c3715SXin Li     x = A.inverse() * b;
119*bf2c3715SXin Li     VERIFY_IS_EQUAL(x.rows(), 0);
120*bf2c3715SXin Li     VERIFY_IS_EQUAL(x.cols(), 1);
121*bf2c3715SXin Li   }
122*bf2c3715SXin Li }
123*bf2c3715SXin Li 
EIGEN_DECLARE_TEST(inverse)124*bf2c3715SXin Li EIGEN_DECLARE_TEST(inverse)
125*bf2c3715SXin Li {
126*bf2c3715SXin Li   int s = 0;
127*bf2c3715SXin Li   for(int i = 0; i < g_repeat; i++) {
128*bf2c3715SXin Li     CALL_SUBTEST_1( inverse(Matrix<double,1,1>()) );
129*bf2c3715SXin Li     CALL_SUBTEST_2( inverse(Matrix2d()) );
130*bf2c3715SXin Li     CALL_SUBTEST_3( inverse(Matrix3f()) );
131*bf2c3715SXin Li     CALL_SUBTEST_4( inverse(Matrix4f()) );
132*bf2c3715SXin Li     CALL_SUBTEST_4( inverse(Matrix<float,4,4,DontAlign>()) );
133*bf2c3715SXin Li 
134*bf2c3715SXin Li     s = internal::random<int>(50,320);
135*bf2c3715SXin Li     CALL_SUBTEST_5( inverse(MatrixXf(s,s)) );
136*bf2c3715SXin Li     TEST_SET_BUT_UNUSED_VARIABLE(s)
137*bf2c3715SXin Li     CALL_SUBTEST_5( inverse_zerosized<float>() );
138*bf2c3715SXin Li     CALL_SUBTEST_5( inverse(MatrixXf(0, 0)) );
139*bf2c3715SXin Li     CALL_SUBTEST_5( inverse(MatrixXf(1, 1)) );
140*bf2c3715SXin Li 
141*bf2c3715SXin Li     s = internal::random<int>(25,100);
142*bf2c3715SXin Li     CALL_SUBTEST_6( inverse(MatrixXcd(s,s)) );
143*bf2c3715SXin Li     TEST_SET_BUT_UNUSED_VARIABLE(s)
144*bf2c3715SXin Li 
145*bf2c3715SXin Li     CALL_SUBTEST_7( inverse(Matrix4d()) );
146*bf2c3715SXin Li     CALL_SUBTEST_7( inverse(Matrix<double,4,4,DontAlign>()) );
147*bf2c3715SXin Li 
148*bf2c3715SXin Li     CALL_SUBTEST_8( inverse(Matrix4cd()) );
149*bf2c3715SXin Li   }
150*bf2c3715SXin Li }
151