Double_t eps = 1.e-12;

cout << "Perform the fit  y = c0 + c1 * x in four different ways" << endl;

const Int_t nrVar  = 2;
const Int_t nrPnts = 4;

Double_t ax[] = {0.0,1.0,2.0,3.0};
Double_t ay[] = {1.4,1.5,3.7,4.1};
Double_t ae[] = {0.5,0.2,1.0,0.5};

TVectorD x; x.Use(nrPnts,ax);
TVectorD y; y.Use(nrPnts,ay);
TVectorD e; e.Use(nrPnts,ae);

TMatrixD A(nrPnts,nrVar);
TMatrixDColumn(A,0) = 1.0;
TMatrixDColumn(A,1) = x;

cout << " - 1. solve through Normal Equations" << endl;

const TVectorD c_norm = NormalEqn(A,y,e);

cout << " - 2. solve through SVD" << endl;

TMatrixD Aw = A;
TVectorD yw = y;
for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
   TMatrixDRow(Aw,irow) *= 1/e(irow);
   yw(irow) /= e(irow);
}

TDecompSVD svd(Aw);
Bool_t ok;
const TVectorD c_svd = svd.Solve(yw,ok);

cout << " - 3. solve with pseudo inverse" << endl;

const TMatrixD pseudo1  = svd.Invert();
TVectorD c_pseudo1 = yw;
c_pseudo1 *= pseudo1;

cout << " - 4. solve with pseudo inverse, calculated brute force" << endl;

TMatrixDSym AtA(TMatrixDSym::kAtA,Aw);
const TMatrixD pseudo2 = AtA.Invert() * Aw.T();
TVectorD c_pseudo2 = yw;
c_pseudo2 *= pseudo2;

cout << " - 5. Minuit through TGraph" << endl;

TGraphErrors *gr = new TGraphErrors(nrPnts,ax,ay,0,ae);
TF1 *f1 = new TF1("f1","pol1",0,5);
gr->Fit("f1","Q");
TVectorD c_graph(nrVar);
c_graph(0) = f1->GetParameter(0);
c_graph(1) = f1->GetParameter(1);

Bool_t same = kTRUE;
same &= VerifyVectorIdentity(c_norm,c_svd,0,eps);
same &= VerifyVectorIdentity(c_norm,c_pseudo1,0,eps);
same &= VerifyVectorIdentity(c_norm,c_pseudo2,0,eps);
same &= VerifyVectorIdentity(c_norm,c_graph,0,eps);
if (same)
   cout << " All solutions are the same within tolerance of " << eps << endl;
else
   cout << " Some solutions differ more than the allowed tolerance of " << eps << endl;