core/vnl/algo/vnl_generalized_schur.cxx
Go to the documentation of this file.
00001 // This is core/vnl/algo/vnl_generalized_schur.cxx
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005 //:
00006 // \file
00007 // \author fsm
00008 
00009 #include "vnl_generalized_schur.h"
00010 
00011 #include <vcl_iostream.h>
00012 #include <vcl_cassert.h>
00013 
00014 #include <vnl/vnl_vector.h>
00015 
00016 #include <vnl/algo/vnl_netlib.h> // dgges_()
00017 
00018 VCL_DEFINE_SPECIALIZATION
00019 bool vnl_generalized_schur(vnl_matrix<double> *A,
00020                            vnl_matrix<double> *B,
00021                            vnl_vector<double> *alphar,
00022                            vnl_vector<double> *alphai,
00023                            vnl_vector<double> *beta,
00024                            vnl_matrix<double> *L,
00025                            vnl_matrix<double> *R)
00026 {
00027   assert(A->cols() == A->cols());
00028   assert(A->cols() == B->rows());
00029   assert(A->cols() == B->cols());
00030 
00031   long n = A->rows();
00032   assert(alphar!=0); alphar->set_size(n);    alphar->fill(0);
00033   assert(alphai!=0); alphai->set_size(n);    alphai->fill(0);
00034   assert(beta!=0);   beta  ->set_size(n);    beta  ->fill(0);
00035   assert(L!=0);      L     ->set_size(n, n); L     ->fill(0);
00036   assert(R!=0);      R     ->set_size(n, n); R     ->fill(0);
00037 
00038   long sdim = 0;
00039   long lwork = 1000 + (8*n + 16);
00040   double *work = new double[lwork];
00041   long info = 0;
00042   A->inplace_transpose();
00043   B->inplace_transpose();
00044   v3p_netlib_dgges_ ("V", "V",
00045                      "N",
00046                      0,
00047                      &n,
00048                      A->data_block(), &n,
00049                      B->data_block(), &n,
00050                      &sdim,
00051                      alphar->data_block(),
00052                      alphai->data_block(),
00053                      beta->data_block(),
00054                      L->data_block(), &n,
00055                      R->data_block(), &n,
00056                      &work[0], &lwork,
00057                      0,
00058                      &info, 1, 1, 1);
00059   A->inplace_transpose();
00060   B->inplace_transpose();
00061   L->inplace_transpose();
00062   R->inplace_transpose();
00063   delete [] work;
00064 
00065   if (info == 0) {
00066     // ok
00067     return true;
00068   }
00069   else
00070   {
00071     // These return codes are taken from dgges.f:
00072     //*          = 0:  successful exit
00073     //*          < 0:  if INFO = -i, the i-th argument had an illegal value.
00074     //*          = 1,...,N:
00075     //*                The QZ iteration failed.  (A,B) are not in Schur
00076     //*                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
00077     //*                be correct for j=INFO+1,...,N.
00078     //*          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
00079     //*                =N+2: after reordering, roundoff changed values of
00080     //*                      some complex eigenvalues so that leading
00081     //*                      eigenvalues in the Generalized Schur form no
00082     //*                      longer satisfy DELZTG=.TRUE.  This could also
00083     //*                      be caused due to scaling.
00084     //*                =N+3: reordering failed in DTGSEN.
00085     vcl_cerr << __FILE__ ": info = " << info << ", something went wrong:\n";
00086     if (info < 0) {
00087       vcl_cerr << __FILE__ ": (internal error) the " << (-info) << "th argument had an illegal value\n";
00088     }
00089     else if (1 <= info && info <= n) {
00090       vcl_cerr << __FILE__ ": the QZ iteration failed, but the last " << (n - info) << " eigenvalues may be correct\n";
00091     }
00092     else if (info == n+1) {
00093       vcl_cerr << __FILE__ ": something went wrong in DHGEQZ\n";
00094     }
00095     else if (info == n+2) {
00096       vcl_cerr << __FILE__ ": roundoff error -- maybe due to poor scaling\n";
00097     }
00098     else if (info == n+3) {
00099       vcl_cerr << __FILE__ ": reordering failed in DTGSEN\n";
00100     }
00101     else {
00102       vcl_cerr << __FILE__ ": unknown error\n";
00103     }
00104     return false;
00105   }
00106 }