void solveLinear()
{
#ifdef __CINT__
  gSystem->Load("libMatrix");
#endif
  
  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};

  // Make the vectors 'Use" the data : they are not copied, the vector data
  // pointer is just set appropriately

  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;

  
  // -------- 1.  Normal equations
  const TVectorD c_norm = NormalEqn(A,y,e);


  // -------- 2. SVD - numerically  preferred method
  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);

  // ------- 3. TLinearFitter obviously...
  TLinearFitter* lf =  new TLinearFitter(1);
  lf->SetFormula("1 ++ x");
  TVectorD params;
  lf->AssignData( nrPnts,nrVar,ax,ay,ae );
  lf->Eval();
  lf->GetParameters(params);
  delete lf;


  cout << " - 1. Normal Equations"
       <<"  par[1] "<<c_norm(0)<<"   par[2] "<<c_norm(1)<<endl; 
  cout << " - 2. SVD             " 
       <<"  par[1] "<<c_svd(0)<<"   par[2] "<<c_svd(1)<< endl;
  cout << " - 3. TLinearFitter   " 
       <<"  par[1] "<<params(0)<<"   par[2] "<< params(1) << endl;

}