1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#include "main.h"
11#include <limits>
12#include <Eigen/Eigenvalues>
13
14template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTime)
15{
16  typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
17  typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
18
19  // Test basic functionality: T is triangular and A = U T U*
20  for(int counter = 0; counter < g_repeat; ++counter) {
21    MatrixType A = MatrixType::Random(size, size);
22    ComplexSchur<MatrixType> schurOfA(A);
23    VERIFY_IS_EQUAL(schurOfA.info(), Success);
24    ComplexMatrixType U = schurOfA.matrixU();
25    ComplexMatrixType T = schurOfA.matrixT();
26    for(int row = 1; row < size; ++row) {
27      for(int col = 0; col < row; ++col) {
28	VERIFY(T(row,col) == (typename MatrixType::Scalar)0);
29      }
30    }
31    VERIFY_IS_APPROX(A.template cast<ComplexScalar>(), U * T * U.adjoint());
32  }
33
34  // Test asserts when not initialized
35  ComplexSchur<MatrixType> csUninitialized;
36  VERIFY_RAISES_ASSERT(csUninitialized.matrixT());
37  VERIFY_RAISES_ASSERT(csUninitialized.matrixU());
38  VERIFY_RAISES_ASSERT(csUninitialized.info());
39
40  // Test whether compute() and constructor returns same result
41  MatrixType A = MatrixType::Random(size, size);
42  ComplexSchur<MatrixType> cs1;
43  cs1.compute(A);
44  ComplexSchur<MatrixType> cs2(A);
45  VERIFY_IS_EQUAL(cs1.info(), Success);
46  VERIFY_IS_EQUAL(cs2.info(), Success);
47  VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
48  VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
49
50  // Test maximum number of iterations
51  ComplexSchur<MatrixType> cs3;
52  cs3.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A);
53  VERIFY_IS_EQUAL(cs3.info(), Success);
54  VERIFY_IS_EQUAL(cs3.matrixT(), cs1.matrixT());
55  VERIFY_IS_EQUAL(cs3.matrixU(), cs1.matrixU());
56  cs3.setMaxIterations(1).compute(A);
57  VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success);
58  VERIFY_IS_EQUAL(cs3.getMaxIterations(), 1);
59
60  MatrixType Atriangular = A;
61  Atriangular.template triangularView<StrictlyLower>().setZero();
62  cs3.setMaxIterations(1).compute(Atriangular); // triangular matrices do not need any iterations
63  VERIFY_IS_EQUAL(cs3.info(), Success);
64  VERIFY_IS_EQUAL(cs3.matrixT(), Atriangular.template cast<ComplexScalar>());
65  VERIFY_IS_EQUAL(cs3.matrixU(), ComplexMatrixType::Identity(size, size));
66
67  // Test computation of only T, not U
68  ComplexSchur<MatrixType> csOnlyT(A, false);
69  VERIFY_IS_EQUAL(csOnlyT.info(), Success);
70  VERIFY_IS_EQUAL(cs1.matrixT(), csOnlyT.matrixT());
71  VERIFY_RAISES_ASSERT(csOnlyT.matrixU());
72
73  if (size > 1)
74  {
75    // Test matrix with NaN
76    A(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
77    ComplexSchur<MatrixType> csNaN(A);
78    VERIFY_IS_EQUAL(csNaN.info(), NoConvergence);
79  }
80}
81
82void test_schur_complex()
83{
84  CALL_SUBTEST_1(( schur<Matrix4cd>() ));
85  CALL_SUBTEST_2(( schur<MatrixXcf>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4)) ));
86  CALL_SUBTEST_3(( schur<Matrix<std::complex<float>, 1, 1> >() ));
87  CALL_SUBTEST_4(( schur<Matrix<float, 3, 3, Eigen::RowMajor> >() ));
88
89  // Test problem size constructors
90  CALL_SUBTEST_5(ComplexSchur<MatrixXf>(10));
91}
92