Go to the documentation of this file.00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006
00007
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>
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
00067 return true;
00068 }
00069 else
00070 {
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
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 }