#include <iostream>
#include <cmath>
#include <cfloat>
#include "Math/SMatrix.h"
#include "Math/SVector.h"

using CovSMatrixD = ROOT::Math::SMatrix<double,3, 3, ROOT::Math::MatRepSym<double, 3>>;

#define _FAILING_

void mtst()
{
  std::cout << "initial symmetric matrix" << std::endl;
  CovSMatrixD m0;
  for (int i=3;i--;) {
    m0(i,i) = 1.+i;
    for (int j=i;j--;) {
      m0(i,j) = 1./((i+3)*(j+3));
    }
  }
  std::cout << m0 << std::endl;

  std::cout << "its inverse" << std::endl;
  int status = 0;
  auto m0inv = m0.Inverse(status);
  std::cout << m0inv << std::endl;

  std::cout << "product with its inverse (auto)" << std::endl;
  auto mmi = m0*m0inv;
  std::cout << mmi << std::endl;  

#ifdef _FAILING_
  std::cout << "product with its inverse as symmetric matrix" << std::endl;
  CovSMatrixD mmi1 = m0*m0inv;
  //CovSMatrixD mmi1 = m0;
  //mmi1 *= m0inv;
  std::cout << mmi1 << std::endl;  
#endif
  
}