{
TMatrixD m(3,3);
TMatrixD n(3,3);

m[0][0] = -2600.61;
m[0][1] = -3009.4;
m[0][2] = -33.82;

m[1][0] = 1068.99;
m[1][1] = -2125.77;
m[1][2] = -72.79;

m[2][0] = -38087.3;
m[2][1] = 5404.86;
m[2][2] = -150.3;

n = m;
Double_t det;
m.Invert(&det);

m.Print();
n.Print();
(m*n).Print();
(n*m).Print();


// 4x4 matrix
Double_t x1=5351;
Double_t x2=506;
Double_t y1=419;
Double_t y2=5125;

TMatrixD A (4,4);
TMatrixD B (4,4);
TMatrixD C (4,4);   
   
A[0][0]=x1*x1;
A[0][1]=x1;
A[0][2]=y1*y1;
A[0][3]=y1;

A[1][0]=x2*x2;
A[1][1]=x2;
A[1][2]=y2*y2;
A[1][3]=y2;

A[2][0]=x1*x1;
A[2][1]=x1;
A[2][2]=-y2*y2;
A[2][3]=-y2;

A[3][0]=x2*x2;
A[3][1]=x2;
A[3][2]=-y1*y1;
A[3][3]=-y1;


B = A;

// Standard Inversion, Determinant is using a LU decomposition
const Double_t det = A.Determinant();
cout << "det: " << det <<endl;

if (det != 0.0) {
  A.Invert();
  A.Print();
  B.Print();
  (A*B).Print();
  (B*A).Print();
} else {
  // In case of a singular/near singular matrix SVD is the 
  // right decomposition choice

  cout << "A is singular or near singular" <<endl;
  TDecompSVD svd(A);
  const TVectorD singularVal = svd.GetSig();

  // We can set the the threshold for the singular values with SetTol
  // Those smaller than tol*max(singularVal) are set to 0 for
  // Inversion, Determinant, equation solving
  svd.SetTol(1.0e-20);
  singularVal.Print();

  // Let's calculate the matrix condition number, large numbers
  // number indicate  trouble. In case of SVD it is simply max(singular)/min(singular)
  cout << "condition: " <<  svd.Condition() << endl;

  svd.Invert(A);
  A.Print();
  B.Print();
  (A*B).Print();
  (B*A).Print();
}
}