1 // This is core/vnl/algo/vnl_generalized_schur.cxx
2 //:
3 // \file
4 // \author fsm
5
6 #include <iostream>
7 #include <cassert>
8 #include "vnl_generalized_schur.h"
9
10 #include <vnl/algo/vnl_netlib.h> // dgges_()
11
12 template <>
vnl_generalized_schur(vnl_matrix<double> * A,vnl_matrix<double> * B,vnl_vector<double> * alphar,vnl_vector<double> * alphai,vnl_vector<double> * beta,vnl_matrix<double> * L,vnl_matrix<double> * R)13 bool vnl_generalized_schur(vnl_matrix<double> *A,
14 vnl_matrix<double> *B,
15 vnl_vector<double> *alphar,
16 vnl_vector<double> *alphai,
17 vnl_vector<double> *beta,
18 vnl_matrix<double> *L,
19 vnl_matrix<double> *R)
20 {
21 assert(A->cols() == A->cols());
22 assert(A->cols() == B->rows());
23 assert(A->cols() == B->cols());
24
25 long n = A->rows();
26 assert(alphar!=nullptr); alphar->set_size(n); alphar->fill(0);
27 assert(alphai!=nullptr); alphai->set_size(n); alphai->fill(0);
28 assert(beta!=nullptr); beta ->set_size(n); beta ->fill(0);
29 assert(L!=nullptr); L ->set_size(n, n); L ->fill(0);
30 assert(R!=nullptr); R ->set_size(n, n); R ->fill(0);
31
32 long sdim = 0;
33 long lwork = 1000 + (8*n + 16);
34 auto *work = new double[lwork];
35 long info = 0;
36 A->inplace_transpose();
37 B->inplace_transpose();
38 v3p_netlib_dgges_ ("V", "V",
39 "N",
40 nullptr,
41 &n,
42 A->data_block(), &n,
43 B->data_block(), &n,
44 &sdim,
45 alphar->data_block(),
46 alphai->data_block(),
47 beta->data_block(),
48 L->data_block(), &n,
49 R->data_block(), &n,
50 &work[0], &lwork,
51 nullptr,
52 &info, 1, 1, 1);
53 A->inplace_transpose();
54 B->inplace_transpose();
55 L->inplace_transpose();
56 R->inplace_transpose();
57 delete [] work;
58
59 if (info == 0) {
60 // ok
61 return true;
62 }
63 else
64 {
65 // These return codes are taken from dgges.f:
66 //* = 0: successful exit
67 //* < 0: if INFO = -i, the i-th argument had an illegal value.
68 //* = 1,...,N:
69 //* The QZ iteration failed. (A,B) are not in Schur
70 //* form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
71 //* be correct for j=INFO+1,...,N.
72 //* > N: =N+1: other than QZ iteration failed in DHGEQZ.
73 //* =N+2: after reordering, roundoff changed values of
74 //* some complex eigenvalues so that leading
75 //* eigenvalues in the Generalized Schur form no
76 //* longer satisfy DELZTG=.TRUE. This could also
77 //* be caused due to scaling.
78 //* =N+3: reordering failed in DTGSEN.
79 std::cerr << __FILE__ ": info = " << info << ", something went wrong:\n";
80 if (info < 0) {
81 std::cerr << __FILE__ ": (internal error) the " << (-info) << "th argument had an illegal value\n";
82 }
83 else if (1 <= info && info <= n) {
84 std::cerr << __FILE__ ": the QZ iteration failed, but the last " << (n - info) << " eigenvalues may be correct\n";
85 }
86 else if (info == n+1) {
87 std::cerr << __FILE__ ": something went wrong in DHGEQZ\n";
88 }
89 else if (info == n+2) {
90 std::cerr << __FILE__ ": roundoff error -- maybe due to poor scaling\n";
91 }
92 else if (info == n+3) {
93 std::cerr << __FILE__ ": reordering failed in DTGSEN\n";
94 }
95 else {
96 std::cerr << __FILE__ ": unknown error\n";
97 }
98 return false;
99 }
100 }
101