vnl_generalized_schur.cxx

InsightToolkit-1.4.0(有大量的优化算法程序)

// This is vxl/vnl/algo/vnl_generalized_schur.cxx
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma implementation
#endif
//:
// \file
// \author fsm

#include "vnl_generalized_schur.h"

#include <vcl_iostream.h>
#include <vcl_cassert.h>

#include <vnl/vnl_vector.h>

#include "vnl_netlib.h" // dgges_()

VCL_DEFINE_SPECIALIZATION
bool 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)
{
  assert(A->cols() == A->cols());
  assert(A->cols() == B->rows());
  assert(A->cols() == B->cols());

  int n = A->rows();
  assert(alphar!=0); alphar->resize(n);    alphar->fill(0);
  assert(alphai!=0); alphai->resize(n);    alphai->fill(0);
  assert(beta!=0);   beta  ->resize(n);    beta  ->fill(0);
  assert(L!=0);      L     ->resize(n, n); L     ->fill(0);
  assert(R!=0);      R     ->resize(n, n); R     ->fill(0);

  int sdim = 0;
  int lwork = 1000 + (8*n + 16);
  double *work = new double[lwork];
  int info = 0;
  A->inplace_transpose();
  B->inplace_transpose();
  dgges_ ("V", "V",
          "N",
          0,
          &n,
          A->data_block(), &n,
          B->data_block(), &n,
          &sdim,
          alphar->data_block(),
          alphai->data_block(),
          beta->data_block(),
          L->data_block(), &n,
          R->data_block(), &n,
          &work[0], &lwork,
          0,
          &info);
  A->inplace_transpose();
  B->inplace_transpose();
  L->inplace_transpose();
  R->inplace_transpose();
  delete [] work;

  if (info == 0) {
    // ok
    return true;
  }
  else {
    // These return codes are taken from dgges.f:
    //*          = 0:  successful exit
    //*          < 0:  if INFO = -i, the i-th argument had an illegal value.
    //*          = 1,...,N:
    //*                The QZ iteration failed.  (A,B) are not in Schur
    //*                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
    //*                be correct for j=INFO+1,...,N.
    //*          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
    //*                =N+2: after reordering, roundoff changed values of
    //*                      some complex eigenvalues so that leading
    //*                      eigenvalues in the Generalized Schur form no
    //*                      longer satisfy DELZTG=.TRUE.  This could also
    //*                      be caused due to scaling.
    //*                =N+3: reordering failed in DTGSEN.
    vcl_cerr << __FILE__ ": info = " << info << ", something went wrong:\n";
    if (info < 0) {
      vcl_cerr << __FILE__ ": (internal error) the " << (-info) << "th argument had an illegal value\n";
    }
    else if (1 <= info && info <= n) {
      vcl_cerr << __FILE__ ": the QZ iteration failed, but the last " << (n - info) << " eigenvalues may be correct\n";
    }
    else if (info == n+1) {
      vcl_cerr << __FILE__ ": something went wrong in DHGEQZ\n";
    }
    else if (info == n+2) {
      vcl_cerr << __FILE__ ": roundoff error -- maybe due to poor scaling\n";
    }
    else if (info == n+3) {
      vcl_cerr << __FILE__ ": reordering failed in DTGSEN\n";
    }
    else {
      vcl_cerr << __FILE__ ": unknown error\n";
    }
    return false;
  }
}