Eigen  3.3.3
Householder.h
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 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_HOUSEHOLDER_H
00012 #define EIGEN_HOUSEHOLDER_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017 template<int n> struct decrement_size
00018 {
00019   enum {
00020     ret = n==Dynamic ? n : n-1
00021   };
00022 };
00023 }
00024 
00041 template<typename Derived>
00042 void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
00043 {
00044   VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
00045   makeHouseholder(essentialPart, tau, beta);
00046 }
00047 
00063 template<typename Derived>
00064 template<typename EssentialPart>
00065 void MatrixBase<Derived>::makeHouseholder(
00066   EssentialPart& essential,
00067   Scalar& tau,
00068   RealScalar& beta) const
00069 {
00070   using std::sqrt;
00071   using numext::conj;
00072   
00073   EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
00074   VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
00075   
00076   RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
00077   Scalar c0 = coeff(0);
00078   const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
00079 
00080   if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol)
00081   {
00082     tau = RealScalar(0);
00083     beta = numext::real(c0);
00084     essential.setZero();
00085   }
00086   else
00087   {
00088     beta = sqrt(numext::abs2(c0) + tailSqNorm);
00089     if (numext::real(c0)>=RealScalar(0))
00090       beta = -beta;
00091     essential = tail / (c0 - beta);
00092     tau = conj((beta - c0) / beta);
00093   }
00094 }
00095 
00111 template<typename Derived>
00112 template<typename EssentialPart>
00113 void MatrixBase<Derived>::applyHouseholderOnTheLeft(
00114   const EssentialPart& essential,
00115   const Scalar& tau,
00116   Scalar* workspace)
00117 {
00118   if(rows() == 1)
00119   {
00120     *this *= Scalar(1)-tau;
00121   }
00122   else if(tau!=Scalar(0))
00123   {
00124     Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
00125     Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
00126     tmp.noalias() = essential.adjoint() * bottom;
00127     tmp += this->row(0);
00128     this->row(0) -= tau * tmp;
00129     bottom.noalias() -= tau * essential * tmp;
00130   }
00131 }
00132 
00148 template<typename Derived>
00149 template<typename EssentialPart>
00150 void MatrixBase<Derived>::applyHouseholderOnTheRight(
00151   const EssentialPart& essential,
00152   const Scalar& tau,
00153   Scalar* workspace)
00154 {
00155   if(cols() == 1)
00156   {
00157     *this *= Scalar(1)-tau;
00158   }
00159   else if(tau!=Scalar(0))
00160   {
00161     Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
00162     Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
00163     tmp.noalias() = right * essential.conjugate();
00164     tmp += this->col(0);
00165     this->col(0) -= tau * tmp;
00166     right.noalias() -= tau * tmp * essential.transpose();
00167   }
00168 }
00169 
00170 } // end namespace Eigen
00171 
00172 #endif // EIGEN_HOUSEHOLDER_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends