fn_eig.hpp
Go to the documentation of this file.
1 // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
2 // Copyright (C) 2008-2011 Conrad Sanderson
3 // Copyright (C) 2009 Edmund Highcock
4 // Copyright (C) 2011 Stanislav Funiak
5 //
6 // This file is part of the Armadillo C++ library.
7 // It is provided without any warranty of fitness
8 // for any purpose. You can redistribute this file
9 // and/or modify it under the terms of the GNU
10 // Lesser General Public License (LGPL) as published
11 // by the Free Software Foundation, either version 3
12 // of the License or (at your option) any later version.
13 // (see http://www.opensource.org/licenses for more info)
14 
15 
18 
19 
20 //
21 // symmetric/hermitian matrices
22 //
23 
24 
26 template<typename T1>
27 inline
28 bool
29 eig_sym
30  (
34  )
35  {
38 
39  // unwrap_check not used as T1::elem_type and T1::pod_type may not be the same.
40  // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::eig_sym() makes a copy of X
41 
42  const bool status = auxlib::eig_sym(eigval, X);
43 
44  if(status == false)
45  {
46  eigval.reset();
47  arma_bad("eig_sym(): failed to converge", false);
48  }
49 
50  return status;
51  }
52 
53 
54 
56 template<typename T1>
57 inline
59 eig_sym
60  (
63  )
64  {
67 
69  const bool status = auxlib::eig_sym(out, X);
70 
71  if(status == false)
72  {
73  out.reset();
74  arma_bad("eig_sym(): failed to converge");
75  }
76 
77  return out;
78  }
79 
80 
81 
83 template<typename T1>
84 inline
85 bool
86 eig_sym
87  (
92  )
93  {
96 
97  arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_sym(): eigval is an alias of eigvec" );
98 
99  const bool status = auxlib::eig_sym(eigval, eigvec, X);
100 
101  if(status == false)
102  {
103  eigval.reset();
104  eigvec.reset();
105  arma_bad("eig_sym(): failed to converge", false);
106  }
107 
108  return status;
109  }
110 
111 
112 
113 //
114 // general matrices
115 //
116 
117 
118 
120 template<typename T1>
121 inline
122 bool
123 eig_gen
124  (
125  Col< std::complex<typename T1::pod_type> >& eigval,
126  Mat<typename T1::elem_type>& l_eigvec,
127  Mat<typename T1::elem_type>& r_eigvec,
130  )
131  {
133  arma_ignore(junk);
134 
136  (
137  ((&l_eigvec) == (&r_eigvec)),
138  "eig_gen(): l_eigvec is an alias of r_eigvec"
139  );
140 
142  (
143  (
144  (((void*)(&eigval)) == ((void*)(&l_eigvec)))
145  ||
146  (((void*)(&eigval)) == ((void*)(&r_eigvec)))
147  ),
148  "eig_gen(): eigval is an alias of l_eigvec or r_eigvec"
149  );
150 
151  const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b');
152 
153  if(status == false)
154  {
155  eigval.reset();
156  l_eigvec.reset();
157  r_eigvec.reset();
158  arma_bad("eig_gen(): failed to converge", false);
159  }
160 
161  return status;
162  }
163 
164 
165 
169 template<typename eT, typename T1>
170 inline
171 bool
172 eig_gen
173  (
174  Col< std::complex<eT> >& eigval,
175  Mat< std::complex<eT> >& eigvec,
176  const Base<eT, T1>& X,
177  const char side = 'r',
179  )
180  {
182  arma_ignore(junk);
183 
184  //std::cout << "real" << std::endl;
185 
186  arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
187 
188  Mat<eT> dummy_eigvec;
189  Mat<eT> tmp_eigvec;
190 
191  bool status;
192 
193  switch(side)
194  {
195  case 'r':
196  status = auxlib::eig_gen(eigval, dummy_eigvec, tmp_eigvec, X, side);
197  break;
198 
199  case 'l':
200  status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side);
201  break;
202 
203  default:
204  arma_stop("eig_gen(): parameter 'side' is invalid");
205  status = false;
206  }
207 
208  if(status == false)
209  {
210  eigval.reset();
211  eigvec.reset();
212  arma_bad("eig_gen(): failed to converge", false);
213  }
214  else
215  {
216  const uword n = eigval.n_elem;
217 
218  if(n > 0)
219  {
220  eigvec.set_size(n,n);
221 
222  for(uword j=0; j<n; ++j)
223  {
224  if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) )
225  {
226  // eigvec.col(j) = Mat< std::complex<eT> >( tmp_eigvec.col(j), tmp_eigvec.col(j+1) );
227  // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -tmp_eigvec.col(j+1) );
228 
229  for(uword i=0; i<n; ++i)
230  {
231  eigvec.at(i,j) = std::complex<eT>( tmp_eigvec.at(i,j), tmp_eigvec.at(i,j+1) );
232  eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eigvec.at(i,j+1) );
233  }
234 
235  ++j;
236  }
237  else
238  {
239  // eigvec.col(i) = tmp_eigvec.col(i);
240 
241  for(uword i=0; i<n; ++i)
242  {
243  eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0));
244  }
245 
246  }
247  }
248  }
249  }
250 
251  return status;
252  }
253 
254 
255 
259 template<typename T, typename T1>
260 inline
261 bool
262 eig_gen
263  (
264  Col<std::complex<T> >& eigval,
265  Mat<std::complex<T> >& eigvec,
266  const Base<std::complex<T>, T1>& X,
267  const char side = 'r',
269  )
270  {
272  arma_ignore(junk);
273 
274  //std::cout << "complex" << std::endl;
275 
276  arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
277 
278  Mat< std::complex<T> > dummy_eigvec;
279 
280  bool status;
281 
282  switch(side)
283  {
284  case 'r':
285  status = auxlib::eig_gen(eigval, dummy_eigvec, eigvec, X, side);
286  break;
287 
288  case 'l':
289  status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side);
290  break;
291 
292  default:
293  arma_stop("eig_gen(): parameter 'side' is invalid");
294  status = false;
295  }
296 
297  if(status == false)
298  {
299  eigval.reset();
300  eigvec.reset();
301  arma_bad("eig_gen(): failed to converge", false);
302  }
303 
304  return status;
305  }
306 
307 
308 
310 
void set_size(const uword in_elem)
change the matrix to have user specified dimensions (data is not preserved)
Definition: Mat_meat.hpp:4211
static bool eig_gen(Col< std::complex< T > > &eigval, Mat< T > &l_eigvec, Mat< T > &r_eigvec, const Base< T, T1 > &X, const char side)
arma_inline const T1 & conj(const Base< typename T1::pod_type, T1 > &A)
Definition: fn_elem.hpp:430
const uword n_elem
number of elements in the matrix (read-only)
Definition: Mat_bones.hpp:31
static bool eig_sym(Col< eT > &eigval, const Base< eT, T1 > &X)
immediate eigenvalues of a symmetric real matrix using LAPACK
u32 uword
Definition: typedef.hpp:85
Class for column vectors (matrices with only one column)
Definition: Col_bones.hpp:20
#define arma_debug_check
Definition: debug.hpp:1084
arma_inline arma_warn_unused eT & at(const uword i)
linear element accessor (treats the matrix as a vector); no bounds check.
Definition: Mat_meat.hpp:3692
#define arma_ignore(variable)
bool eig_gen(Col< std::complex< typename T1::pod_type > > &eigval, Mat< typename T1::elem_type > &l_eigvec, Mat< typename T1::elem_type > &r_eigvec, const Base< typename T1::elem_type, T1 > &X, const typename arma_blas_type_only< typename T1::elem_type >::result *junk=0)
Eigenvalues and eigenvectors (both left and right) of general real/complex square matrix X...
Definition: fn_eig.hpp:124
void reset()
Definition: Mat_meat.hpp:4553
void arma_cold arma_stop(const T1 &x)
print a message to get_stream_err1() and/or throw a logic_error exception
Definition: debug.hpp:98
#define arma_extra_debug_sigprint
Definition: debug.hpp:1116
void arma_cold arma_bad(const T1 &x, const bool hurl=true)
print a message to get_stream_err2() and/or throw a run-time error exception
Definition: debug.hpp:150
Dense matrix class.
bool eig_sym(Col< typename T1::pod_type > &eigval, const Base< typename T1::elem_type, T1 > &X, const typename arma_blas_type_only< typename T1::elem_type >::result *junk=0)
Eigenvalues of real/complex symmetric/hermitian matrix X.
Definition: fn_eig.hpp:30


armadillo_matrix
Author(s):
autogenerated on Fri Apr 16 2021 02:31:57