1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file RealSchur.hpp
19 \author Manh Ha NGUYEN
20 \since 22 Jan 2013
21 \date 7 August 2013
22
23 \brief Definition of the xlifepp::RealSchur class
24
25 Class deals with real Schur decomposition.
26 */
27
28 // This file is adapted from Eigen, a lightweight C++ template library
29 // for linear algebra.
30 //
31 // Copyright (C) 2009 Claire Maurice
32 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
33 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
34 //
35 // This Source Code Form is subject to the terms of the Mozilla
36 // Public License v. 2.0. If a copy of the MPL was not distributed
37 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
38
39 #ifndef EIGEN_REAL_SCHUR_HPP
40 #define EIGEN_REAL_SCHUR_HPP
41
42 #include "utils.h"
43 #include "../utils/Jacobi.hpp"
44 #include "../utils/VectorEigenDense.hpp"
45 #include "./HessenbergDecomposition.hpp"
46 #include "./HouseholderQR.hpp"
47
48 namespace xlifepp
49 {
50
51 namespace internalEigenSolver
52 {
53 template<typename MatrixType>
54 void swapRealSchurInPlace(MatrixType& matrixT, MatrixType& matrixQ, unsigned int ifst, unsigned int ilst);
55 template<typename MatrixType>
56 MatrixType directSwapping(MatrixType& matT, unsigned int p, unsigned int q, real_t gamma);
57 }
58
59 /*!
60 \class RealSchur
61
62 \brief Performs a real Schur decomposition of a square matrix
63
64 \tparam _MatrixType the type of the matrix of which we are computing the
65 real Schur decomposition; this is expected to be an instantiation of the
66 Matrix class template.
67
68 Given a real square matrix A, this class computes the real Schur
69 decomposition: \f$ A = U T U^T \f$ where U is a real orthogonal matrix and
70 T is a real quasi-triangular matrix. An orthogonal matrix is a matrix whose
71 inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular
72 matrix is a block-triangular matrix whose diagonal consists of 1-by-1
73 blocks and 2-by-2 blocks with complex eigenvalues. The eigenvalues of the
74 blocks on the diagonal of T are the same as the eigenvalues of the matrix
75 A, and thus the real Schur decomposition is used in EigenSolver to compute
76 the eigendecomposition of a matrix.
77
78 Call the function compute() to compute the real Schur decomposition of a
79 given matrix. Alternatively, you can use the RealSchur(const MatrixType&, bool)
80 constructor which computes the real Schur decomposition at construction
81 time. Once the decomposition is computed, you can use the matrixU() and
82 matrixT() functions to retrieve the matrices U and T in the decomposition.
83
84 \note The implementation is adapted from
85 <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
86 Their code is based on EISPACK.
87
88 */
89 template<typename _MatrixType>
90 class RealSchur
91 {
92 public:
93 typedef _MatrixType MatrixType;
94 typedef typename MatrixType::type_t Scalar;
95 typedef typename NumTraits<Scalar>::ComplexScalar ComplexScalar;
96
97 typedef VectorEigenDense<ComplexScalar> EigenvalueType;
98 typedef VectorEigenDense<Scalar> ColumnVectorType;
99
100 /*!
101 \brief Default constructor.
102
103 \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed.
104
105 The default constructor is useful in cases in which the user intends to
106 perform decompositions via compute(). The \p size parameter is only
107 used as a hint. It is not an error to give a wrong \p size, but it may
108 impair performance.
109
110 \sa compute() for an example.
111 */
RealSchur(Index size)112 RealSchur(Index size)
113 : matT_(size, size),
114 matU_(size, size),
115 hess_(size),
116 maxIterations_(40),
117 isInitialized_(false),
118 matUisUptodate_(false),
119 info_(_noConvergence)
120 { }
121
122 /*!
123 \brief Constructor; computes real Schur decomposition of given matrix.
124
125 \param[in] matrix Square matrix whose Schur decomposition is to be computed.
126 \param[in] computeU If true, both T and U are computed; if false, only T is computed.
127
128 This constructor calls compute() to compute the Schur decomposition.
129 */
RealSchur(const MatrixType & matrix,bool computeU=true)130 RealSchur(const MatrixType& matrix, bool computeU = true)
131 : matT_(matrix.numOfRows(), matrix.numOfCols()),
132 matU_(matrix.numOfRows(), matrix.numOfCols()),
133 hess_(matrix.numOfRows()),
134 maxIterations_(40),
135 isInitialized_(false),
136 matUisUptodate_(false),
137 info_(_noConvergence)
138 {
139 compute(matrix, computeU);
140 }
141
142 /*!
143 \brief Returns the orthogonal matrix in the Schur decomposition.
144
145 \returns A const reference to the matrix U.
146
147 \pre Either the constructor RealSchur(const MatrixType&, bool) or the
148 member function compute(const MatrixType&, bool) has been called before
149 to compute the Schur decomposition of a matrix, and \p computeU was set
150 to true (the default value).
151
152 \sa RealSchur(const MatrixType&, bool) for an example
153 */
matrixU() const154 const MatrixType& matrixU() const
155 {
156 if (!isInitialized_) { error("eigensolver_not_initialized", "RealSchur"); }
157 if (!matUisUptodate_) { error("eigensolver_matrix_not_computed", "U", "Real"); }
158 return matU_;
159 }
160
161 /*!
162 \brief Returns the quasi-triangular matrix in the Schur decomposition.
163
164 \returns A const reference to the matrix T.
165
166 \pre Either the constructor RealSchur(const MatrixType&, bool) or the
167 member function compute(const MatrixType&, bool) has been called before
168 to compute the Schur decomposition of a matrix.
169 */
matrixT() const170 const MatrixType& matrixT() const
171 {
172 if (!isInitialized_) { error("eigensolver_not_initialized", "RealSchur"); }
173 return matT_;
174 }
175
176 /*!
177 \brief Computes Schur decomposition of given matrix.
178
179 \param[in] matrix Square matrix whose Schur decomposition is to be computed.
180 \param[in] computeU If true, both T and U are computed; if false, only T is computed.
181 \returns Reference to \c *this
182
183 The Schur decomposition is computed by first reducing the matrix to
184 Hessenberg form using the class HessenbergDecomposition. The Hessenberg
185 matrix is then reduced to triangular form by performing Francis QR
186 iterations with implicit double shift. The cost of computing the Schur
187 decomposition depends on the number of iterations; as a rough guide, it
188 may be taken to be \f$25n^3\f$ flops if \a computeU is true and
189 \f$10n^3\f$ flops if \a computeU is false.
190 */
191 RealSchur& compute(const MatrixType& matrix, bool computeU = true);
192
193 /*!
194 \brief Reports whether previous computation was successful.
195
196 \returns \c Success if computation was succesful, \c NoConvergence otherwise.
197 */
info() const198 ComputationInfo info() const
199 {
200 if (!isInitialized_) { error("eigensolver_not_initialized", "RealSchur"); }
201 return info_;
202 }
203
204 /*!
205 \brief Reorder block diagonal of quasi-triagular from row ifst to row with index ilst
206 */
207 void swapSchur(unsigned int ifst, unsigned int ilst, bool computeU);
208
209
210 private:
211 MatrixType matT_;
212 MatrixType matU_;
213 HessenbergDecomposition<MatrixType> hess_;
214
215 /*!
216 \brief Maximum number of iterations.
217
218 Maximum number of iterations allowed for an eigenvalue to converge.
219 Default value = 40 (taken from LAPACK)
220 */
221 int maxIterations_;
222
223 bool isInitialized_; //!< true if initialized
224 bool matUisUptodate_;
225 ComputationInfo info_;
226
227
228 typedef VectorEigenDense<Scalar> Vector3s;
229
230 Scalar computeNormOfT();
231 Index findSmallSubdiagEntry(Index iu, Scalar norm);
232 void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
233 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
234 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
235 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector);
236 };
237
238
239 template<typename MatrixType>
compute(const MatrixType & matrix,bool computeU)240 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
241 {
242 trace_p->push("RealSchur::compute");
243 if (matrix.numOfCols() != matrix.numOfRows()) { matrix.nonSquare("Computing real schur", matrix.numOfRows(), matrix.numOfCols()); }
244
245 // Step 1. Reduce to Hessenberg form
246 hess_.compute(matrix);
247 matT_ = hess_.matrixH();
248 if (computeU)
249 { matU_ = hess_.matrixQ(); }
250
251 // Step 2. Reduce to real Schur form
252 //m_workspaceVector.resize(matT_.numOfCols());
253 //Scalar* workspace = &m_workspaceVector.coeffRef(0);
254
255 // The matrix matT_ is divided in three parts.
256 // Rows 0,...,il-1 are decoupled from the rest because matT_(il,il-1) is zero.
257 // Rows il,...,iu is the part we are working on (the active window).
258 // Rows iu+1,...,end are already brought in triangular form.
259 Index iu = matT_.numOfCols() - 1;
260 Index iter = 0; // iteration count for current eigenvalue
261 Index totalIter = 0; // iteration count for whole matrix
262 Scalar exshift(0); // sum of exceptional shifts
263 Scalar norm = computeNormOfT();
264
265 if(norm != 0)
266 {
267 while (iu >= 0)
268 {
269 Index il = findSmallSubdiagEntry(iu, norm);
270
271 // Check for convergence
272 if (il == iu) // One root found
273 {
274 matT_.coeffRef(iu, iu) = matT_.coeff(iu, iu) + exshift;
275 if (iu > 0)
276 { matT_.coeffRef(iu, iu - 1) = Scalar(0); }
277 iu--;
278 iter = 0;
279 }
280 else if (il == iu - 1) // Two roots found
281 {
282 splitOffTwoRows(iu, computeU, exshift);
283 iu -= 2;
284 iter = 0;
285 }
286 else // No convergence yet
287 {
288 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
289 Vector3s firstHouseholderVector(3, 0.0), shiftInfo(3, 0.0);
290 computeShift(iu, iter, exshift, shiftInfo);
291 iter = iter + 1;
292 totalIter = totalIter + 1;
293 if (totalIter > maxIterations_ * matrix.numOfCols()) { break; }
294 Index im;
295 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
296 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector);
297 }
298 }
299 }
300 if(totalIter <= maxIterations_ * matrix.numOfCols())
301 { info_ = _success; }
302 else
303 { info_ = _noConvergence; }
304
305 isInitialized_ = true;
306 matUisUptodate_ = computeU;
307
308 trace_p->pop();
309 return *this;
310 }
311
312 /*! \internal Computes and returns vector L1 norm of T */
313 template<typename MatrixType>
computeNormOfT()314 typename MatrixType::type_t RealSchur<MatrixType>::computeNormOfT()
315 {
316 const Index size = matT_.numOfCols();
317 // FIXME to be efficient the following would requires a triangular reduxion code
318 // Scalar norm = matT_.upper().cwiseAbs().sum()
319 // + matT_.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
320 Scalar norm(0);
321 for (Index j = 0; j < size; ++j)
322 //norm += matT_.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
323 { norm += matT_.blockRow(j, std::max(j - 1, Index(0)), size - std::max(j - 1, Index(0))).sumAbs(); }
324 return norm;
325 }
326
327 /*! \internal Look for single small sub-diagonal element and returns its index */
328 template<typename MatrixType>
findSmallSubdiagEntry(Index iu,Scalar norm)329 Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
330 {
331 Index res = iu;
332 while (res > 0)
333 {
334 Scalar s = std::abs(matT_.coeff(res - 1, res - 1)) + std::abs(matT_.coeff(res, res));
335 if (s == 0.0)
336 { s = norm; }
337 if (std::abs(matT_.coeff(res, res - 1)) < (NumTraits<Scalar>::epsilon() * s))
338 { break; }
339 res--;
340 }
341 return res;
342 }
343
344 /*! \internal Update T given that rows iu-1 and iu decouple from the rest. */
345 template<typename MatrixType>
splitOffTwoRows(Index iu,bool computeU,Scalar exshift)346 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
347 {
348 const Index size = matT_.numOfCols();
349 const Index row = matT_.numOfRows();
350
351 // The eigenvalues of the 2x2 matrix [a b; c d] are
352 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
353 Scalar p = Scalar(0.5) * (matT_.coeff(iu - 1, iu - 1) - matT_.coeff(iu, iu));
354 Scalar q = p * p + matT_.coeff(iu, iu - 1) * matT_.coeff(iu - 1, iu); // q = tr^2 / 4 - det = discr/4
355 matT_.coeffRef(iu, iu) += exshift;
356 matT_.coeffRef(iu - 1, iu - 1) += exshift;
357
358 if (q >= Scalar(0)) // Two real eigenvalues
359 {
360 Scalar z = std::sqrt(std::abs(q));
361 JacobiRotation<Scalar> rot;
362 if (p >= Scalar(0))
363 { rot.makeGivens(p + z, matT_.coeff(iu, iu - 1)); }
364 else
365 { rot.makeGivens(p - z, matT_.coeff(iu, iu - 1)); }
366
367 // matT_.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
368 MatrixType rightCols(matT_, 0, iu - 1, row, size - iu + 1);
369 rightCols.applyOnTheLeft(iu - 1, iu, rot.adjoint());
370 matT_.replace(rightCols, 0, iu - 1, row, size - iu + 1);
371 // matT_.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
372 MatrixType topRows(matT_, 0, 0, iu + 1, size);
373 topRows.applyOnTheRight(iu - 1, iu, rot);
374 matT_.replace(topRows, 0, 0, iu + 1, size);
375
376 matT_.coeffRef(iu, iu - 1) = Scalar(0);
377
378 if (computeU)
379 { matU_.applyOnTheRight(iu - 1, iu, rot); }
380 }
381
382 if (iu > 1)
383 { matT_.coeffRef(iu - 1, iu - 2) = Scalar(0); }
384 }
385
386 /*! \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */
387 template<typename MatrixType>
computeShift(Index iu,Index iter,Scalar & exshift,Vector3s & shiftInfo)388 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
389 {
390 shiftInfo.coeffRef(0) = matT_.coeff(iu, iu);
391 shiftInfo.coeffRef(1) = matT_.coeff(iu - 1, iu - 1);
392 shiftInfo.coeffRef(2) = matT_.coeff(iu, iu - 1) * matT_.coeff(iu - 1, iu);
393
394 // Wilkinson's original ad hoc shift
395 if (iter == 10)
396 {
397 exshift += shiftInfo.coeff(0);
398 for (Index i = 0; i <= iu; ++i)
399 { matT_.coeffRef(i, i) -= shiftInfo.coeff(0); }
400 Scalar s = std::abs(matT_.coeff(iu, iu - 1)) + std::abs(matT_.coeff(iu - 1, iu - 2));
401 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
402 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
403 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
404 }
405
406 // MATLAB's new ad hoc shift
407 if (iter == 30)
408 {
409 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
410 s = s * s + shiftInfo.coeff(2);
411 if (s > Scalar(0))
412 {
413 s = std::sqrt(s);
414 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
415 { s = -s; }
416 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
417 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
418 exshift += s;
419 for (Index i = 0; i <= iu; ++i)
420 { matT_.coeffRef(i, i) -= s; }
421 //shiftInfo.setConstant(Scalar(0.964));
422 shiftInfo.set(Scalar(0.964));
423 }
424 }
425 }
426
427 /*! \internal Compute index im at which Francis QR step starts and the first Householder vector. */
428 template<typename MatrixType>
initFrancisQRStep(Index il,Index iu,const Vector3s & shiftInfo,Index & im,Vector3s & firstHouseholderVector)429 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
430 {
431 Vector3s& v = firstHouseholderVector; // alias to save typing
432
433 for (im = iu - 2; im >= il; --im)
434 {
435 const Scalar Tmm = matT_.coeff(im, im);
436 const Scalar r = shiftInfo.coeff(0) - Tmm;
437 const Scalar s = shiftInfo.coeff(1) - Tmm;
438 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / matT_.coeff(im + 1, im) + matT_.coeff(im, im + 1);
439 v.coeffRef(1) = matT_.coeff(im + 1, im + 1) - Tmm - r - s;
440 v.coeffRef(2) = matT_.coeff(im + 2, im + 1);
441 if (im == il)
442 {
443 break;
444 }
445 const Scalar lhs = matT_.coeff(im, im - 1) * (std::abs(v.coeff(1)) + std::abs(v.coeff(2)));
446 const Scalar rhs = v.coeff(0) * (std::abs(matT_.coeff(im - 1, im - 1)) + std::abs(Tmm) + std::abs(matT_.coeff(im + 1, im + 1)));
447 if (std::abs(lhs) < (NumTraits<Scalar>::epsilon() * rhs))
448 {
449 break;
450 }
451 }
452 }
453
454 /*! \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
455 template<typename MatrixType>
performFrancisQRStep(Index il,Index im,Index iu,bool computeU,const Vector3s & firstHouseholderVector)456 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector)
457 {
458 if (im < il) { error("is_lesser", im, il); }
459 if (im > iu-2) { error("is_greater", im, iu-2); }
460 const Index size = matT_.numOfCols();
461
462 for (Index k = im; k <= iu - 2; ++k)
463 {
464 bool firstIteration = (k == im);
465
466 Vector3s v(3);
467 if (firstIteration)
468 { v = firstHouseholderVector; }
469 else
470 //v = matT_.template block<3,1>(k,k-1);
471 { v = matT_.blockCol(k, k - 1, 3); }
472
473 Scalar tau, beta;
474 ColumnVectorType ess(2);
475 //Matrix<Scalar, 2, 1> ess;
476 v.makeHouseHolder(ess, tau, beta);
477
478 if (beta != Scalar(0)) // if v is not zero
479 {
480 if (firstIteration && k > il)
481 { matT_.coeffRef(k, k - 1) = -matT_.coeff(k, k - 1); }
482 else if (!firstIteration)
483 { matT_.coeffRef(k, k - 1) = beta; }
484
485 // These Householder transformations form the O(n^3) part of the algorithm
486 //matT_.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
487 MatrixType block1(matT_, k, k, 3, size - k);
488 block1.applyHouseholderOnTheLeft(ess, tau);
489 matT_.replace(block1, k, k, 3, size - k);
490 //matT_.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
491 MatrixType block2(matT_, 0, k, std::min(iu, k + 3) + 1, 3);
492 block2.applyHouseholderOnTheRight(ess, tau);
493 matT_.replace(block2, 0, k, std::min(iu, k + 3) + 1, 3);
494 if (computeU)
495 {
496 //matU_.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
497 MatrixType blockU(matU_, 0, k, size, 3);
498 blockU.applyHouseholderOnTheRight(ess, tau);
499 matU_.replace(blockU, 0, k, size, 3);
500 }
501 }
502 }
503
504 //Matrix<Scalar, 2, 1> v = matT_.template block<2,1>(iu-1, iu-2);
505 ColumnVectorType v = matT_.blockCol(iu - 1, iu - 2, 2);
506 Scalar tau, beta;
507 ColumnVectorType ess(1);
508 //Matrix<Scalar, 1, 1> ess;
509 v.makeHouseHolder(ess, tau, beta);
510
511 if (beta != Scalar(0)) // if v is not zero
512 {
513 matT_.coeffRef(iu - 1, iu - 2) = beta;
514 // matT_.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
515 // matT_.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
516 MatrixType block1(matT_, iu - 1, iu - 1, 2, size - iu + 1);
517 block1.applyHouseholderOnTheLeft(ess, tau);
518 matT_.replace(block1, iu - 1, iu - 1, 2, size - iu + 1);
519 MatrixType block2(matT_, 0, iu - 1, iu + 1, 2);
520 block2.applyHouseholderOnTheRight(ess, tau);
521 matT_.replace(block2, 0, iu - 1, iu + 1, 2);
522
523 if (computeU)
524 {
525 //matU_.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
526 MatrixType blockU(matU_, 0, iu - 1, size, 2);
527 blockU.applyHouseholderOnTheRight(ess, tau);
528 matU_.replace(blockU, 0, iu - 1, size, 2);
529 }
530 }
531
532 // clean up pollution due to round-off errors
533 for (Index i = im + 2; i <= iu; ++i)
534 {
535 matT_.coeffRef(i, i - 2) = Scalar(0);
536 if (i > im + 2)
537 { matT_.coeffRef(i, i - 3) = Scalar(0); }
538 }
539 }
540
541 template<typename MatrixType>
swapSchur(unsigned int ifst,unsigned int ilst,bool computeU)542 void RealSchur<MatrixType>::swapSchur(unsigned int ifst, unsigned int ilst, bool computeU)
543 {
544 internalEigenSolver::swapRealSchurInPlace(matT_,matU_,ifst, ilst);
545 }
546
547 namespace internalEigenSolver {
548 /*!
549 * \brief Given the quasi-triangular matrix T and orthogonal matrix Q obtained from
550 * the real Schur decomposition, this function reorders the eigenvalues appearing on the (block) diagonal of matrix T
551 * The diagonal element of block of T with row index ifst is moved to row ilst
552 *
553 * This is implemented from "On Swapping Diagonal Blocks in Real Schur Form"
554 * \param[in,out] matrixT quasi-triangular matrix of real Schur decomposition
555 * \param[in,out] matrixQ orthogonal of real Schur decomposition
556 * \param[in] ifst index of row needs moving
557 * \param[in] ilst index of destination row
558 */
559 template<typename MatrixType>
swapRealSchurInPlace(MatrixType & matrixT,MatrixType & matrixQ,unsigned int ifst,unsigned int ilst)560 void swapRealSchurInPlace(MatrixType& matrixT, MatrixType& matrixQ, unsigned int ifst, unsigned int ilst)
561 {
562 if (ifst != ilst) {
563 bool moveDown = (ifst < ilst);
564 unsigned int size = matrixT.numOfRows();
565 unsigned int idxf = std::min(ifst,ilst);
566 unsigned int idxl = std::max(ifst,ilst);//moveDown ? (ilst-ifst+1):(ifst-ilst+1);
567 unsigned int idx = idxf;
568
569 std::vector<unsigned int> blockSize;
570 while (idx <= idxl) {
571 if ((idx != (size-1)) && (NumTraits<real_t>::prec() < std::abs(matrixT.coeff(idx+1,idx)))) {
572 blockSize.push_back(2);
573 idx += 2;
574 } else {
575 blockSize.push_back(1);
576 ++idx;
577 }
578 }
579
580 int blockNum = blockSize.size();
581
582
583 // Very naive value, just to make sure gamma <= 1
584 real_t gamma = (std::abs(matrixT.coeff(idxf+1,idxf+1)+matrixT.coeff(idxf,idxf)))/(std::abs(matrixT.coeff(idxf+1,idxf+1)+std::abs(matrixT.coeff(idxf,idxf))));
585
586 if (moveDown) {
587 idx = idxf;
588 for (int i = 0; i < blockNum-1; ++i) {
589 JacobiRotation<real_t> rot;
590 if ( 2 == (blockSize[i]+blockSize[i+1])) {
591 if (NumTraits<real_t>::epsilon() < (std::abs(matrixT.coeff(idx+1,idx+1)- matrixT.coeff(idx,idx)))) {
592 rot.makeGivens(matrixT.coeff(idx,idx+1), (matrixT.coeff(idx+1,idx+1) - matrixT.coeff(idx,idx)));
593 matrixT.applyOnTheRight(idx,idx+1,rot);
594 matrixT.applyOnTheLeft(idx,idx+1,rot.adjoint());
595 matrixQ.applyOnTheRight(idx,idx+1,rot);
596 }
597 matrixT.coeffRef(idx+1,idx) = real_t(0);
598
599 ++idx;
600 }
601 else {
602 MatrixType T(matrixT,idx,idx,blockSize[i]+blockSize[i+1],blockSize[i]+blockSize[i+1]);
603 MatrixType Q = directSwapping(T, blockSize[i], blockSize[i+1], gamma);
604 MatrixType tempTcol(matrixT,0,idx,size,blockSize[i]+blockSize[i+1]);
605 MatrixType tmp(tempTcol);
606 multMatMat(tempTcol,Q,tmp);
607 matrixT.replace(tmp,0,idx,size,blockSize[i]+blockSize[i+1]);
608
609 MatrixType tempQ(matrixQ,0,idx,size,blockSize[i]+blockSize[i+1]);
610 multMatMat(tempQ,Q,tmp);
611 matrixQ.replace(tmp,0,idx,size,blockSize[i]+blockSize[i+1]);
612
613 MatrixType tempTrow(matrixT,idx,0,blockSize[i]+blockSize[i+1],size);
614 tmp.reshape(blockSize[i]+blockSize[i+1],size);
615 multMatMat(transpose(Q),tempTrow,tmp);
616 matrixT.replace(tmp,idx,0,blockSize[i]+blockSize[i+1],size);
617
618 if (1 == blockSize[i+1]) {
619 matrixT.coeffRef(idx+1,idx) = real_t(0);
620 matrixT.coeffRef(idx+2,idx) = real_t(0);
621 } else {
622 matrixT.coeffRef(idx+2,idx) = real_t(0);
623 matrixT.coeffRef(idx+2,idx+1) = real_t(0);
624 if (2 == blockSize[i]) {
625 matrixT.coeffRef(idx+3,idx) = real_t(0);
626 matrixT.coeffRef(idx+3,idx+1) = real_t(0);
627 }
628 }
629 idx += 2;
630
631 }
632 std::swap(blockSize[i],blockSize[i+1]);
633 }
634 }
635 else {
636 idx = idxl;
637 for (int i = blockNum-2; i >=0 ; --i) {
638 JacobiRotation<real_t> rot;
639 if ( 2 == (blockSize[i]+blockSize[i+1])) {
640 --idx;
641 if (NumTraits<real_t>::epsilon() < (std::abs(matrixT.coeff(idx+1,idx+1)- matrixT.coeff(idx,idx)))) {
642 rot.makeGivens(matrixT.coeff(idx,idx+1), (matrixT.coeff(idx+1,idx+1) - matrixT.coeff(idx,idx)));
643 matrixT.applyOnTheRight(idx,idx+1,rot);
644 matrixT.applyOnTheLeft(idx,idx+1,rot.adjoint());
645 matrixQ.applyOnTheRight(idx,idx+1,rot);
646 }
647 matrixT.coeffRef(idx+1,idx) = real_t(0);
648 }
649 else {
650 if (2 == blockSize[i]) idx -= 2;
651 else --idx;
652 MatrixType T(matrixT,idx,idx,blockSize[i]+blockSize[i+1],blockSize[i]+blockSize[i+1]);
653 MatrixType Q = directSwapping(T, blockSize[i], blockSize[i+1], gamma);
654 MatrixType tempTcol(matrixT,0,idx,size,blockSize[i]+blockSize[i+1]);
655 MatrixType tmp(tempTcol);
656 multMatMat(tempTcol,Q,tmp);
657 matrixT.replace(tmp,0,idx,size,blockSize[i]+blockSize[i+1]);
658
659 MatrixType tempQ(matrixQ,0,idx,size,blockSize[i]+blockSize[i+1]);
660 multMatMat(tempQ,Q,tmp);
661 matrixQ.replace(tmp,0,idx,size,blockSize[i]+blockSize[i+1]);
662
663 MatrixType tempTrow(matrixT,idx,0,blockSize[i]+blockSize[i+1],size);
664 tmp.reshape(blockSize[i]+blockSize[i+1],size);
665 multMatMat(transpose(Q),tempTrow,tmp);
666 matrixT.replace(tmp,idx,0,blockSize[i]+blockSize[i+1],size);
667 if (1 == blockSize[i+1]) {
668 matrixT.coeffRef(idx+1,idx) = real_t(0);
669 matrixT.coeffRef(idx+2,idx) = real_t(0);
670 } else {
671 matrixT.coeffRef(idx+2,idx) = real_t(0);
672 matrixT.coeffRef(idx+2,idx+1) = real_t(0);
673 if (2 == blockSize[i]) {
674 matrixT.coeffRef(idx+3,idx) = real_t(0);
675 matrixT.coeffRef(idx+3,idx+1) = real_t(0);
676 }
677 }
678 }
679 std::swap(blockSize[i],blockSize[i+1]);
680 }
681 }
682 }
683 }
684
685 template<typename MatrixType>
directSwapping(MatrixType & matT,unsigned int p,unsigned int q,real_t gamma)686 MatrixType directSwapping(MatrixType& matT, unsigned int p, unsigned int q, real_t gamma)
687 {
688 MatrixType eyeP(p), eyeQ(q);
689 eyeMatrix(eyeP); eyeMatrix(eyeQ);
690 MatrixType A11(matT,0,0,p,p);
691 MatrixType A12(matT,0,p,p,q);
692 MatrixType A22(matT,p,p,q,q);
693 A22 = transpose(A22);
694 MatrixType Kron = kroneckerProduct(eyeQ,A11);
695 MatrixType Kron2 = kroneckerProduct(A22,eyeP);
696 Kron -= Kron2;
697
698 std::vector<typename MatrixType::Scalar> X(p*q);
699 for (unsigned int j = 0; j < q; ++j) {
700 for (unsigned int i = 0; i < p; ++i) {
701 X[j*p+i] = gamma*A12.coeff(i,j);
702 }
703 }
704
705 typename MatrixType::Scalar piv;
706 number_t row;
707 gaussSolver(Kron,X,piv,row);
708 unsigned int qrSize = p+q;
709 MatrixType QR(qrSize,q);
710 for (unsigned int j = 0; j < q; ++j) {
711 for (unsigned int i = 0; i < p; ++i) {
712 QR.coeffRef(i,j) = -X[j*p+i];
713 }
714 QR.coeffRef(p+j,j) = gamma;
715 }
716
717 HouseholderQR<MatrixType> qr(QR);
718 MatrixType matQ = qr.matrixQ();
719
720 return matQ;
721 }
722
723 } // end namespace internalEigenSolver
724
725 } // end namespace xlifepp
726
727 #endif // EIGEN_REAL_SCHUR_HPP
728
729