Int_t msize=6;

if (msize < 2 || msize > 10) {
   std::cout << "2 <= msize <= 10" <<std::endl;
   return;
}
std::cout << "--------------------------------------------------------" <<std::endl;
std::cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<std::endl;
std::cout << "For each inversion procedure we check the maximum size  " <<std::endl;
std::cout << "of the off-diagonal elements of Inv(A) * A              " <<std::endl;
std::cout << "--------------------------------------------------------" <<std::endl;

TMatrixD H_square = THilbertMatrixD(msize,msize);

std::cout << "1. Use .InvertFast(&det)" <<std::endl;
if (msize > 6)
   std::cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<std::endl;

Double_t det1;
TMatrixD H1 = H_square;
H1.InvertFast(&det1);

TMatrixD U1(H1,TMatrixD::kMult,H_square);
TMatrixDDiag diag1(U1); diag1 = 0.0;
const Double_t U1_max_offdiag = (U1.Abs()).Max();
std::cout << "  Maximum off-diagonal = " << U1_max_offdiag << std::endl;
std::cout << "  Determinant          = " << det1 << std::endl;

std::cout << "2. Use .Invert(&det)" << std::endl;

Double_t det2;
TMatrixD H2 = H_square;
H2.Invert(&det2);

TMatrixD U2(H2,TMatrixD::kMult,H_square);
TMatrixDDiag diag2(U2); diag2 = 0.0;
const Double_t U2_max_offdiag = (U2.Abs()).Max();
std::cout << "  Maximum off-diagonal = " << U2_max_offdiag << std::endl;
std::cout << "  Determinant          = " << det2 << std::endl;

std::cout << "3. Use TDecompLU" << std::endl;

TMatrixD H3 = H_square;
TDecompLU lu(H_square);

lu.Invert(H3);
Double_t d1_lu; Double_t d2_lu;
lu.Det(d1_lu,d2_lu);
Double_t det3 = d1_lu*TMath::Power(2.,d2_lu);

TMatrixD U3(H3,TMatrixD::kMult,H_square);
TMatrixDDiag diag3(U3); diag3 = 0.0;
const Double_t U3_max_offdiag = (U3.Abs()).Max();
std::cout << "  Maximum off-diagonal = " << U3_max_offdiag << std::endl;
std::cout << "  Determinant          = " << det3 << std::endl;

std::cout << "4. Use TDecompSVD on non-square matrix" << std::endl;

TMatrixD H_nsquare = THilbertMatrixD(msize,msize-1);

TDecompSVD svd(H_nsquare);

TMatrixD H4 = svd.Invert();
Double_t d1_svd; Double_t d2_svd;
svd.Det(d1_svd,d2_svd);
Double_t det4 = d1_svd*TMath::Power(2.,d2_svd);

TMatrixD U4(H4,TMatrixD::kMult,H_nsquare);
TMatrixDDiag diag4(U4); diag4 = 0.0;
const Double_t U4_max_offdiag = (U4.Abs()).Max();
std::cout << "  Maximum off-diagonal = " << U4_max_offdiag << std::endl;
std::cout << "  Determinant          = " << det4 << std::endl;