// @(#)root/matrix:$Name:  $:$Id: TMatrixDSymEigen.cxx,v 1.9 2005/02/15 16:17:10 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann  Dec 2003

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TMatrixDSymEigen                                                     //
//                                                                      //
// Eigenvalues and eigenvectors of a real symmetric matrix.             //
//                                                                      //
// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is  //
// diagonal and the eigenvector matrix V is orthogonal. That is, the    //
// diagonal values of D are the eigenvalues, and V*V' = I, where I is   //
// the identity matrix.  The columns of V represent the eigenvectors in //
// the sense that A*V = V*D.                                            //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "TMatrixDSymEigen.h"

ClassImp(TMatrixDSymEigen)

//______________________________________________________________________________
 TMatrixDSymEigen::TMatrixDSymEigen(const TMatrixDSym &a)
{
  Assert(a.IsValid());

  const Int_t nRows  = a.GetNrows();
  const Int_t rowLwb = a.GetRowLwb();

  fEigenValues.ResizeTo(rowLwb,rowLwb+nRows-1);
  fEigenVectors.ResizeTo(a);

  fEigenVectors = a;

  TVectorD offDiag;
  Double_t work[kWorkMax];
  if (nRows > kWorkMax) offDiag.ResizeTo(nRows);
  else                  offDiag.Use(nRows,work);

  // Tridiagonalize matrix
  MakeTridiagonal(fEigenVectors,fEigenValues,offDiag);

  // Make eigenvectors and -values
  MakeEigenVectors(fEigenVectors,fEigenValues,offDiag);
}

//______________________________________________________________________________
 TMatrixDSymEigen::TMatrixDSymEigen(const TMatrixDSymEigen &another)
{
  *this = another;
}

//______________________________________________________________________________
 void TMatrixDSymEigen::MakeTridiagonal(TMatrixD &v,TVectorD &d,TVectorD &e)
{
// This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and
// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.

  Double_t *pV = v.GetMatrixArray();
  Double_t *pD = d.GetMatrixArray();
  Double_t *pE = e.GetMatrixArray();

  const Int_t n = v.GetNrows();

  Int_t i,j,k;
  Int_t off_n1 = (n-1)*n;
  for (j = 0; j < n; j++)
    pD[j] = pV[off_n1+j];

  // Householder reduction to tridiagonal form.

  for (i = n-1; i > 0; i--) {
    const Int_t off_i1 = (i-1)*n;
    const Int_t off_i  = i*n;

    // Scale to avoid under/overflow.

    Double_t scale = 0.0;
    Double_t h = 0.0;
    for (k = 0; k < i; k++)
      scale = scale+TMath::Abs(pD[k]);
    if (scale == 0.0) {
      pE[i] = pD[i-1];
      for (j = 0; j < i; j++) {
        const Int_t off_j = j*n;
        pD[j] = pV[off_i1+j];
        pV[off_i+j] = 0.0;
        pV[off_j+i] = 0.0;
      }
    } else {

     // Generate Householder vector.

      for (k = 0; k < i; k++) {
        pD[k] /= scale;
        h += pD[k]*pD[k];
      }
      Double_t f = pD[i-1];
      Double_t g = TMath::Sqrt(h);
      if (f > 0)
        g = -g;
      pE[i]   = scale*g;
      h       = h-f*g;
      pD[i-1] = f-g;
      for (j = 0; j < i; j++)
        pE[j] = 0.0;

      // Apply similarity transformation to remaining columns.

      for (j = 0; j < i; j++) {
        const Int_t off_j = j*n;
        f = pD[j];
        pV[off_j+i] = f;
        g = pE[j]+pV[off_j+j]*f;
        for (k = j+1; k <= i-1; k++) {
          const Int_t off_k = k*n;
          g += pV[off_k+j]*pD[k];
          pE[k] += pV[off_k+j]*f;
        }
        pE[j] = g;
      }
      f = 0.0;
      for (j = 0; j < i; j++) {
        pE[j] /= h;
        f += pE[j]*pD[j];
      }
      Double_t hh = f/(h+h);
      for (j = 0; j < i; j++)
        pE[j] -= hh*pD[j];
      for (j = 0; j < i; j++) {
        f = pD[j];
        g = pE[j];
        for (k = j; k <= i-1; k++) {
          const Int_t off_k = k*n;
          pV[off_k+j] -= (f*pE[k]+g*pD[k]);
        }
        pD[j] = pV[off_i1+j];
        pV[off_i+j] = 0.0;
      }
    }
    pD[i] = h;
  }

  // Accumulate transformations.

  for (i = 0; i < n-1; i++) {
    const Int_t off_i  = i*n;
    pV[off_n1+i] = pV[off_i+i];
    pV[off_i+i] = 1.0;
    Double_t h = pD[i+1];
    if (h != 0.0) {
      for (k = 0; k <= i; k++) {
        const Int_t off_k = k*n;
        pD[k] = pV[off_k+i+1]/h;
      }
      for (j = 0; j <= i; j++) {
        Double_t g = 0.0;
        for (k = 0; k <= i; k++) {
          const Int_t off_k = k*n;
          g += pV[off_k+i+1]*pV[off_k+j];
        }
        for (k = 0; k <= i; k++) {
          const Int_t off_k = k*n;
          pV[off_k+j] -= g*pD[k];
        }
      }
    }
    for (k = 0; k <= i; k++) {
      const Int_t off_k = k*n;
      pV[off_k+i+1] = 0.0;
    }
  }
  for (j = 0; j < n; j++) {
    pD[j] = pV[off_n1+j];
    pV[off_n1+j] = 0.0;
  }
  pV[off_n1+n-1] = 1.0;
  pE[0] = 0.0;
}

