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;
98 m_isInitialized(false),
99 m_matUisUptodate(false),
112 template<
typename InputType>
114 : m_matT(matrix.rows(),matrix.cols()),
115 m_matU(matrix.rows(),matrix.cols()),
116 m_hess(matrix.rows()),
117 m_isInitialized(false),
118 m_matUisUptodate(false),
121 compute(matrix.
derived(), computeU);
140 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
141 eigen_assert(m_matUisUptodate &&
"The matrix U has not been computed during the ComplexSchur decomposition.");
164 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
190 template<
typename InputType>
210 template<
typename HessMatrixType,
typename OrthMatrixType>
211 ComplexSchur& computeFromHessenberg(
const HessMatrixType& matrixH,
const OrthMatrixType& matrixQ,
bool computeU=
true);
219 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
230 m_maxIters = maxIters;
245 static const int m_maxIterationsPerRow = 30;
248 ComplexMatrixType m_matT, m_matU;
251 bool m_isInitialized;
252 bool m_matUisUptodate;
256 bool subdiagonalEntryIsNeglegible(Index i);
257 ComplexScalar computeShift(Index iu, Index iter);
258 void reduceToTriangularForm(
bool computeU);
259 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType,
NumTraits<Scalar>::IsComplex>;
265 template<typename MatrixType>
266 inline bool
ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
268 RealScalar d = numext::norm1(m_matT.
coeff(i,i)) + numext::norm1(m_matT.
coeff(i+1,i+1));
269 RealScalar sd = numext::norm1(m_matT.
coeff(i+1,i));
272 m_matT.
coeffRef(i+1,i) = ComplexScalar(0);
280 template<
typename MatrixType>
284 if (iter == 10 || iter == 20)
287 return abs(numext::real(m_matT.
coeff(iu,iu-1))) +
abs(numext::real(m_matT.
coeff(iu-1,iu-2)));
293 RealScalar normt = t.
cwiseAbs().sum();
296 ComplexScalar b = t.
coeff(0,1) * t.
coeff(1,0);
297 ComplexScalar c = t.
coeff(0,0) - t.
coeff(1,1);
298 ComplexScalar disc =
sqrt(c*c + RealScalar(4)*b);
299 ComplexScalar det = t.
coeff(0,0) * t.
coeff(1,1) - b;
300 ComplexScalar trace = t.
coeff(0,0) + t.
coeff(1,1);
301 ComplexScalar eival1 = (trace + disc) / RealScalar(2);
302 ComplexScalar eival2 = (trace - disc) / RealScalar(2);
303 RealScalar eival1_norm = numext::norm1(eival1);
304 RealScalar eival2_norm = numext::norm1(eival2);
307 if(eival1_norm > eival2_norm)
308 eival2 = det / eival1;
309 else if(eival2_norm!=RealScalar(0))
310 eival1 = det / eival2;
313 if(numext::norm1(eival1-t.
coeff(1,1)) < numext::norm1(eival2-t.
coeff(1,1)))
314 return normt * eival1;
316 return normt * eival2;
320 template<
typename MatrixType>
321 template<
typename InputType>
324 m_matUisUptodate =
false;
325 eigen_assert(matrix.
cols() == matrix.
rows());
327 if(matrix.
cols() == 1)
329 m_matT = matrix.
derived().template cast<ComplexScalar>();
330 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
332 m_isInitialized =
true;
333 m_matUisUptodate = computeU;
337 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*
this, matrix.
derived(), computeU);
338 computeFromHessenberg(m_matT, m_matU, computeU);
342 template<
typename MatrixType>
343 template<
typename HessMatrixType,
typename OrthMatrixType>
349 reduceToTriangularForm(computeU);
355 template<
typename MatrixType,
bool IsComplex>
356 struct complex_schur_reduce_to_hessenberg
362 _this.m_matT = _this.m_hess.
matrixH();
363 if(computeU) _this.m_matU = _this.m_hess.
matrixQ();
367 template<
typename MatrixType>
368 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
376 _this.m_matT = _this.m_hess.
matrixH().template cast<ComplexScalar>();
380 MatrixType Q = _this.m_hess.
matrixQ();
381 _this.m_matU = Q.template cast<ComplexScalar>();
389 template<
typename MatrixType>
392 Index maxIters = m_maxIters;
394 maxIters = m_maxIterationsPerRow * m_matT.rows();
400 Index iu = m_matT.cols() - 1;
410 if(!subdiagonalEntryIsNeglegible(iu-1))
break;
421 if(totalIter > maxIters)
break;
425 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
434 ComplexScalar shift = computeShift(iu, iter);
437 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.
adjoint());
438 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
439 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
441 for(Index i=il+1 ; i<iu ; i++)
444 m_matT.
coeffRef(i+1,i-1) = ComplexScalar(0);
445 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.
adjoint());
446 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
447 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
451 if(totalIter <= maxIters)
456 m_isInitialized =
true;
457 m_matUisUptodate = computeU;
462 #endif // EIGEN_COMPLEX_SCHUR_H ComplexSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:113
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:228
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Rotation given by a cosine-sine pair.
Definition: Jacobi.h:34
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Derived & derived()
Definition: EigenBase.h:46
ComplexSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU=true)
Compute Schur decomposition from a given Hessenberg matrix.
std::complex< RealScalar > ComplexScalar
Complex scalar type for _MatrixType.
Definition: ComplexSchur.h:74
ComplexSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
Definition: EigenBase.h:29
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > ComplexMatrixType
Type for the matrices in the Schur decomposition.
Definition: ComplexSchur.h:81
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: ComplexSchur.h:64
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:175
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:152
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
Definition: HessenbergDecomposition.h:262
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
ComplexSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
Definition: Constants.h:442
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:152
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:138
Eigen::Index Index
Definition: ComplexSchur.h:66
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Definition: HessenbergDecomposition.h:234
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
const int Dynamic
Definition: Constants.h:22
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:51
ComputationInfo
Definition: Constants.h:440
const CwiseAbsReturnType cwiseAbs() const
Definition: MatrixBase.h:34
Definition: Constants.h:446
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162
JacobiRotation adjoint() const
Definition: Jacobi.h:67
void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:162