vnl_real_eigensystem.cxx

Go to the documentation of this file.

// This is core/vnl/algo/vnl_real_eigensystem.cxx
//:
// \file
//
// \author Andrew W. Fitzgibbon, Oxford RRG
// \date   23 Jan 97

//-----------------------------------------------------------------------------

#include <iostream>
#include <cassert>
#include "vnl_real_eigensystem.h"
#include <vnl/vnl_fortran_copy.h>
#include <vnl/algo/vnl_netlib.h> // rg_()

//: Extract eigensystem of non-symmetric matrix M, using the EISPACK routine rg.
//  Should probably switch to using LAPACK's dgeev to avoid transposing.
vnl_real_eigensystem::vnl_real_eigensystem(vnl_matrix<double> const & M):
  Vreal(M.rows(), M.columns()),
  V(M.rows(), M.columns()),
  D(M.rows())
{
  long n = M.rows();
  assert(n == (int)(M.columns()));

  vnl_fortran_copy<double> mf(M);

  vnl_vector<double> wr(n);
  vnl_vector<double> wi(n);
  vnl_vector<long> iv1(n);
  vnl_vector<double> fv1(n);
  vnl_matrix<double> devout(n, n);

  long ierr = 0;
  long matz = 1;
  v3p_netlib_rg_(&n, &n, mf,
                 wr.data_block(), wi.data_block(),
                 &matz, devout.data_block(),
                 iv1.data_block(), fv1.data_block(),
                 &ierr);

  if (ierr != 0) {
    std::cerr << " *** vnl_real_eigensystem: Failed on " << ierr << "th eigenvalue\n"
             << M << std::endl;
  }

  // Copy out eigenvalues and eigenvectors
  for (int c = 0; c < n; ++c) {
    D(c,c) = std::complex<double>(wr[c], wi[c]);
    if (wi[c] != 0) {
      // Complex - copy conjugates and inc c.
      D(c+1, c+1) = std::complex<double>(wr[c], -wi[c]);
      for (int r = 0; r < n; ++r) {
        V(r, c) = std::complex<double>(devout(c,r), devout(c+1,r));
        V(r, c+1) = std::complex<double>(devout(c,r), -devout(c+1,r));
      }

      ++c;
    }
    else
      for (int r = 0; r < n; ++r) {
        V(r, c) = std::complex<double>(devout(c,r), 0);
        Vreal(r,c) = devout(c,r);
      }
  }
}