re_vision: triangular.cpp Source File

Go to the documentation of this file.
00001 // This file is triangularView of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "main.h"
00026 
00027 
00028 
00029 template<typename MatrixType> void triangular_square(const MatrixType& m)
00030 {
00031   typedef typename MatrixType::Scalar Scalar;
00032   typedef typename NumTraits<Scalar>::Real RealScalar;
00033   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00034 
00035   RealScalar largerEps = 10*test_precision<RealScalar>();
00036 
00037   typename MatrixType::Index rows = m.rows();
00038   typename MatrixType::Index cols = m.cols();
00039 
00040   MatrixType m1 = MatrixType::Random(rows, cols),
00041              m2 = MatrixType::Random(rows, cols),
00042              m3(rows, cols),
00043              m4(rows, cols),
00044              r1(rows, cols),
00045              r2(rows, cols),
00046              mzero = MatrixType::Zero(rows, cols),
00047              mones = MatrixType::Ones(rows, cols),
00048              identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00049                               ::Identity(rows, rows),
00050              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00051                               ::Random(rows, rows);
00052   VectorType v1 = VectorType::Random(rows),
00053              v2 = VectorType::Random(rows),
00054              vzero = VectorType::Zero(rows);
00055 
00056   MatrixType m1up = m1.template triangularView<Upper>();
00057   MatrixType m2up = m2.template triangularView<Upper>();
00058 
00059   if (rows*cols>1)
00060   {
00061     VERIFY(m1up.isUpperTriangular());
00062     VERIFY(m2up.transpose().isLowerTriangular());
00063     VERIFY(!m2.isLowerTriangular());
00064   }
00065 
00066 //   VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2);
00067 
00068   // test overloaded operator+=
00069   r1.setZero();
00070   r2.setZero();
00071   r1.template triangularView<Upper>() +=  m1;
00072   r2 += m1up;
00073   VERIFY_IS_APPROX(r1,r2);
00074 
00075   // test overloaded operator=
00076   m1.setZero();
00077   m1.template triangularView<Upper>() = m2.transpose() + m2;
00078   m3 = m2.transpose() + m2;
00079   VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1);
00080 
00081   // test overloaded operator=
00082   m1.setZero();
00083   m1.template triangularView<Lower>() = m2.transpose() + m2;
00084   VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
00085 
00086   VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(),
00087                    m3.conjugate().template triangularView<Lower>().toDenseMatrix());
00088 
00089   m1 = MatrixType::Random(rows, cols);
00090   for (int i=0; i<rows; ++i)
00091     while (internal::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random<Scalar>();
00092 
00093   Transpose<MatrixType> trm4(m4);
00094   // test back and forward subsitution with a vector as the rhs
00095   m3 = m1.template triangularView<Upper>();
00096   VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps));
00097   m3 = m1.template triangularView<Lower>();
00098   VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps));
00099   m3 = m1.template triangularView<Upper>();
00100   VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps));
00101   m3 = m1.template triangularView<Lower>();
00102   VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps));
00103 
00104   // test back and forward subsitution with a matrix as the rhs
00105   m3 = m1.template triangularView<Upper>();
00106   VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps));
00107   m3 = m1.template triangularView<Lower>();
00108   VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps));
00109   m3 = m1.template triangularView<Upper>();
00110   VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps));
00111   m3 = m1.template triangularView<Lower>();
00112   VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps));
00113 
00114   // check M * inv(L) using in place API
00115   m4 = m3;
00116   m3.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4);
00117   VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>()));
00118 
00119   // check M * inv(U) using in place API
00120   m3 = m1.template triangularView<Upper>();
00121   m4 = m3;
00122   m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4);
00123   VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>()));
00124 
00125   // check solve with unit diagonal
00126   m3 = m1.template triangularView<UnitUpper>();
00127   VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));
00128 
00129 //   VERIFY((  m1.template triangularView<Upper>()
00130 //           * m2.template triangularView<Upper>()).isUpperTriangular());
00131 
00132   // test swap
00133   m1.setOnes();
00134   m2.setZero();
00135   m2.template triangularView<Upper>().swap(m1);
00136   m3.setZero();
00137   m3.template triangularView<Upper>().setOnes();
00138   VERIFY_IS_APPROX(m2,m3);
00139 
00140 }
00141 
00142 
00143 template<typename MatrixType> void triangular_rect(const MatrixType& m)
00144 {
00145   typedef const typename MatrixType::Index Index;
00146   typedef typename MatrixType::Scalar Scalar;
00147   typedef typename NumTraits<Scalar>::Real RealScalar;
00148   enum { Rows =  MatrixType::RowsAtCompileTime, Cols =  MatrixType::ColsAtCompileTime };
00149   typedef Matrix<Scalar, Rows, 1> VectorType;
00150   typedef Matrix<Scalar, Rows, Rows> RMatrixType;
00151   
00152 
00153   Index rows = m.rows();
00154   Index cols = m.cols();
00155 
00156   MatrixType m1 = MatrixType::Random(rows, cols),
00157              m2 = MatrixType::Random(rows, cols),
00158              m3(rows, cols),
00159              m4(rows, cols),
00160              r1(rows, cols),
00161              r2(rows, cols),
00162              mzero = MatrixType::Zero(rows, cols),
00163              mones = MatrixType::Ones(rows, cols);
00164   RMatrixType identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00165                               ::Identity(rows, rows),
00166               square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00167                               ::Random(rows, rows);
00168   VectorType v1 = VectorType::Random(rows),
00169              v2 = VectorType::Random(rows),
00170              vzero = VectorType::Zero(rows);
00171 
00172   MatrixType m1up = m1.template triangularView<Upper>();
00173   MatrixType m2up = m2.template triangularView<Upper>();
00174 
00175   if (rows*cols>1)
00176   {
00177     VERIFY(m1up.isUpperTriangular());
00178     VERIFY(m2up.transpose().isLowerTriangular());
00179     VERIFY(!m2.isLowerTriangular());
00180   }
00181 
00182   // test overloaded operator+=
00183   r1.setZero();
00184   r2.setZero();
00185   r1.template triangularView<Upper>() +=  m1;
00186   r2 += m1up;
00187   VERIFY_IS_APPROX(r1,r2);
00188 
00189   // test overloaded operator=
00190   m1.setZero();
00191   m1.template triangularView<Upper>() = 3 * m2;
00192   m3 = 3 * m2;
00193   VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
00194 
00195 
00196   m1.setZero();
00197   m1.template triangularView<Lower>() = 3 * m2;
00198   VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
00199 
00200   m1.setZero();
00201   m1.template triangularView<StrictlyUpper>() = 3 * m2;
00202   VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
00203 
00204 
00205   m1.setZero();
00206   m1.template triangularView<StrictlyLower>() = 3 * m2;
00207   VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
00208   m1.setRandom();
00209   m2 = m1.template triangularView<Upper>();
00210   VERIFY(m2.isUpperTriangular());
00211   VERIFY(!m2.isLowerTriangular());
00212   m2 = m1.template triangularView<StrictlyUpper>();
00213   VERIFY(m2.isUpperTriangular());
00214   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
00215   m2 = m1.template triangularView<UnitUpper>();
00216   VERIFY(m2.isUpperTriangular());
00217   m2.diagonal().array() -= Scalar(1);
00218   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
00219   m2 = m1.template triangularView<Lower>();
00220   VERIFY(m2.isLowerTriangular());
00221   VERIFY(!m2.isUpperTriangular());
00222   m2 = m1.template triangularView<StrictlyLower>();
00223   VERIFY(m2.isLowerTriangular());
00224   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
00225   m2 = m1.template triangularView<UnitLower>();
00226   VERIFY(m2.isLowerTriangular());
00227   m2.diagonal().array() -= Scalar(1);
00228   VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
00229   // test swap
00230   m1.setOnes();
00231   m2.setZero();
00232   m2.template triangularView<Upper>().swap(m1);
00233   m3.setZero();
00234   m3.template triangularView<Upper>().setOnes();
00235   VERIFY_IS_APPROX(m2,m3);
00236 }
00237 
00238 void bug_159()
00239 {
00240   Matrix3d m = Matrix3d::Random().triangularView<Lower>(); 
00241 }
00242 
00243 void test_triangular()
00244 {
00245   for(int i = 0; i < g_repeat ; i++)
00246   {
00247     int r = internal::random<int>(2,20); EIGEN_UNUSED_VARIABLE(r);
00248     int c = internal::random<int>(2,20); EIGEN_UNUSED_VARIABLE(c);
00249 
00250     CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
00251     CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
00252     CALL_SUBTEST_3( triangular_square(Matrix3d()) );
00253     CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
00254     CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
00255     CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
00256 
00257     CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
00258     CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
00259     CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
00260     CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
00261     CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
00262   }
00263   
00264   CALL_SUBTEST_1( bug_159() );
00265 }