1 // This is core/vnl/algo/vnl_complex_eigensystem.cxx
2 #include <iostream>
3 #include <cassert>
4 #include "vnl_complex_eigensystem.h"
5 // \author fsm
6
7 #include <vnl/vnl_matlab_print.h>
8 #include <vnl/vnl_complexify.h>
9 #include <vnl/algo/vnl_netlib.h> // zgeev_()
10
compute(vnl_matrix<std::complex<double>> const & A,bool right,bool left)11 void vnl_complex_eigensystem::compute(vnl_matrix<std::complex<double> > const & A,
12 bool right,
13 bool left)
14 {
15 A.assert_size(N, N);
16
17 A.assert_finite();
18 assert(! A.is_zero());
19
20 if (right)
21 R.set_size(N, N);
22 if (left)
23 L.set_size(N, N);
24
25 //
26 // Remember that fortran matrices and C matrices are transposed
27 // relative to each other. Moreover, the documentation for zgeev
28 // says that left eigenvectors u satisfy u^h A = lambda u^h,
29 // where ^h denotes adjoint (conjugate transpose).
30 // So we pass our left eigenvector storage as their right
31 // eigenvector storage and vice versa.
32 // But then we also have to conjugate our R after calling the routine.
33 //
34 vnl_matrix<std::complex<double> > tmp(A);
35
36 long work_space=10*N;
37 vnl_vector<std::complex<double> > work(work_space);
38
39 long rwork_space=2*N;
40 vnl_vector<double> rwork(rwork_space);
41
42 long info;
43 long tmpN = N;
44 v3p_netlib_zgeev_(
45 right ? "V" : "N", // jobvl
46 left ? "V" : "N", // jobvr
47 &tmpN, // n
48 tmp.data_block(), // a
49 &tmpN, // lda
50 W.data_block(), // w
51 right ? R.data_block() : nullptr, // vl
52 &tmpN, // ldvl
53 left ? L.data_block() : nullptr, // vr
54 &tmpN, // ldvr
55 work.data_block(), // work
56 &work_space, // lwork
57 rwork.data_block(), // rwork
58 &info, // info
59 1, 1);
60 assert(tmpN == int(N));
61
62 if (right) {
63 // conjugate all elements of R :
64 for (unsigned int i=0;i<N;i++)
65 for (unsigned int j=0;j<N;j++)
66 R(i,j) = std::conj( R(i,j) );
67 }
68
69 if (info == 0) {
70 // success
71 }
72 else if (info < 0) {
73 std::cerr << __FILE__ ": info = " << info << std::endl
74 << __FILE__ ": " << (-info) << "th argument has illegal value\n";
75 assert(false);
76 }
77 else /* if (info > 0) */ {
78 std::cerr << __FILE__ ": info = " << info << std::endl
79 << __FILE__ ": QR algorithm failed to compute all eigenvalues.\n";
80 vnl_matlab_print(std::cerr, A, "A", vnl_matlab_print_format_long);
81 assert(false);
82 }
83 }
84
85 //--------------------------------------------------------------------------------
86
87 //
vnl_complex_eigensystem(vnl_matrix<std::complex<double>> const & A,bool right,bool left)88 vnl_complex_eigensystem::vnl_complex_eigensystem(vnl_matrix<std::complex<double> > const &A,
89 bool right,
90 bool left)
91 : N(A.rows())
92 // L and R are intentionally not initialized.
93 , W(N)
94 {
95 compute(A, right, left);
96 }
97
98 //
vnl_complex_eigensystem(vnl_matrix<double> const & A_real,vnl_matrix<double> const & A_imag,bool right,bool left)99 vnl_complex_eigensystem::vnl_complex_eigensystem(vnl_matrix<double> const &A_real,
100 vnl_matrix<double> const &A_imag,
101 bool right,
102 bool left)
103 : N(A_real.rows())
104 // L and R are intentionally not initialized.
105 , W(N)
106 {
107 A_real.assert_size(N,N);
108 A_imag.assert_size(N,N);
109
110 vnl_matrix<std::complex<double> > A(N,N);
111 vnl_complexify(A_real.begin(), A_imag.begin(), A.begin(), A.size());
112
113 compute(A, right, left);
114 }
115