// $Id: CovMatrix5Diagonalize.cc,v 1.8 2002/10/22 17:52:36 sachs Exp $
//
// Compute Eigenvalues and optionally Eigenvectors of a 5x5
//	Positive definite symmetric matrix.
//	Return true if the computation was successful.
//	Return false if the matrix could not be diagonalized.
//
// David Sachs, Fermilab - March 2002
//
// Based on LAPACK routine dspev.f

#include "CovMatrices/CovMatrix5.h"
#include "CovMatrixEigCommon.h"

using namespace CovMatrices::LAPACK;

bool CovMatrices::CovMatrix5::diagonalize ( Vector5& eigenvalues, Vector5 eigenvectors[] ) const
{

  // Compute Eigenvalues and optionally Eigenvectors 
  // by using C++ versions of proper LAPACK routines
  //
  // Arguments are:
  //
  //  Vector5& eigenvalues:	OUTPUT eignvalues
  //
  //  Vector5 eigenvectors[5]:	OUTPUT eigenvectors
  //	Supply a null pointer to compute Eigenvectors only
  //
  // Return value is true if computation was successful
  //
  
  double in[15];	// Copy of original matrixC
  double val[5];	// Temporary to hold Eigenvalues
  double vec[25];	// Temporary to hold Eigenvectors
  double e[4], t[4];	// Temporary for tridagonal form
  
  memcpy(in, &M, sizeof(M));	// Copy input matrix

  // Convert Matrix to TriDiagonal form
  dsptrd(5, in, val, e, t);

  // For eigenvalues only call dsteqr
  // For eigenvectors call dopgtr and dsteqr
  if (!eigenvectors)
  {
    if (dsteqr(5, val, e) != 0) return false;
    eigenvalues = val;
    return true;
  }
  else
  {
    dopgtr( 5, in, t, vec);
    if (dsteqr(5, val, e, vec) !=0) return false;

    // store eigenvector and eigenvalues
    eigenvalues = val;
    eigenvectors[0] = vec;
    eigenvectors[1] = vec+5;
    eigenvectors[2] = vec+10;
    eigenvectors[3] = vec+15;
    eigenvectors[4] = vec+20;
    return true;
  }
}