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;
94 ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
98 m_isInitialized(false),
99 m_matUisUptodate(false),
113 : m_matT(matrix.rows(),matrix.cols()),
114 m_matU(matrix.rows(),matrix.cols()),
115 m_hess(matrix.rows()),
116 m_isInitialized(false),
117 m_matUisUptodate(false),
120 compute(matrix, computeU);
139 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
140 eigen_assert(m_matUisUptodate &&
"The matrix U has not been computed during the ComplexSchur decomposition.");
163 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
189 ComplexSchur& compute(
const MatrixType& matrix,
bool computeU =
true);
208 template<
typename HessMatrixType,
typename OrthMatrixType>
209 ComplexSchur& computeFromHessenberg(
const HessMatrixType& matrixH,
const OrthMatrixType& matrixQ,
bool computeU=
true);
217 eigen_assert(m_isInitialized &&
"ComplexSchur is not initialized.");
228 m_maxIters = maxIters;
243 static const int m_maxIterationsPerRow = 30;
246 ComplexMatrixType m_matT, m_matU;
249 bool m_isInitialized;
250 bool m_matUisUptodate;
254 bool subdiagonalEntryIsNeglegible(Index i);
255 ComplexScalar computeShift(Index iu, Index iter);
256 void reduceToTriangularForm(
bool computeU);
257 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType,
NumTraits<Scalar>::IsComplex>;
263 template<typename MatrixType>
264 inline bool
ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
266 RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
267 RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
270 m_matT.coeffRef(i+1,i) = ComplexScalar(0);
278 template<
typename MatrixType>
282 if (iter == 10 || iter == 20)
285 return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
291 RealScalar normt = t.
cwiseAbs().sum();
294 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
295 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
296 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
297 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
298 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
299 ComplexScalar eival1 = (trace + disc) / RealScalar(2);
300 ComplexScalar eival2 = (trace - disc) / RealScalar(2);
302 if(numext::norm1(eival1) > numext::norm1(eival2))
303 eival2 = det / eival1;
305 eival1 = det / eival2;
308 if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
309 return normt * eival1;
311 return normt * eival2;
315 template<
typename MatrixType>
318 m_matUisUptodate =
false;
319 eigen_assert(matrix.cols() == matrix.rows());
321 if(matrix.cols() == 1)
323 m_matT = matrix.template cast<ComplexScalar>();
324 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
326 m_isInitialized =
true;
327 m_matUisUptodate = computeU;
331 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*
this, matrix, computeU);
332 computeFromHessenberg(m_matT, m_matU, computeU);
336 template<
typename MatrixType>
337 template<
typename HessMatrixType,
typename OrthMatrixType>
343 reduceToTriangularForm(computeU);
349 template<
typename MatrixType,
bool IsComplex>
350 struct complex_schur_reduce_to_hessenberg
356 _this.m_matT = _this.m_hess.
matrixH();
357 if(computeU) _this.m_matU = _this.m_hess.
matrixQ();
361 template<
typename MatrixType>
362 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
370 _this.m_matT = _this.m_hess.
matrixH().template cast<ComplexScalar>();
374 MatrixType Q = _this.m_hess.
matrixQ();
375 _this.m_matU = Q.template cast<ComplexScalar>();
383 template<
typename MatrixType>
386 Index maxIters = m_maxIters;
388 maxIters = m_maxIterationsPerRow * m_matT.rows();
394 Index iu = m_matT.cols() - 1;
404 if(!subdiagonalEntryIsNeglegible(iu-1))
break;
415 if(totalIter > maxIters)
break;
419 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
428 ComplexScalar shift = computeShift(iu, iter);
430 rot.
makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
431 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.
adjoint());
432 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
433 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
435 for(Index i=il+1 ; i<iu ; i++)
437 rot.
makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
438 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
439 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.
adjoint());
440 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
441 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
445 if(totalIter <= maxIters)
450 m_isInitialized =
true;
451 m_matUisUptodate = computeU;
456 #endif // EIGEN_COMPLEX_SCHUR_H ComplexSchur(const MatrixType &matrix, bool computeU=true)
Constructor; computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:112
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > ComplexMatrixType
Type for the matrices in the Schur decomposition.
Definition: ComplexSchur.h:81
ComplexSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:316
ComplexSchur(Index size=RowsAtCompileTime==Dynamic?1:RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:228
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
ComplexSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU=true)
Compute Schur decomposition from a given Hessenberg matrix.
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:226
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition: Jacobi.h:148
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
Definition: HessenbergDecomposition.h:260
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: ComplexSchur.h:64
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Definition: HessenbergDecomposition.h:232
Definition: Eigen_Colamd.h:50
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:233
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:137
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:161
std::complex< RealScalar > ComplexScalar
Complex scalar type for _MatrixType.
Definition: ComplexSchur.h:74
Definition: Constants.h:380
Definition: Constants.h:376
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:215
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:51
JacobiRotation adjoint() const
Definition: Jacobi.h:62
ComputationInfo
Definition: Constants.h:374
const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > > cwiseAbs() const
Definition: MatrixBase.h:22
HessenbergDecomposition & compute(const MatrixType &matrix)
Computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:150