ground/gcs/src/libs/eigen/test/triangular.cpp

   1 // This file is triangularView of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
   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
  12
  13
  14 template<typename MatrixType> void triangular_square(const MatrixType& m)
  15 {
  16   typedef typename MatrixType::Scalar Scalar;
  17   typedef typename NumTraits<Scalar>::Real RealScalar;
  18   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  19
  20   RealScalar largerEps = 10*test_precision<RealScalar>();
  21
  22   typename MatrixType::Index rows = m.rows();
  23   typename MatrixType::Index cols = m.cols();
  24
  25   MatrixType m1 = MatrixType::Random(rows, cols),
  26              m2 = MatrixType::Random(rows, cols),
  27              m3(rows, cols),
  28              m4(rows, cols),
  29              r1(rows, cols),
  30              r2(rows, cols);
  31   VectorType v2 = VectorType::Random(rows);
  32
  33   MatrixType m1up = m1.template triangularView<Upper>();
  34   MatrixType m2up = m2.template triangularView<Upper>();
  35
  36   if (rows*cols>1)
  37   {
  38     VERIFY(m1up.isUpperTriangular());
  39     VERIFY(m2up.transpose().isLowerTriangular());
  40     VERIFY(!m2.isLowerTriangular());
  41   }
  42
  43 //   VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2);
  44
  45   // test overloaded operator+=
  46   r1.setZero();
  47   r2.setZero();
  48   r1.template triangularView<Upper>() +=  m1;
  49   r2 += m1up;
  50   VERIFY_IS_APPROX(r1,r2);
  51
  52   // test overloaded operator=
  53   m1.setZero();
  54   m1.template triangularView<Upper>() = m2.transpose() + m2;
  55   m3 = m2.transpose() + m2;
  56   VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1);
  57
  58   // test overloaded operator=
  59   m1.setZero();
  60   m1.template triangularView<Lower>() = m2.transpose() + m2;
  61   VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
  62
  63   VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(),
  64                    m3.conjugate().template triangularView<Lower>().toDenseMatrix());
  65
  66   m1 = MatrixType::Random(rows, cols);
  67   for (int i=0; i<rows; ++i)
  68     while (numext::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random<Scalar>();
  69
  70   Transpose<MatrixType> trm4(m4);
  71   // test back and forward subsitution with a vector as the rhs
  72   m3 = m1.template triangularView<Upper>();
  73   VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps));
  74   m3 = m1.template triangularView<Lower>();
  75   VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps));
  76   m3 = m1.template triangularView<Upper>();
  77   VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps));
  78   m3 = m1.template triangularView<Lower>();
  79   VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps));
  80
  81   // test back and forward subsitution with a matrix as the rhs
  82   m3 = m1.template triangularView<Upper>();
  83   VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps));
  84   m3 = m1.template triangularView<Lower>();
  85   VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps));
  86   m3 = m1.template triangularView<Upper>();
  87   VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps));
  88   m3 = m1.template triangularView<Lower>();
  89   VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps));
  90
  91   // check M * inv(L) using in place API
  92   m4 = m3;
  93   m1.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4);
  94   VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Lower>(), m3);
  95
  96   // check M * inv(U) using in place API
  97   m3 = m1.template triangularView<Upper>();
  98   m4 = m3;
  99   m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4);
 100   VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Upper>(), m3);
 101
 102   // check solve with unit diagonal
 103   m3 = m1.template triangularView<UnitUpper>();
 104   VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));
 105
 106 //   VERIFY((  m1.template triangularView<Upper>()
 107 //           * m2.template triangularView<Upper>()).isUpperTriangular());
 108
 109   // test swap
 110   m1.setOnes();
 111   m2.setZero();
 112   m2.template triangularView<Upper>().swap(m1);
 113   m3.setZero();
 114   m3.template triangularView<Upper>().setOnes();
 115   VERIFY_IS_APPROX(m2,m3);
 116
 117 }
 118
 119
 120 template<typename MatrixType> void triangular_rect(const MatrixType& m)
 121 {
 122   typedef const typename MatrixType::Index Index;
 123   typedef typename MatrixType::Scalar Scalar;
 124   typedef typename NumTraits<Scalar>::Real RealScalar;
 125   enum { Rows =  MatrixType::RowsAtCompileTime, Cols =  MatrixType::ColsAtCompileTime };
 126
 127   Index rows = m.rows();
 128   Index cols = m.cols();
 129
 130   MatrixType m1 = MatrixType::Random(rows, cols),
 131              m2 = MatrixType::Random(rows, cols),
 132              m3(rows, cols),
 133              m4(rows, cols),
 134              r1(rows, cols),
 135              r2(rows, cols);
 136
 137   MatrixType m1up = m1.template triangularView<Upper>();
 138   MatrixType m2up = m2.template triangularView<Upper>();
 139
 140   if (rows>1 && cols>1)
 141   {
 142     VERIFY(m1up.isUpperTriangular());
 143     VERIFY(m2up.transpose().isLowerTriangular());
 144     VERIFY(!m2.isLowerTriangular());
 145   }
 146
 147   // test overloaded operator+=
 148   r1.setZero();
 149   r2.setZero();
 150   r1.template triangularView<Upper>() +=  m1;
 151   r2 += m1up;
 152   VERIFY_IS_APPROX(r1,r2);
 153
 154   // test overloaded operator=
 155   m1.setZero();
 156   m1.template triangularView<Upper>() = 3 * m2;
 157   m3 = 3 * m2;
 158   VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
 159
 160
 161   m1.setZero();
 162   m1.template triangularView<Lower>() = 3 * m2;
 163   VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
 164
 165   m1.setZero();
 166   m1.template triangularView<StrictlyUpper>() = 3 * m2;
 167   VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
 168
 169
 170   m1.setZero();
 171   m1.template triangularView<StrictlyLower>() = 3 * m2;
 172   VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
 173   m1.setRandom();
 174   m2 = m1.template triangularView<Upper>();
 175   VERIFY(m2.isUpperTriangular());
 176   VERIFY(!m2.isLowerTriangular());
 177   m2 = m1.template triangularView<StrictlyUpper>();
 178   VERIFY(m2.isUpperTriangular());
 179   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 180   m2 = m1.template triangularView<UnitUpper>();
 181   VERIFY(m2.isUpperTriangular());
 182   m2.diagonal().array() -= Scalar(1);
 183   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 184   m2 = m1.template triangularView<Lower>();
 185   VERIFY(m2.isLowerTriangular());
 186   VERIFY(!m2.isUpperTriangular());
 187   m2 = m1.template triangularView<StrictlyLower>();
 188   VERIFY(m2.isLowerTriangular());
 189   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 190   m2 = m1.template triangularView<UnitLower>();
 191   VERIFY(m2.isLowerTriangular());
 192   m2.diagonal().array() -= Scalar(1);
 193   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 194   // test swap
 195   m1.setOnes();
 196   m2.setZero();
 197   m2.template triangularView<Upper>().swap(m1);
 198   m3.setZero();
 199   m3.template triangularView<Upper>().setOnes();
 200   VERIFY_IS_APPROX(m2,m3);
 201 }
 202
 203 void bug_159()
 204 {
 205   Matrix3d m = Matrix3d::Random().triangularView<Lower>();
 206   EIGEN_UNUSED_VARIABLE(m)
 207 }
 208
 209 void test_triangular()
 210 {
 211   int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20);
 212   for(int i = 0; i < g_repeat ; i++)
 213   {
 214     int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r)
 215     int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c)
 216
 217     CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
 218     CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
 219     CALL_SUBTEST_3( triangular_square(Matrix3d()) );
 220     CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
 221     CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
 222     CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
 223
 224     CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
 225     CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
 226     CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
 227     CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
 228     CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
 229   }
 230
 231   CALL_SUBTEST_1( bug_159() );
 232 }