HessenbergDecomposition.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H
27 #define EIGEN_HESSENBERGDECOMPOSITION_H
28 
29 namespace Eigen {
30 
31 namespace internal {
32 
33 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
34 template<typename MatrixType>
35 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
36 {
37  typedef MatrixType ReturnType;
38 };
39 
40 }
41 
72 template<typename _MatrixType> class HessenbergDecomposition
73 {
74  public:
75 
77  typedef _MatrixType MatrixType;
78 
79  enum {
80  Size = MatrixType::RowsAtCompileTime,
83  MaxSize = MatrixType::MaxRowsAtCompileTime,
85  };
86 
88  typedef typename MatrixType::Scalar Scalar;
89  typedef typename MatrixType::Index Index;
90 
98 
101 
102  typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
103 
116  : m_matrix(size,size),
117  m_temp(size),
118  m_isInitialized(false)
119  {
120  if(size>1)
121  m_hCoeffs.resize(size-1);
122  }
123 
134  : m_matrix(matrix),
135  m_temp(matrix.rows()),
136  m_isInitialized(false)
137  {
138  if(matrix.rows()<2)
139  {
140  m_isInitialized = true;
141  return;
142  }
143  m_hCoeffs.resize(matrix.rows()-1,1);
144  _compute(m_matrix, m_hCoeffs, m_temp);
145  m_isInitialized = true;
146  }
147 
166  {
167  m_matrix = matrix;
168  if(matrix.rows()<2)
169  {
170  m_isInitialized = true;
171  return *this;
172  }
173  m_hCoeffs.resize(matrix.rows()-1,1);
174  _compute(m_matrix, m_hCoeffs, m_temp);
175  m_isInitialized = true;
176  return *this;
177  }
178 
193  {
194  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
195  return m_hCoeffs;
196  }
197 
227  const MatrixType& packedMatrix() const
228  {
229  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
230  return m_matrix;
231  }
232 
248  {
249  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
250  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
251  .setLength(m_matrix.rows() - 1)
252  .setShift(1);
253  }
254 
276  {
277  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
278  return MatrixHReturnType(*this);
279  }
280 
281  private:
282 
284  typedef typename NumTraits<Scalar>::Real RealScalar;
285  static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
286 
287  protected:
292 };
293 
306 template<typename MatrixType>
307 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
308 {
309  assert(matA.rows()==matA.cols());
310  Index n = matA.rows();
311  temp.resize(n);
312  for (Index i = 0; i<n-1; ++i)
313  {
314  // let's consider the vector v = i-th column starting at position i+1
315  Index remainingSize = n-i-1;
316  RealScalar beta;
317  Scalar h;
318  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
319  matA.col(i).coeffRef(i+1) = beta;
320  hCoeffs.coeffRef(i) = h;
321 
322  // Apply similarity transformation to remaining columns,
323  // i.e., compute A = H A H'
324 
325  // A = H A
326  matA.bottomRightCorner(remainingSize, remainingSize)
327  .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
328 
329  // A = A H'
330  matA.rightCols(remainingSize)
331  .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
332  }
333 }
334 
335 namespace internal {
336 
352 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
353 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
354 {
355  typedef typename MatrixType::Index Index;
356  public:
361  HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
362 
368  template <typename ResultType>
369  inline void evalTo(ResultType& result) const
370  {
371  result = m_hess.packedMatrix();
372  Index n = result.rows();
373  if (n>2)
374  result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
375  }
376 
377  Index rows() const { return m_hess.packedMatrix().rows(); }
378  Index cols() const { return m_hess.packedMatrix().cols(); }
379 
380  protected:
381  const HessenbergDecomposition<MatrixType>& m_hess;
382 };
383 
384 } // end namespace internal
385 
386 } // end namespace Eigen
387 
388 #endif // EIGEN_HESSENBERGDECOMPOSITION_H