#include <iostream>
#include "TDecompSVD.h"
#include "TMath.h"

void testTDecompSVD()
{

  ////////////////////////////////////////////////
  // use a nice simple diagonal matrix for testing
  //
  //               25           0           0
  // m    =         0          25           0
  //                0           0          25
  //
  // det(m) = 15625
  //
  //             0.04           0           0
  // m^-1 =         0        0.04           0
  //                0           0        0.04
  ////////////////////////////////////////////////

  TMatrixD m(3,3);
  m.Zero();
  m(0,0) = 25.;
  m(1,1) = 25.;
  m(2,2) = 25.;

  m.Print();

  TDecompSVD svd(m);

  ////////////////////////////////////////////////
  // STEP 1: Decompose...
  ////////////////////////////////////////////////
  cout << "svd.Decompose() = " << svd.Decompose() << endl;

  ////////////////////////////////////////////////
  // STEP 2: Calculate determinant...
  ////////////////////////////////////////////////
  Double_t d1, d2;
  svd.Det(d1, d2);
  Double_t det = d1*TMath::Power(2., d2);
  cout << "Determinant     = " << det << endl;

  ////////////////////////////////////////////////
  // STEP 3: Get condition number...
  ////////////////////////////////////////////////
  cout << "svd.Condition() = " << svd.Condition() << endl;

  ////////////////////////////////////////////////
  // STEP 4: Invert...
  ////////////////////////////////////////////////
  TMatrixD mInv = svd.Invert();
  mInv.Print();

}