vnl_generalized_schur.cxx

Go to the documentation of this file.

// This is core/vnl/algo/vnl_generalized_schur.cxx
//:
// \file
// \author fsm

#include <iostream>
#include <cassert>
#include "vnl_generalized_schur.h"

#include <vnl/algo/vnl_netlib.h> // dgges_()

template <>
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());

  long n = A->rows();
  assert(alphar!=nullptr); alphar->set_size(n);    alphar->fill(0);
  assert(alphai!=nullptr); alphai->set_size(n);    alphai->fill(0);
  assert(beta!=nullptr);   beta  ->set_size(n);    beta  ->fill(0);
  assert(L!=nullptr);      L     ->set_size(n, n); L     ->fill(0);
  assert(R!=nullptr);      R     ->set_size(n, n); R     ->fill(0);

  long sdim = 0;
  long lwork = 1000 + (8*n + 16);
  auto *work = new double[lwork];
  long info = 0;
  A->inplace_transpose();
  B->inplace_transpose();
  v3p_netlib_dgges_ ("V", "V",
                     "N",
                     nullptr,
                     &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,
                     nullptr,
                     &info, 1, 1, 1);
  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.
    std::cerr << __FILE__ ": info = " << info << ", something went wrong:\n";
    if (info < 0) {
      std::cerr << __FILE__ ": (internal error) the " << (-info) << "th argument had an illegal value\n";
    }
    else if (1 <= info && info <= n) {
      std::cerr << __FILE__ ": the QZ iteration failed, but the last " << (n - info) << " eigenvalues may be correct\n";
    }
    else if (info == n+1) {
      std::cerr << __FILE__ ": something went wrong in DHGEQZ\n";
    }
    else if (info == n+2) {
      std::cerr << __FILE__ ": roundoff error -- maybe due to poor scaling\n";
    }
    else if (info == n+3) {
      std::cerr << __FILE__ ": reordering failed in DTGSEN\n";
    }
    else {
      std::cerr << __FILE__ ": unknown error\n";
    }
    return false;
  }
}