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

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2009-2010 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/QR>
  12
  13 template<typename MatrixType> void householder(const MatrixType& m)
  14 {
  15   typedef typename MatrixType::Index Index;
  16   static bool even = true;
  17   even = !even;
  18   /* this test covers the following files:
  19      Householder.h
  20   */
  21   Index rows = m.rows();
  22   Index cols = m.cols();
  23
  24   typedef typename MatrixType::Scalar Scalar;
  25   typedef typename NumTraits<Scalar>::Real RealScalar;
  26   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  27   typedef Matrix<Scalar, internal::decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
  28   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
  29   typedef Matrix<Scalar, Dynamic, MatrixType::ColsAtCompileTime> HBlockMatrixType;
  30   typedef Matrix<Scalar, Dynamic, 1> HCoeffsVectorType;
  31
  32   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType;
  33
  34   Matrix<Scalar, EIGEN_SIZE_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp((std::max)(rows,cols));
  35   Scalar* tmp = &_tmp.coeffRef(0,0);
  36
  37   Scalar beta;
  38   RealScalar alpha;
  39   EssentialVectorType essential;
  40
  41   VectorType v1 = VectorType::Random(rows), v2;
  42   v2 = v1;
  43   v1.makeHouseholder(essential, beta, alpha);
  44   v1.applyHouseholderOnTheLeft(essential,beta,tmp);
  45   VERIFY_IS_APPROX(v1.norm(), v2.norm());
  46   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(v1.tail(rows-1).norm(), v1.norm());
  47   v1 = VectorType::Random(rows);
  48   v2 = v1;
  49   v1.applyHouseholderOnTheLeft(essential,beta,tmp);
  50   VERIFY_IS_APPROX(v1.norm(), v2.norm());
  51
  52   MatrixType m1(rows, cols),
  53              m2(rows, cols);
  54
  55   v1 = VectorType::Random(rows);
  56   if(even) v1.tail(rows-1).setZero();
  57   m1.colwise() = v1;
  58   m2 = m1;
  59   m1.col(0).makeHouseholder(essential, beta, alpha);
  60   m1.applyHouseholderOnTheLeft(essential,beta,tmp);
  61   VERIFY_IS_APPROX(m1.norm(), m2.norm());
  62   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
  63   VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m1(0,0)), numext::real(m1(0,0)));
  64   VERIFY_IS_APPROX(numext::real(m1(0,0)), alpha);
  65
  66   v1 = VectorType::Random(rows);
  67   if(even) v1.tail(rows-1).setZero();
  68   SquareMatrixType m3(rows,rows), m4(rows,rows);
  69   m3.rowwise() = v1.transpose();
  70   m4 = m3;
  71   m3.row(0).makeHouseholder(essential, beta, alpha);
  72   m3.applyHouseholderOnTheRight(essential,beta,tmp);
  73   VERIFY_IS_APPROX(m3.norm(), m4.norm());
  74   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
  75   VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m3(0,0)), numext::real(m3(0,0)));
  76   VERIFY_IS_APPROX(numext::real(m3(0,0)), alpha);
  77
  78   // test householder sequence on the left with a shift
  79
  80   Index shift = internal::random<Index>(0, std::max<Index>(rows-2,0));
  81   Index brows = rows - shift;
  82   m1.setRandom(rows, cols);
  83   HBlockMatrixType hbm = m1.block(shift,0,brows,cols);
  84   HouseholderQR<HBlockMatrixType> qr(hbm);
  85   m2 = m1;
  86   m2.block(shift,0,brows,cols) = qr.matrixQR();
  87   HCoeffsVectorType hc = qr.hCoeffs().conjugate();
  88   HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc);
  89   hseq.setLength(hc.size()).setShift(shift);
  90   VERIFY(hseq.length() == hc.size());
  91   VERIFY(hseq.shift() == shift);
  92
  93   MatrixType m5 = m2;
  94   m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
  95   VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
  96   m3 = hseq;
  97   VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying
  98
  99   SquareMatrixType hseq_mat = hseq;
 100   SquareMatrixType hseq_mat_conj = hseq.conjugate();
 101   SquareMatrixType hseq_mat_adj = hseq.adjoint();
 102   SquareMatrixType hseq_mat_trans = hseq.transpose();
 103   SquareMatrixType m6 = SquareMatrixType::Random(rows, rows);
 104   VERIFY_IS_APPROX(hseq_mat.adjoint(),    hseq_mat_adj);
 105   VERIFY_IS_APPROX(hseq_mat.conjugate(),  hseq_mat_conj);
 106   VERIFY_IS_APPROX(hseq_mat.transpose(),  hseq_mat_trans);
 107   VERIFY_IS_APPROX(hseq_mat * m6,             hseq_mat * m6);
 108   VERIFY_IS_APPROX(hseq_mat.adjoint() * m6,   hseq_mat_adj * m6);
 109   VERIFY_IS_APPROX(hseq_mat.conjugate() * m6, hseq_mat_conj * m6);
 110   VERIFY_IS_APPROX(hseq_mat.transpose() * m6, hseq_mat_trans * m6);
 111   VERIFY_IS_APPROX(m6 * hseq_mat,             m6 * hseq_mat);
 112   VERIFY_IS_APPROX(m6 * hseq_mat.adjoint(),   m6 * hseq_mat_adj);
 113   VERIFY_IS_APPROX(m6 * hseq_mat.conjugate(), m6 * hseq_mat_conj);
 114   VERIFY_IS_APPROX(m6 * hseq_mat.transpose(), m6 * hseq_mat_trans);
 115
 116   // test householder sequence on the right with a shift
 117
 118   TMatrixType tm2 = m2.transpose();
 119   HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc);
 120   rhseq.setLength(hc.size()).setShift(shift);
 121   VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
 122   m3 = rhseq;
 123   VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying
 124 }
 125
 126 void test_householder()
 127 {
 128   for(int i = 0; i < g_repeat; i++) {
 129     CALL_SUBTEST_1( householder(Matrix<double,2,2>()) );
 130     CALL_SUBTEST_2( householder(Matrix<float,2,3>()) );
 131     CALL_SUBTEST_3( householder(Matrix<double,3,5>()) );
 132     CALL_SUBTEST_4( householder(Matrix<float,4,4>()) );
 133     CALL_SUBTEST_5( householder(MatrixXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
 134     CALL_SUBTEST_6( householder(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
 135     CALL_SUBTEST_7( householder(MatrixXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
 136     CALL_SUBTEST_8( householder(Matrix<double,1,1>()) );
 137   }
 138 }