Householder.h
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) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_HOUSEHOLDER_H
00027 #define EIGEN_HOUSEHOLDER_H
00028 
00029 namespace Eigen { 
00030 
00031 namespace internal {
00032 template<int n> struct decrement_size
00033 {
00034   enum {
00035     ret = n==Dynamic ? n : n-1
00036   };
00037 };
00038 }
00039 
00056 template<typename Derived>
00057 void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
00058 {
00059   VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
00060   makeHouseholder(essentialPart, tau, beta);
00061 }
00062 
00078 template<typename Derived>
00079 template<typename EssentialPart>
00080 void MatrixBase<Derived>::makeHouseholder(
00081   EssentialPart& essential,
00082   Scalar& tau,
00083   RealScalar& beta) const
00084 {
00085   EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
00086   VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
00087   
00088   RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
00089   Scalar c0 = coeff(0);
00090 
00091   if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0))
00092   {
00093     tau = RealScalar(0);
00094     beta = internal::real(c0);
00095     essential.setZero();
00096   }
00097   else
00098   {
00099     beta = internal::sqrt(internal::abs2(c0) + tailSqNorm);
00100     if (internal::real(c0)>=RealScalar(0))
00101       beta = -beta;
00102     essential = tail / (c0 - beta);
00103     tau = internal::conj((beta - c0) / beta);
00104   }
00105 }
00106 
00122 template<typename Derived>
00123 template<typename EssentialPart>
00124 void MatrixBase<Derived>::applyHouseholderOnTheLeft(
00125   const EssentialPart& essential,
00126   const Scalar& tau,
00127   Scalar* workspace)
00128 {
00129   if(rows() == 1)
00130   {
00131     *this *= Scalar(1)-tau;
00132   }
00133   else
00134   {
00135     Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
00136     Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
00137     tmp.noalias() = essential.adjoint() * bottom;
00138     tmp += this->row(0);
00139     this->row(0) -= tau * tmp;
00140     bottom.noalias() -= tau * essential * tmp;
00141   }
00142 }
00143 
00159 template<typename Derived>
00160 template<typename EssentialPart>
00161 void MatrixBase<Derived>::applyHouseholderOnTheRight(
00162   const EssentialPart& essential,
00163   const Scalar& tau,
00164   Scalar* workspace)
00165 {
00166   if(cols() == 1)
00167   {
00168     *this *= Scalar(1)-tau;
00169   }
00170   else
00171   {
00172     Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
00173     Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
00174     tmp.noalias() = right * essential.conjugate();
00175     tmp += this->col(0);
00176     this->col(0) -= tau * tmp;
00177     right.noalias() -= tau * tmp * essential.transpose();
00178   }
00179 }
00180 
00181 } // end namespace Eigen
00182 
00183 #endif // EIGEN_HOUSEHOLDER_H