12 #ifndef EIGEN_COMPLEX_SCHUR_H
13 #define EIGEN_COMPLEX_SCHUR_H
15 #include "./HessenbergDecomposition.h"
20 template<
typename MatrixType,
bool IsComplex>
struct complex_schur_reduce_to_hessenberg;
54 typedef _MatrixType MatrixType;
56 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
57 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
58 Options = MatrixType::Options,
59 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
60 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64 typedef typename MatrixType::Scalar
Scalar;
66 typedef typename MatrixType::Index Index;
98 m_isInitialized(false),
99 m_matUisUptodate(false)
112 : m_matT(matrix.rows(),matrix.cols()),
113 m_matU(matrix.rows(),matrix.cols()),
114 m_hess(matrix.rows()),
115 m_isInitialized(false),
116 m_matUisUptodate(false)
137 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
138 eigen_assert(m_matUisUptodate &&
"The matrix U has not been computed during the ComplexSchur decomposition.");
161 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
192 eigen_assert(m_isInitialized &&
"RealSchur is not initialized.");
206 bool m_isInitialized;
207 bool m_matUisUptodate;
210 bool subdiagonalEntryIsNeglegible(Index i);
212 void reduceToTriangularForm(
bool computeU);
213 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType,
NumTraits<
Scalar>::IsComplex>;
219 template<typename MatrixType>
220 inline bool
ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
222 RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
223 RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
234 template<
typename MatrixType>
237 if (iter == 10 || iter == 20)
240 return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
245 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
246 RealScalar normt = t.cwiseAbs().sum();
257 if(internal::norm1(eival1) > internal::norm1(eival2))
258 eival2 = det / eival1;
260 eival1 = det / eival2;
263 if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
264 return normt * eival1;
266 return normt * eival2;
270 template<
typename MatrixType>
273 m_matUisUptodate =
false;
274 eigen_assert(matrix.cols() == matrix.rows());
276 if(matrix.cols() == 1)
278 m_matT = matrix.template cast<ComplexScalar>();
279 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
281 m_isInitialized =
true;
282 m_matUisUptodate = computeU;
286 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*
this, matrix, computeU);
287 reduceToTriangularForm(computeU);
294 template<
typename MatrixType,
bool IsComplex>
295 struct complex_schur_reduce_to_hessenberg
301 _this.m_matT = _this.m_hess.
matrixH();
302 if(computeU) _this.m_matU = _this.m_hess.
matrixQ();
306 template<
typename MatrixType>
307 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
309 static void run(ComplexSchur<MatrixType>& _this,
const MatrixType& matrix,
bool computeU)
315 _this.m_hess.compute(matrix);
316 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
320 MatrixType Q = _this.m_hess.matrixQ();
321 _this.m_matU = Q.template cast<ComplexScalar>();
329 template<
typename MatrixType>
330 void ComplexSchur<MatrixType>::reduceToTriangularForm(
bool computeU)
336 Index iu = m_matT.cols() - 1;
346 if(!subdiagonalEntryIsNeglegible(iu-1))
break;
361 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
370 ComplexScalar shift = computeShift(iu, iter);
371 JacobiRotation<ComplexScalar> rot;
372 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
373 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
374 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
375 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
377 for(Index i=il+1 ; i<iu ; i++)
379 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
381 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
382 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
383 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
392 m_isInitialized =
true;
393 m_matUisUptodate = computeU;
398 #endif // EIGEN_COMPLEX_SCHUR_H