Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
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 }
00182
00183 #endif // EIGEN_HOUSEHOLDER_H