IT++ Logo
schur.cpp
Go to the documentation of this file.
1
29#ifndef _MSC_VER
30# include <itpp/config.h>
31#else
32# include <itpp/config_msvc.h>
33#endif
34
35#if defined(HAVE_LAPACK)
36# include <itpp/base/algebra/lapack.h>
37#endif
38
40
41
42namespace itpp
43{
44
45#if defined(HAVE_LAPACK)
46
47bool schur(const mat &A, mat &U, mat &T)
48{
49 it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square");
50
51 char jobvs = 'V';
52 char sort = 'N';
53 int info;
54 int n = A.rows();
55 int lda = n;
56 int ldvs = n;
57 int lwork = 3 * n; // This may be choosen better!
58 int sdim = 0;
59 vec wr(n);
60 vec wi(n);
61 vec work(lwork);
62
63 T.set_size(lda, n, false);
64 U.set_size(ldvs, n, false);
65
66 T = A; // The routine overwrites input matrix with eigenvectors
67
68 dgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, wr._data(), wi._data(),
69 U._data(), &ldvs, work._data(), &lwork, 0, &info);
70
71 return (info == 0);
72}
73
74
75bool schur(const cmat &A, cmat &U, cmat &T)
76{
77 it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square");
78
79 char jobvs = 'V';
80 char sort = 'N';
81 int info;
82 int n = A.rows();
83 int lda = n;
84 int ldvs = n;
85 int lwork = 2 * n; // This may be choosen better!
86 int sdim = 0;
87 vec rwork(n);
88 cvec w(n);
89 cvec work(lwork);
90
91 T.set_size(lda, n, false);
92 U.set_size(ldvs, n, false);
93
94 T = A; // The routine overwrites input matrix with eigenvectors
95
96 zgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, w._data(), U._data(),
97 &ldvs, work._data(), &lwork, rwork._data(), 0, &info);
98
99 return (info == 0);
100}
101
102#else
103
104bool schur(const mat &A, mat &U, mat &T)
105{
106 it_error("LAPACK library is needed to use schur() function");
107 return false;
108}
109
110
111bool schur(const cmat &A, cmat &U, cmat &T)
112{
113 it_error("LAPACK library is needed to use schur() function");
114 return false;
115}
116
117#endif // HAVE_LAPACK
118
119mat schur(const mat &A)
120{
121 mat U, T;
122 schur(A, U, T);
123 return T;
124}
125
126
127cmat schur(const cmat &A)
128{
129 cmat U, T;
130 schur(A, U, T);
131 return T;
132}
133
134} // namespace itpp
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
bool schur(const mat &A, mat &U, mat &T)
Schur decomposition of a real matrix.
Definition: schur.cpp:104
itpp namespace
Definition: itmex.h:37
Definitions of Schur decomposition functions.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.5