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 
42 namespace itpp
43 {
44 
45 #if defined(HAVE_LAPACK)
46 
47 bool 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 
75 bool 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 
104 bool 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 
111 bool 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 
119 mat schur(const mat &A)
120 {
121  mat U, T;
122  schur(A, U, T);
123  return T;
124 }
125 
126 
127 cmat schur(const cmat &A)
128 {
129  cmat U, T;
130  schur(A, U, T);
131  return T;
132 }
133 
134 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:19 for IT++ by Doxygen 1.8.2