//______________________________________________________________________________
 void TMatrixDSymEigen::MakeEigenVectors(TMatrixD &v,TVectorD &d,TVectorD &e)
{
// Symmetric tridiagonal QL algorithm.
// This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and
// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.

  Double_t *pV = v.GetMatrixArray();
  Double_t *pD = d.GetMatrixArray();
  Double_t *pE = e.GetMatrixArray();

  const Int_t n = v.GetNrows();

  Int_t i,j,k,l;
  for (i = 1; i < n; i++)
    pE[i-1] = pE[i];
  pE[n-1] = 0.0;

  Double_t f = 0.0;
  Double_t tst1 = 0.0;
  Double_t eps = TMath::Power(2.0,-52.0);
  for (l = 0; l < n; l++) {

    // Find small subdiagonal element

    tst1 = TMath::Max(tst1,TMath::Abs(pD[l])+TMath::Abs(pE[l]));
    Int_t m = l;

    // Original while-loop from Java code
    while (m < n) {
      if (TMath::Abs(pE[m]) <= eps*tst1) {
        break;
      }
      m++;
    }

    // If m == l, pD[l] is an eigenvalue,
    // otherwise, iterate.

    if (m > l) {
      Int_t iter = 0;
      do {
        if (iter++ == 30) {  // (check iteration count here.)
          Error("MakeEigenVectors","too many iterations");
          break;
        }

        // Compute implicit shift

        Double_t g = pD[l];
        Double_t p = (pD[l+1]-g)/(2.0*pE[l]);
        Double_t r = TMath::Hypot(p,1.0);
        if (p < 0)
          r = -r;
        pD[l] = pE[l]/(p+r);
        pD[l+1] = pE[l]*(p+r);
        Double_t dl1 = pD[l+1];
        Double_t h = g-pD[l];
        for (i = l+2; i < n; i++)
          pD[i] -= h;
        f = f+h;

        // Implicit QL transformation.

        p = pD[m];
        Double_t c = 1.0;
        Double_t c2 = c;
        Double_t c3 = c;
        Double_t el1 = pE[l+1];
        Double_t s = 0.0;
        Double_t s2 = 0.0;
        for (i = m-1; i >= l; i--) {
          c3 = c2;
          c2 = c;
          s2 = s;
          g = c*pE[i];
          h = c*p;
          r = TMath::Hypot(p,pE[i]);
          pE[i+1] = s*r;
          s = pE[i]/r;
          c = p/r;
          p = c*pD[i]-s*g;
          pD[i+1] = h+s*(c*g+s*pD[i]);

          // Accumulate transformation.

          for (k = 0; k < n; k++) {
            const Int_t off_k = k*n;
            h = pV[off_k+i+1];
            pV[off_k+i+1] = s*pV[off_k+i]+c*h;
            pV[off_k+i]   = c*pV[off_k+i]-s*h;
          }
        }
        p = -s*s2*c3*el1*pE[l]/dl1;
        pE[l] = s*p;
        pD[l] = c*p;

        // Check for convergence.

      } while (TMath::Abs(pE[l]) > eps*tst1);
    }
    pD[l] = pD[l]+f;
    pE[l] = 0.0;
  }

  // Sort eigenvalues and corresponding vectors.

  for (i = 0; i < n-1; i++) {
    Int_t k = i;
    Double_t p = pD[i];
    for (j = i+1; j < n; j++) {
      if (pD[j] > p) {
        k = j;
        p = pD[j];
      }
    }
    if (k != i) {
      pD[k] = pD[i];
      pD[i] = p;
      for (j = 0; j < n; j++) {
        const Int_t off_j = j*n;
        p = pV[off_j+i];
        pV[off_j+i] = pV[off_j+k];
        pV[off_j+k] = p;
      }
    }
  }
}

//______________________________________________________________________________
TMatrixDSymEigen &TMatrixDSymEigen::operator=(const TMatrixDSymEigen &source)
{
  if (this != &source) {
    fEigenVectors.ResizeTo(source.fEigenVectors);
    fEigenValues.ResizeTo(source.fEigenValues);
  }
  return *this;
}