1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <limits> 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Eigenvalues> 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void verifyIsQuasiTriangular(const MatrixType& T) 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = T.cols(); 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check T is lower Hessenberg 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int row = 2; row < size; ++row) { 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col = 0; col < row - 1; ++col) { 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(T(row,col) == Scalar(0)); 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check that any non-zero on the subdiagonal is followed by a zero and is 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // part of a 2x2 diagonal block with imaginary eigenvalues. 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int row = 1; row < size; ++row) { 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (T(row,row-1) != Scalar(0)) { 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(row == size-1 || T(row+1,row) == 0); 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tr = T(row-1,row-1) + T(row,row); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar det = T(row-1,row-1) * T(row,row) - T(row-1,row) * T(row,row-1); 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(4 * det > tr * tr); 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTime) 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test basic functionality: T is quasi-triangular and A = U T U* 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int counter = 0; counter < g_repeat; ++counter) { 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType A = MatrixType::Random(size, size); 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> schurOfA(A); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(schurOfA.info(), Success); 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType U = schurOfA.matrixU(); 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType T = schurOfA.matrixT(); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath verifyIsQuasiTriangular(T); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(A, U * T * U.transpose()); 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test asserts when not initialized 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> rsUninitialized; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(rsUninitialized.matrixT()); 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(rsUninitialized.matrixU()); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(rsUninitialized.info()); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test whether compute() and constructor returns same result 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType A = MatrixType::Random(size, size); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> rs1; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rs1.compute(A); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> rs2(A); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rs1.info(), Success); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rs2.info(), Success); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT()); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU()); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Test maximum number of iterations 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealSchur<MatrixType> rs3; 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rs3.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A); 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.info(), Success); 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.matrixT(), rs1.matrixT()); 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.matrixU(), rs1.matrixU()); 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (size > 2) { 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rs3.setMaxIterations(1).compute(A); 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.info(), NoConvergence); 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.getMaxIterations(), 1); 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType Atriangular = A; 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Atriangular.template triangularView<StrictlyLower>().setZero(); 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rs3.setMaxIterations(1).compute(Atriangular); // triangular matrices do not need any iterations 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.info(), Success); 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.matrixT(), Atriangular); 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(rs3.matrixU(), MatrixType::Identity(size, size)); 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test computation of only T, not U 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> rsOnlyT(A, false); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rsOnlyT.info(), Success); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rs1.matrixT(), rsOnlyT.matrixT()); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(rsOnlyT.matrixU()); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (size > 2) 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test matrix with NaN 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A(0,0) = std::numeric_limits<typename MatrixType::Scalar>::quiet_NaN(); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> rsNaN(A); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(rsNaN.info(), NoConvergence); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_schur_real() 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(( schur<Matrix4f>() )); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2(( schur<MatrixXd>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4)) )); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3(( schur<Matrix<float, 1, 1> >() )); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4(( schur<Matrix<double, 3, 3, Eigen::RowMajor> >() )); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test problem size constructors 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_5(RealSchur<MatrixXf>(10)); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 113