householder.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
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 #include <Eigen/QR>
00027 
00028 template<typename MatrixType> void householder(const MatrixType& m)
00029 {
00030   typedef typename MatrixType::Index Index;
00031   static bool even = true;
00032   even = !even;
00033   /* this test covers the following files:
00034      Householder.h
00035   */
00036   Index rows = m.rows();
00037   Index cols = m.cols();
00038 
00039   typedef typename MatrixType::Scalar Scalar;
00040   typedef typename NumTraits<Scalar>::Real RealScalar;
00041   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00042   typedef Matrix<Scalar, internal::decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
00043   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
00044   typedef Matrix<Scalar, Dynamic, MatrixType::ColsAtCompileTime> HBlockMatrixType;
00045   typedef Matrix<Scalar, Dynamic, 1> HCoeffsVectorType;
00046 
00047   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RightSquareMatrixType;
00048   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic> VBlockMatrixType;
00049   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType;
00050   
00051   Matrix<Scalar, EIGEN_SIZE_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
00052   Scalar* tmp = &_tmp.coeffRef(0,0);
00053 
00054   Scalar beta;
00055   RealScalar alpha;
00056   EssentialVectorType essential;
00057 
00058   VectorType v1 = VectorType::Random(rows), v2;
00059   v2 = v1;
00060   v1.makeHouseholder(essential, beta, alpha);
00061   v1.applyHouseholderOnTheLeft(essential,beta,tmp);
00062   VERIFY_IS_APPROX(v1.norm(), v2.norm());
00063   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(v1.tail(rows-1).norm(), v1.norm());
00064   v1 = VectorType::Random(rows);
00065   v2 = v1;
00066   v1.applyHouseholderOnTheLeft(essential,beta,tmp);
00067   VERIFY_IS_APPROX(v1.norm(), v2.norm());
00068 
00069   MatrixType m1(rows, cols),
00070              m2(rows, cols);
00071 
00072   v1 = VectorType::Random(rows);
00073   if(even) v1.tail(rows-1).setZero();
00074   m1.colwise() = v1;
00075   m2 = m1;
00076   m1.col(0).makeHouseholder(essential, beta, alpha);
00077   m1.applyHouseholderOnTheLeft(essential,beta,tmp);
00078   VERIFY_IS_APPROX(m1.norm(), m2.norm());
00079   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
00080   VERIFY_IS_MUCH_SMALLER_THAN(internal::imag(m1(0,0)), internal::real(m1(0,0)));
00081   VERIFY_IS_APPROX(internal::real(m1(0,0)), alpha);
00082 
00083   v1 = VectorType::Random(rows);
00084   if(even) v1.tail(rows-1).setZero();
00085   SquareMatrixType m3(rows,rows), m4(rows,rows);
00086   m3.rowwise() = v1.transpose();
00087   m4 = m3;
00088   m3.row(0).makeHouseholder(essential, beta, alpha);
00089   m3.applyHouseholderOnTheRight(essential,beta,tmp);
00090   VERIFY_IS_APPROX(m3.norm(), m4.norm());
00091   if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
00092   VERIFY_IS_MUCH_SMALLER_THAN(internal::imag(m3(0,0)), internal::real(m3(0,0)));
00093   VERIFY_IS_APPROX(internal::real(m3(0,0)), alpha);
00094 
00095   // test householder sequence on the left with a shift
00096 
00097   Index shift = internal::random<Index>(0, std::max<Index>(rows-2,0));
00098   Index brows = rows - shift;
00099   m1.setRandom(rows, cols);
00100   HBlockMatrixType hbm = m1.block(shift,0,brows,cols);
00101   HouseholderQR<HBlockMatrixType> qr(hbm);
00102   m2 = m1;
00103   m2.block(shift,0,brows,cols) = qr.matrixQR();
00104   HCoeffsVectorType hc = qr.hCoeffs().conjugate();
00105   HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc);
00106   hseq.setLength(hc.size()).setShift(shift);
00107   VERIFY(hseq.length() == hc.size());
00108   VERIFY(hseq.shift() == shift);
00109 
00110   MatrixType m5 = m2;
00111   m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
00112   VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
00113   m3 = hseq;
00114   VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying
00115 
00116   // test householder sequence on the right with a shift
00117 
00118   TMatrixType tm2 = m2.transpose();
00119   HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc);
00120   rhseq.setLength(hc.size()).setShift(shift);
00121   VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
00122   m3 = rhseq;
00123   VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying
00124 }
00125 
00126 void test_householder()
00127 {
00128   for(int i = 0; i < g_repeat; i++) {
00129     CALL_SUBTEST_1( householder(Matrix<double,2,2>()) );
00130     CALL_SUBTEST_2( householder(Matrix<float,2,3>()) );
00131     CALL_SUBTEST_3( householder(Matrix<double,3,5>()) );
00132     CALL_SUBTEST_4( householder(Matrix<float,4,4>()) );
00133     CALL_SUBTEST_5( householder(MatrixXd(10,12)) );
00134     CALL_SUBTEST_6( householder(MatrixXcf(16,17)) );
00135     CALL_SUBTEST_7( householder(MatrixXf(25,7)) );
00136     CALL_SUBTEST_8( householder(Matrix<double,1,1>()) );
00137   }
00138 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:46