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))<RealScalar(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 substitution 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   m1.setRandom();
 118   m3 = m1.template triangularView<Upper>();
 119   Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> m5(cols, internal::random<int>(1,20));  m5.setRandom();
 120   Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> m6(internal::random<int>(1,20), rows);  m6.setRandom();
 121   VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3*m5);
 122   VERIFY_IS_APPROX(m6*m1.template triangularView<Upper>(), m6*m3);
 123
 124   m1up = m1.template triangularView<Upper>();
 125   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
 126   VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
 127   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
 128   VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
 129
 130   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal());
 131
 132 }
 133
 134
 135 template<typename MatrixType> void triangular_rect(const MatrixType& m)
 136 {
 137   typedef const typename MatrixType::Index Index;
 138   typedef typename MatrixType::Scalar Scalar;
 139   typedef typename NumTraits<Scalar>::Real RealScalar;
 140   enum { Rows =  MatrixType::RowsAtCompileTime, Cols =  MatrixType::ColsAtCompileTime };
 141
 142   Index rows = m.rows();
 143   Index cols = m.cols();
 144
 145   MatrixType m1 = MatrixType::Random(rows, cols),
 146              m2 = MatrixType::Random(rows, cols),
 147              m3(rows, cols),
 148              m4(rows, cols),
 149              r1(rows, cols),
 150              r2(rows, cols);
 151
 152   MatrixType m1up = m1.template triangularView<Upper>();
 153   MatrixType m2up = m2.template triangularView<Upper>();
 154
 155   if (rows>1 && cols>1)
 156   {
 157     VERIFY(m1up.isUpperTriangular());
 158     VERIFY(m2up.transpose().isLowerTriangular());
 159     VERIFY(!m2.isLowerTriangular());
 160   }
 161
 162   // test overloaded operator+=
 163   r1.setZero();
 164   r2.setZero();
 165   r1.template triangularView<Upper>() +=  m1;
 166   r2 += m1up;
 167   VERIFY_IS_APPROX(r1,r2);
 168
 169   // test overloaded operator=
 170   m1.setZero();
 171   m1.template triangularView<Upper>() = 3 * m2;
 172   m3 = 3 * m2;
 173   VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
 174
 175
 176   m1.setZero();
 177   m1.template triangularView<Lower>() = 3 * m2;
 178   VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
 179
 180   m1.setZero();
 181   m1.template triangularView<StrictlyUpper>() = 3 * m2;
 182   VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
 183
 184
 185   m1.setZero();
 186   m1.template triangularView<StrictlyLower>() = 3 * m2;
 187   VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
 188   m1.setRandom();
 189   m2 = m1.template triangularView<Upper>();
 190   VERIFY(m2.isUpperTriangular());
 191   VERIFY(!m2.isLowerTriangular());
 192   m2 = m1.template triangularView<StrictlyUpper>();
 193   VERIFY(m2.isUpperTriangular());
 194   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 195   m2 = m1.template triangularView<UnitUpper>();
 196   VERIFY(m2.isUpperTriangular());
 197   m2.diagonal().array() -= Scalar(1);
 198   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 199   m2 = m1.template triangularView<Lower>();
 200   VERIFY(m2.isLowerTriangular());
 201   VERIFY(!m2.isUpperTriangular());
 202   m2 = m1.template triangularView<StrictlyLower>();
 203   VERIFY(m2.isLowerTriangular());
 204   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 205   m2 = m1.template triangularView<UnitLower>();
 206   VERIFY(m2.isLowerTriangular());
 207   m2.diagonal().array() -= Scalar(1);
 208   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
 209   // test swap
 210   m1.setOnes();
 211   m2.setZero();
 212   m2.template triangularView<Upper>().swap(m1);
 213   m3.setZero();
 214   m3.template triangularView<Upper>().setOnes();
 215   VERIFY_IS_APPROX(m2,m3);
 216 }
 217
 218 void bug_159()
 219 {
 220   Matrix3d m = Matrix3d::Random().triangularView<Lower>();
 221   EIGEN_UNUSED_VARIABLE(m)
 222 }
 223
 224 void test_triangular()
 225 {
 226   int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20);
 227   for(int i = 0; i < g_repeat ; i++)
 228   {
 229     int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r)
 230     int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c)
 231
 232     CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
 233     CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
 234     CALL_SUBTEST_3( triangular_square(Matrix3d()) );
 235     CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
 236     CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
 237     CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
 238
 239     CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
 240     CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
 241     CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
 242     CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
 243     CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
 244   }
 245
 246   CALL_SUBTEST_1( bug_159() );
 247 }