1 // This is core/vnl/algo/vnl_generalized_schur.h
2 #ifndef vnl_generalized_schur_h_
3 #define vnl_generalized_schur_h_
4 //:
5 // \file
6 // \brief Solves the generalized eigenproblem det(t A - s B) = 0.
7 // \author fsm, Oxford RRG
8 // \date 2 Oct 2001
9
10 #include <algorithm>
11 #include <vnl/vnl_matrix.h>
12 #include <vnl/vnl_vector.h>
13 #include <vnl/algo/vnl_algo_export.h>
14
15 //:
16 // For a \e real scalar type T, this function uses orthogonal
17 // matrices L, R to reduce the (square) matrices A, B to generalized
18 // (real) Schur form. This means that B is upper triangular and A is
19 // block upper triangular with blocks of size at most 2x2 such that
20 // the 2x2 blocks B corresponding to 2x2 blocks of A are diagonal.
21 // E.g.:
22 // \verbatim
23 // [ * * * * * ]
24 // [ * * * * ]
25 // A <- L^* A R = [ * * * * ]
26 // [ * * ]
27 // [ * * ]
28 //
29 // [ * * * * * ]
30 // [ * * * ]
31 // B <- L^* B R = [ * * * ]
32 // [ * ]
33 // [ * ]
34 // \endverbatim
35 //
36 // In addition, the function computes the generalized eigenvalues
37 // (alphar(k) + i alphai(k) : beta(k) for k = 0, 1, 2,...
38 template <class T>
39 bool vnl_generalized_schur(vnl_matrix<T> *A,
40 vnl_matrix<T> *B,
41 vnl_vector<T> *alphar,
42 vnl_vector<T> *alphai,
43 vnl_vector<T> *beta,
44 vnl_matrix<T> *L,
45 vnl_matrix<T> *R);
46
47 template <>
48 VNL_ALGO_EXPORT bool vnl_generalized_schur(vnl_matrix<double> *A,
49 vnl_matrix<double> *B,
50 vnl_vector<double> *alphar,
51 vnl_vector<double> *alphai,
52 vnl_vector<double> *beta,
53 vnl_matrix<double> *L,
54 vnl_matrix<double> *R);
55
56 #ifdef _MSC_VER
57 # include <vcl_msvc_warnings.h>
58 #endif
59
60 template <class T>
vnl_generalized_schur_convert_cast(double a)61 T vnl_generalized_schur_convert_cast(double a) { return static_cast<T>(a); }
62
63 template <class T>
vnl_generalized_schur(vnl_matrix<T> * A,vnl_matrix<T> * B,vnl_vector<T> * alphar,vnl_vector<T> * alphai,vnl_vector<T> * beta,vnl_matrix<T> * L,vnl_matrix<T> * R)64 inline bool vnl_generalized_schur(vnl_matrix<T> *A,
65 vnl_matrix<T> *B,
66 vnl_vector<T> *alphar,
67 vnl_vector<T> *alphai,
68 vnl_vector<T> *beta,
69 vnl_matrix<T> *L,
70 vnl_matrix<T> *R)
71 {
72 vnl_matrix<double> A_(A->rows(), A->cols());
73 vnl_matrix<double> B_(B->rows(), B->cols());
74 std::copy(A->begin(), A->end(), A_.begin());
75 std::copy(B->begin(), B->end(), B_.begin());
76
77 vnl_vector<double> alphar_;
78 vnl_vector<double> alphai_;
79 vnl_vector<double> beta_;
80 vnl_matrix<double> L_;
81 vnl_matrix<double> R_;
82
83 if (! vnl_generalized_schur/*<double>*/(&A_, &B_, &alphar_, &alphai_, &beta_, &L_, &R_))
84 return false;
85
86 std::transform(A_.begin(), A_.end(), A->begin(), vnl_generalized_schur_convert_cast<T>);
87 std::transform(B_.begin(), B_.end(), B->begin(), vnl_generalized_schur_convert_cast<T>);
88
89 alphar->set_size(alphar_.size());
90 std::transform(alphar_.begin(), alphar_.end(), alphar->begin(), vnl_generalized_schur_convert_cast<T>);
91
92 alphai->set_size(alphai_.size());
93 std::transform(alphai_.begin(), alphai_.end(), alphai->begin(), vnl_generalized_schur_convert_cast<T>);
94
95 beta ->set_size(beta_ .size());
96 std::transform(beta_ .begin(), beta_ .end(), beta ->begin(), vnl_generalized_schur_convert_cast<T>);
97
98 L->set_size(L_.rows(), L_.cols());
99 std::transform(L_.begin(), L_.end(), L->begin(), vnl_generalized_schur_convert_cast<T>);
100
101 R->set_size(R_.rows(), R_.cols());
102 std::transform(R_.begin(), R_.end(), R->begin(), vnl_generalized_schur_convert_cast<T>);
103
104 return true;
105 }
106
107 #endif // vnl_generalized_schur_h_
108