Added AliMillePede2 + related matrix and solver classes
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Jan 2009 14:35:57 +0000 (14:35 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Jan 2009 14:35:57 +0000 (14:35 +0000)
STEER/AliMatrixSparse.cxx [new file with mode: 0644]
STEER/AliMatrixSparse.h [new file with mode: 0644]
STEER/AliMatrixSq.cxx [new file with mode: 0644]
STEER/AliMatrixSq.h [new file with mode: 0644]
STEER/AliMillePede2.cxx [new file with mode: 0644]
STEER/AliMillePede2.h [new file with mode: 0644]
STEER/AliMinResSolve.cxx [new file with mode: 0644]
STEER/AliMinResSolve.h [new file with mode: 0644]
STEER/AliSymMatrix.cxx [new file with mode: 0644]
STEER/AliSymMatrix.h [new file with mode: 0644]

diff --git a/STEER/AliMatrixSparse.cxx b/STEER/AliMatrixSparse.cxx
new file mode 100644 (file)
index 0000000..88a0a5f
--- /dev/null
@@ -0,0 +1,261 @@
+#include "AliMatrixSparse.h"
+
+//___________________________________________________________
+AliVectorSparse::AliVectorSparse()
+  : fNElems(0),fIndex(0),fElems(0) {}
+
+//___________________________________________________________
+AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
+  : fNElems(src.fNElems),fIndex(0),fElems(0)
+{
+  fIndex = new UShort_t[fNElems];
+  fElems = new Double_t[fNElems];
+  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
+  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
+}
+
+//___________________________________________________________
+void AliVectorSparse::Clear()
+{
+  delete[] fIndex; fIndex = 0;
+  delete[] fElems; fElems = 0;
+  fNElems = 0;
+}
+
+//___________________________________________________________
+Float_t AliMatrixSparse::GetDensity() const
+{
+  // get fraction of non-zero elements
+  Int_t nel = 0;
+  for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
+  int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
+  return float(nel)/den;
+}
+
+//___________________________________________________________
+AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
+{
+  if (&src==this) return *this;
+  Clear();
+  fNElems = src.fNElems;
+  fIndex = new UShort_t[fNElems];
+  fElems = new Double_t[fNElems];
+  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
+  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
+  //
+  return *this;
+}
+
+//___________________________________________________________
+Double_t AliVectorSparse::FindIndex(Int_t ind) const
+{
+  // return an element with given index
+  //printf("V: findindex\n");
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else return fElems[mid];
+  }
+  return 0.0;
+}
+
+//___________________________________________________________
+void AliVectorSparse::Zero(Int_t ind)
+{
+  // set element to 0 if it was already defined
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else {fElems[mid] = 0.; return;}
+  }
+}
+
+//___________________________________________________________
+Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
+{
+  // increment an element with given index
+  //printf("V: findindexAdd\n");
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else return fElems[mid];
+  }
+  // need to insert a new element
+  UShort_t *arrI = new UShort_t[fNElems+1];
+  memcpy(arrI,fIndex,first*sizeof(UShort_t));
+  arrI[first] = ind;
+  memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
+  delete[] fIndex;
+  fIndex = arrI;
+  //
+  Double_t   *arrE = new Double_t[fNElems+1];
+  memcpy(arrE,fElems,first*sizeof(Double_t));
+  arrE[first] = 0;
+  memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
+  delete[] fElems;
+  fElems = arrE;
+  //
+  fNElems++;
+  return fElems[first];
+  //
+}
+
+//__________________________________________________________
+void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
+{
+  if (sz<1) {Clear(); return;}
+    // need to insert a new element
+  UShort_t *arrI = new UShort_t[sz];
+  Double_t *arrE = new Double_t[sz];
+  memset(arrI,0,sz*sizeof(UShort_t));
+  memset(arrE,0,sz*sizeof(Double_t));
+  //
+  if (copy && fIndex) {
+    int cpsz = TMath::Min(fNElems,sz);
+    memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
+    memcpy(arrE,fElems,cpsz*sizeof(Double_t));
+  }
+  delete[] fIndex;
+  delete[] fElems;
+  fIndex = arrI;
+  fElems = arrE;
+  fNElems = sz;
+  //
+}
+
+//__________________________________________________________
+void AliVectorSparse::SortIndices(Bool_t valuesToo)
+{
+  // sort indices in increasing order. Used to fix the row after ILUk decomposition
+  for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
+       UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
+       if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
+      }
+}
+
+//__________________________________________________________
+void AliVectorSparse::Print(Option_t* )  const
+{
+  printf("|");
+  for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
+  printf("|\n");
+}
+
+
+///////////////////////////////////////////////////////////////////////////////////////////
+//___________________________________________________________
+ClassImp(AliMatrixSparse)
+
+//___________________________________________________________
+AliMatrixSparse::AliMatrixSparse(Int_t sz)
+: AliMatrixSq(),fVecs(0)
+{
+  fNcols=fNrows=sz;
+  //
+  fVecs = new AliVectorSparse*[sz];
+  for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
+}
+
+//___________________________________________________________
+AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
+  : AliMatrixSq(src),fVecs(0)
+{
+  fVecs = new AliVectorSparse*[fNcols];
+  for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
+}
+
+//___________________________________________________________
+AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
+{
+  if (ir>=fNcols) {
+    AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
+    for (int i=GetSize();i--;) arrv[i] = fVecs[i];
+    delete fVecs;
+    fVecs = arrv;    
+    for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
+    fNcols = fNrows = ir+1;
+  }
+  return fVecs[ir];
+}
+
+//___________________________________________________________
+AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
+{
+  if (*this == src) return *this;
+  Clear();
+  fNcols = fNrows = src.GetSize();
+  SetSymmetric(src.IsSymmetric());
+  fVecs = new AliVectorSparse*[fNcols];
+  for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
+  return *this;
+}
+
+//___________________________________________________________
+void AliMatrixSparse::Clear(Option_t*) 
+{
+  for (int i=fNcols;i--;) delete GetRow(i);
+  delete[] fVecs;
+  fNcols = fNrows = 0;
+}
+
+//___________________________________________________________
+void AliMatrixSparse::Print(Option_t*)  const
+{
+  printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":"");
+  for (int i=0;i<fNcols;i++) {
+    AliVectorSparse* row = GetRow(i);
+    if (!row->GetNElems()) continue;
+    printf("%3d: ",i); 
+    row->Print();
+  }
+}
+
+//___________________________________________________________
+void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
+{
+  // fill vecOut by matrix*vecIn
+  // vector should be of the same size as the matrix
+  //
+  memset(vecOut,0,GetSize()*sizeof(Double_t));
+  //
+  for (int rw=GetSize();rw--;) {  // loop over rows >>>
+    const AliVectorSparse* rowV = GetRow(rw);
+    Int_t nel  = rowV->GetNElems();
+    if (!nel) continue;
+    //
+    UShort_t* indV = rowV->GetIndices();
+    Double_t* elmV = rowV->GetElems();
+    //
+    if (IsSymmetric()) {
+      // treat diagonal term separately. If filled, it should be the last one
+      if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel];
+      else nel = rowV->GetNElems(); // diag elem was not filled
+      //
+      for (int iel=nel;iel--;) {          // less element retrieval for symmetric case
+       if (elmV[iel]) {        
+         vecOut[rw]        += vecIn[indV[iel]]*elmV[iel];
+         vecOut[indV[iel]] += vecIn[rw]*elmV[iel];
+       }
+      }
+    }
+    else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
+    //
+  } // loop over rows <<<
+  //
+}
+
+//___________________________________________________________
+void AliMatrixSparse::SortIndices(Bool_t valuesToo)
+{
+  // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
+  for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
+}
diff --git a/STEER/AliMatrixSparse.h b/STEER/AliMatrixSparse.h
new file mode 100644 (file)
index 0000000..6d0e0c1
--- /dev/null
@@ -0,0 +1,144 @@
+#ifndef ALIMATRIXSPARSE_H
+#define ALIMATRIXSPARSE_H
+
+#include "AliMatrixSq.h"
+
+
+///////////////////////////////////////////////////////////////////////////////////////
+class AliVectorSparse {
+ public:
+  AliVectorSparse();
+  AliVectorSparse(const AliVectorSparse& src);
+  virtual ~AliVectorSparse() {Clear();}
+  virtual void Print(Option_t* option="")                 const;
+  //
+  Int_t     GetNElems()                                   const {return fNElems;}
+  UShort_t *GetIndices()                                  const {return fIndex;}
+  Double_t *GetElems()                                    const {return fElems;}
+  UShort_t& GetIndex(Int_t i)                                   {return fIndex[i];}
+  Double_t& GetElem(Int_t i)                              const {return fElems[i];}
+  void      Clear();
+  void      Reset()                                             {memset(fElems,0,fNElems*sizeof(Double_t));}
+  void      ReSize(Int_t sz,Bool_t copy=kFALSE);
+  void      SortIndices(Bool_t valuesToo=kFALSE);
+  //
+  AliVectorSparse& operator=(const AliVectorSparse& src);
+  //
+  virtual Double_t         operator()(Int_t ind)         const;
+  virtual Double_t&        operator()(Int_t ind);
+  virtual void             Zero(Int_t ind);
+  Double_t                 FindIndex(Int_t ind)          const;
+  Double_t&                FindIndexAdd(Int_t ind);
+  //
+  Int_t     GetLastIndex()                               const {return fIndex[fNElems-1];}
+  Double_t  GetLastElem()                                const {return fElems[fNElems-1];}
+  Double_t &GetLastElem()                                      {return fElems[fNElems-1];}
+  //
+ protected:
+  Int_t            fNElems;   // 
+  UShort_t*        fIndex;    // Index of stored elems
+  Double_t*        fElems;    // pointer on elements
+};
+
+
+//___________________________________________________
+inline Double_t AliVectorSparse::operator()(Int_t ind) const
+{
+  return FindIndex(ind);
+}
+
+//___________________________________________________
+inline Double_t& AliVectorSparse::operator()(Int_t ind)
+{
+  return FindIndexAdd(ind);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+
+/////////////////////////////////////////////////////////////////////////////////////////////
+// Sparse matrix class
+//
+class AliMatrixSparse : public AliMatrixSq 
+{
+ public:
+  AliMatrixSparse() : fVecs(0) {}
+  AliMatrixSparse(Int_t size);
+  AliMatrixSparse(const AliMatrixSparse &mat);
+  virtual ~AliMatrixSparse() {Clear();}
+  //
+  AliVectorSparse* GetRow(Int_t ir)        const {return (ir<fNcols) ? fVecs[ir] : 0;}
+  AliVectorSparse* GetRowAdd(Int_t ir);
+
+  void Clear(Option_t* option="");
+  void Reset()                            {for (int i=fNcols;i--;) GetRow(i)->Reset();}
+  void Print(Option_t* option="")                      const;
+  AliMatrixSparse& operator=(const AliMatrixSparse& src);
+  Double_t&        operator()(Int_t row,Int_t col);
+  Double_t         operator()(Int_t row,Int_t col)     const;
+  void             Zero(Int_t row,Int_t col);
+  Float_t          GetDensity()                        const;
+  //
+  Double_t         DiagElem(Int_t r)                   const;
+  Double_t&        DiagElem(Int_t r);
+  void             SortIndices(Bool_t valuesToo=kFALSE);
+  //
+  void MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
+  void MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const {MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());}
+  //
+ protected:
+  //
+  AliVectorSparse** fVecs;
+  //
+  ClassDef(AliMatrixSparse,0)
+};
+
+//___________________________________________________
+inline void AliMatrixSparse::Zero(Int_t row,Int_t col)
+{
+  //  set existing element to 0
+  if (IsSymmetric() && col>row) Swap(row,col); 
+  AliVectorSparse* rowv = GetRow(row);
+  if (rowv) rowv->Zero(col);
+}
+
+//___________________________________________________
+inline Double_t AliMatrixSparse::operator()(Int_t row,Int_t col) const
+{
+  //  printf("M: find\n");
+  if (IsSymmetric() && col>row) Swap(row,col); 
+  AliVectorSparse* rowv = GetRow(row);
+  if (!rowv) return 0;
+  return rowv->FindIndex(col);
+}
+
+//___________________________________________________
+inline Double_t& AliMatrixSparse::operator()(Int_t row,Int_t col)
+{
+  //  printf("M: findindexAdd\n");
+  if (IsSymmetric() && col>row) Swap(row,col); 
+  AliVectorSparse* rowv = GetRowAdd(row);
+  return rowv->FindIndexAdd(col);
+}
+
+//___________________________________________________
+inline Double_t AliMatrixSparse::DiagElem(Int_t row) const
+{
+  AliVectorSparse* rowv = GetRow(row);
+  if (!rowv) return 0;
+  if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ? rowv->GetLastElem() : 0.;
+  else return rowv->FindIndex(row);
+  //
+}
+
+//___________________________________________________
+inline Double_t &AliMatrixSparse::DiagElem(Int_t row)
+{
+  AliVectorSparse* rowv = GetRowAdd(row);
+  if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ? 
+                      rowv->GetLastElem() : rowv->FindIndexAdd(row);
+  else return rowv->FindIndexAdd(row);
+  //
+}
+
+#endif
+
diff --git a/STEER/AliMatrixSq.cxx b/STEER/AliMatrixSq.cxx
new file mode 100644 (file)
index 0000000..d315e54
--- /dev/null
@@ -0,0 +1,43 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream>
+//
+#include "TClass.h"
+#include "TMath.h"
+#include "TVectorD.h"
+#include "AliMatrixSq.h"
+//
+
+using namespace std;
+
+ClassImp(AliMatrixSq)
+
+
+
+//___________________________________________________________
+void AliMatrixSq::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
+{
+  // fill vecOut by matrix*vecIn
+  // vector should be of the same size as the matrix
+  for (int i=GetSize();i--;) {
+    vecOut[i] = 0.0;
+    for (int j=GetSize();j--;) vecOut[i] += vecIn[j]*(*this)(i,j);
+  }
+  //
+}
+
+//___________________________________________________________
+void AliMatrixSq::PrintCOO() const
+{
+  // print matrix in COO sparse format
+  //
+  // get number of non-zero elements
+  int nnz = 0;
+  int sz = GetSize();
+  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if (Querry(ir,ic)!=0) nnz++;
+  //
+  printf("%d %d %d\n",sz,sz,nnz);
+  double vl;
+  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if ((vl=Querry(ir,ic))!=0) printf("%d %d %f\n",ir,ic,vl);
+  //
+}
diff --git a/STEER/AliMatrixSq.h b/STEER/AliMatrixSq.h
new file mode 100644 (file)
index 0000000..09f0828
--- /dev/null
@@ -0,0 +1,72 @@
+#ifndef ALIMATRIXSQ_H
+#define ALIMATRIXSQ_H
+
+#include <TMatrixDBase.h>
+#include <TVectorD.h>
+
+class AliMatrixSq : public TMatrixDBase {
+  //
+ public:
+  AliMatrixSq(): fSymmetric(kFALSE) {}
+  AliMatrixSq(const AliMatrixSq &src) : TMatrixDBase(src), fSymmetric(src.fSymmetric) {}
+  virtual ~AliMatrixSq() {}
+  virtual Int_t   GetSize()                            const {return fNcols;}
+  virtual Float_t GetDensity()                         const     = 0;
+  //
+  virtual void  Clear(Option_t* option="")                       = 0;//{Error("Clear","Dummy");}
+  //
+  virtual       Double_t      Querry(Int_t rown, Int_t coln)     const {return operator()(rown,coln);}
+  virtual       Double_t      operator()(Int_t rown, Int_t coln) const = 0;//{Error("(i,j)","Dummy");return 0;}
+  virtual       Double_t&     operator()(Int_t rown, Int_t coln) = 0;//{Error("(i,j)","Dummy");return 0;}
+  //
+  virtual       Double_t      QuerryDiag(Int_t rc)               const {return DiagElem(rc);}
+  virtual       Double_t      DiagElem(Int_t r)                  const = 0;
+  virtual       Double_t&     DiagElem(Int_t r)                  = 0;
+  //
+  virtual void  Print(Option_t* option="")           const       = 0;//{Error("Print","Dummy");}
+  virtual void  Reset()                                          = 0;
+  virtual void  PrintCOO()                           const;          // print in COO format
+  //
+  virtual void  MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
+  virtual void  MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
+  //
+  Bool_t        IsSymmetric()                       const {return fSymmetric;}
+  void          SetSymmetric(Bool_t v=kTRUE)              {fSymmetric = v;}
+  //
+  // ---------------------------------- Dummy methods of MatrixBase
+  virtual       const Double_t   *GetMatrixArray  () const {Error("GetMatrixArray","Dummy"); return 0;};
+  virtual             Double_t   *GetMatrixArray  ()       {Error("GetMatrixArray","Dummy"); return 0;};
+  virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
+  virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
+  virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
+  virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
+  virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
+  virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
+  virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
+  virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
+  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
+  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
+  //
+  virtual void Allocate      (Int_t ,Int_t ,Int_t , Int_t ,Int_t ,Int_t ) 
+  {Error("Allocate","Dummy"); return;}
+  //
+ protected:
+  //
+  void    Swap(int &r,int &c) const {int t=r;r=c;c=t;}
+  //
+ protected:
+  //
+  Bool_t        fSymmetric;     // is the matrix symmetric? Only lower triangle is filled
+  //
+  ClassDef(AliMatrixSq,1) //Square Matrix Class
+};
+
+
+//___________________________________________________________
+inline void AliMatrixSq::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
+{
+  MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
+}
+
+
+#endif
diff --git a/STEER/AliMillePede2.cxx b/STEER/AliMillePede2.cxx
new file mode 100644 (file)
index 0000000..263d87e
--- /dev/null
@@ -0,0 +1,957 @@
+#include "AliMillePede2.h"
+#include "AliLog.h"
+#include <TStopwatch.h>
+
+using namespace std;
+
+ClassImp(AliMillePedeRecord)
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::AliMillePedeRecord() : 
+fSize(0),fIndex(0),fValue(0) {SetBufferSize(0);}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) : 
+  TObject(src),fSize(src.fSize),fIndex(0),fValue(0)
+{
+  fIndex = new int[GetBufferSize()];
+  memcpy(fIndex,src.fIndex,fSize*sizeof(int));
+  fValue = new double[GetBufferSize()];
+  memcpy(fValue,src.fValue,fSize*sizeof(double));
+}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
+{
+  if (this!=&rhs) {
+    Reset();
+    for (int i=0;i<rhs.GetSize();i++) {
+      double val;
+      int    ind;
+      rhs.GetIndexValue(i,ind,val);
+      AddIndexValue(ind,val);
+    }
+  }
+  return *this;
+}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue;}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::Print(const Option_t *) const
+{
+  if (!fSize) {AliInfo("No data"); return;}
+  int cnt=0,point=0;
+  //  
+  while(cnt<fSize) {
+    //
+    double resid = fValue[cnt++];
+    double *derLoc = GetValue()+cnt;
+    int    *indLoc = GetIndex()+cnt;
+    int     nLoc = 0;
+    while(!IsWeight(cnt)) {nLoc++;cnt++;}
+    double weight = GetValue(cnt++);
+    double *derGlo = GetValue()+cnt;
+    int    *indGlo = GetIndex()+cnt;
+    int     nGlo = 0;
+    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
+    //
+    printf("\n*** Point#%2d | Residual = %+.4e | Weight = %+.4e\n",point++,resid,weight);
+    printf("Locals : "); 
+    for (int i=0;i<nLoc;i++) printf("[%5d] %+.4e|",indLoc[i],derLoc[i]); printf("\n");
+    printf("Globals: "); 
+    for (int i=0;i<nGlo;i++) printf("[%5d] %+.4e|",indGlo[i],derGlo[i]); printf("\n");
+    //
+  }
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::Expand(int bfsize)
+{
+  // add extra space 
+  bfsize = TMath::Max(bfsize,GetBufferSize());
+  int *tmpI = new int[bfsize];
+  memcpy(tmpI,fIndex, fSize*sizeof(int));
+  delete fIndex;
+  fIndex = tmpI;
+  //
+  double *tmpD = new double[bfsize];
+  memcpy(tmpD,fValue, fSize*sizeof(double));
+  delete fValue;
+  fValue = tmpD;
+  //
+  SetBufferSize(bfsize);
+}
+
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+ClassImp(AliMillePede2)
+
+Bool_t   AliMillePede2::fgIsMatGloSparse = kFALSE;   // use faster dense matrix by default
+Int_t    AliMillePede2::fgMinResCondType = kTRUE;        // No preconditioner by default
+Double_t AliMillePede2::fgMinResTol      = 1.e-8;   // default tolerance
+Int_t    AliMillePede2::fgMinResMaxIter  = 3000;     // default max number of iterations 
+Int_t    AliMillePede2::fgIterSol        = AliMinResSolve::kSolMinRes; // default iterative solver
+Int_t    AliMillePede2::fgNKrylovV       = 30;       // default number of Krylov vectors to keep
+
+//_____________________________________________________________________________________________
+AliMillePede2::AliMillePede2() 
+: fNLocPar(0),
+  fNGloPar(0),
+  fNGloSize(0),
+//
+  fNLocEquations(0),
+  fIter(0),
+  fMaxIter(10),
+  fNStdDev(3),
+  fNGloConstraints(0),
+  fNLocFits(0),
+  fNLocFitsRejected(0),
+  fNGloFix(0),
+  fGloSolveStatus(gkFailed),
+//
+  fChi2CutFactor(1.),
+  fChi2CutRef(1.),
+  fResCutInit(100.),
+  fResCut(100.),
+  fMinPntValid(0),
+//
+  fProcPnt(0),
+  fVecBLoc(0),
+  fDiagCGlo(0),
+  fVecBGlo(0),
+  fInitPar(0),
+  fDeltaPar(0),
+  fSigmaPar(0),
+  fIsLinear(0),
+  fConstrUsed(0),
+//
+  fGlo2CGlo(0),
+  fCGlo2Glo(0),
+//
+  fMatCLoc(0),
+  fMatCGlo(0),
+  fMatCGloLoc(0),
+//
+  fDataRecFName("/tmp/mp2_data_records.root"),
+  fRecord(0),
+  fDataRecFile(0),  
+  fTreeData(0),
+  //
+  fConstrRecFName("/tmp/mp2_constraints_records.root"),
+  fTreeConstr(0),
+  fConsRecFile(0),
+  fCurrRecDataID(0),
+  fCurrRecConstrID(0)
+{}
+
+//_____________________________________________________________________________________________
+AliMillePede2::AliMillePede2(const AliMillePede2& src) : 
+  TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0),
+  fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLocFits(0),fNLocFitsRejected(0),
+  fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0),
+  fResCut(0),fMinPntValid(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
+  fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
+  fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fDataRecFName(0),fRecord(0),fDataRecFile(0),
+  fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),fCurrRecConstrID(0)
+{printf("Dummy\n");}
+
+//_____________________________________________________________________________________________
+AliMillePede2::~AliMillePede2() 
+{
+  CloseDataRecStorage();
+  CloseConsRecStorage();
+  //
+  delete[] fProcPnt;
+  delete[] fVecBLoc;
+  delete[] fDiagCGlo;
+  delete[] fVecBGlo;
+  delete[] fInitPar;
+  delete[] fDeltaPar;
+  delete[] fSigmaPar;
+  delete[] fGlo2CGlo;
+  delete[] fCGlo2Glo;
+  delete[] fIsLinear;
+  delete[] fConstrUsed;
+  //
+  delete fRecord;
+  delete fMatCLoc;
+  delete fMatCGlo;
+  delete fMatCGloLoc;
+} 
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit)
+{
+  //
+  if (nLoc>0)        fNLocPar = nLoc;
+  if (nGlo>0)        fNGloPar = nGlo;
+  if (lResCutInit>0) fResCutInit = lResCutInit; 
+  if (lResCut>0)     fResCut     = lResCut;
+  if (lNStdDev>0)    fNStdDev    = lNStdDev;
+  //
+  fNGloSize = fNGloPar;
+  //
+  try {
+    //
+    if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
+    else                   fMatCGlo = new AliSymMatrix(fNGloPar);
+    //
+    fMatCLoc      = new AliSymMatrix(fNLocPar);
+    fMatCGloLoc   = new TMatrixD(fNGloPar,fNLocPar);
+    //
+    fProcPnt      = new Int_t[fNGloPar];
+    fVecBLoc      = new Double_t[fNLocPar];
+    fDiagCGlo     = new Double_t[fNGloPar];
+    //
+    fInitPar      = new Double_t[fNGloPar];
+    fDeltaPar     = new Double_t[fNGloPar];
+    fSigmaPar     = new Double_t[fNGloPar];
+    fIsLinear     = new Bool_t[fNGloPar];
+    //
+    fGlo2CGlo     = new Int_t[fNGloPar];
+    fCGlo2Glo     = new Int_t[fNGloPar];
+  }
+  catch(bad_alloc&) {
+    AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar));
+    return 0;
+  }
+  //
+  memset(fVecBLoc   ,0,fNLocPar*sizeof(Double_t));
+  memset(fDiagCGlo  ,0,fNGloPar*sizeof(Double_t));
+  memset(fInitPar   ,0,fNGloPar*sizeof(Double_t));
+  memset(fDeltaPar  ,0,fNGloPar*sizeof(Double_t));
+  memset(fSigmaPar  ,0,fNGloPar*sizeof(Double_t));
+  memset(fProcPnt   ,0,fNGloPar*sizeof(Int_t));
+  //
+  for (int i=fNGloPar;i--;) {
+    fGlo2CGlo[i]=fCGlo2Glo[i] = -1;
+    fIsLinear[i] = kTRUE;
+  }
+  //
+  return 1;
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::ResetConstraints()
+{
+  // reset constraints tree. Allows to redefine the constraints for preprocessed data 
+  CloseConsRecStorage();
+  InitConsRecStorage(kFALSE);
+  fNGloConstraints = 0;
+  AliInfo("Constraints are reset");
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::ImposeDataRecFile(const char* fname)
+{
+  CloseDataRecStorage();
+  SetDataRecFName(fname);
+  return InitDataRecStorage(kTRUE); // open in read mode
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::ImposeConsRecFile(const char* fname)
+{
+  CloseConsRecStorage();
+  SetConsRecFName(fname);
+  return InitConsRecStorage(kTRUE); // open in read mode
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::InitDataRecStorage(Bool_t read)
+{
+  // initialize the buffer for processed measurements records
+  //
+  if (fDataRecFile) {AliInfo("Data Records File is already initialized"); return kFALSE;} 
+  //
+  if (!fRecord) fRecord = new AliMillePedeRecord();
+  //
+  fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate");
+  if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
+  //
+  AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
+  if (read) {
+    fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data");
+    if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
+    fTreeData->SetBranchAddress("Record_Data",&fRecord);
+    AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries()));
+  }
+  else {
+    fTreeData = new TTree("AliMillePedeRecords_Data","Data Records for AliMillePede2");
+    fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99);
+  }
+  fCurrRecDataID = -1;
+  //
+  return kTRUE;
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::InitConsRecStorage(Bool_t read)
+{
+  // initialize the buffer for processed measurements records
+  //
+  if (fConsRecFile) {AliInfo("Constraints Records File is already initialized"); return kFALSE;} 
+  //
+  if (!fRecord) fRecord = new AliMillePedeRecord();
+  //
+  fConsRecFile = TFile::Open(GetConsRecFName(),read ? "":"recreate");
+  if (!fConsRecFile) {AliInfo(Form("Failed to initialize constraints records file %s",GetConsRecFName())); return kFALSE;}
+  //
+  AliInfo(Form("File %s used for constraints records",GetConsRecFName()));
+  if (read) {
+    fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints");
+    if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;}
+    fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord);
+    AliInfo(Form("Found %d constraints records",fTreeConstr->GetEntries()));
+    //
+  }
+  else {
+    //
+    fTreeConstr = new TTree("AliMillePedeRecords_Constraints","Constraints Records for AliMillePede2");
+    fTreeConstr->Branch("Record_Constraints","AliMillePedeRecord",&fRecord,32000,99);
+  }
+  fCurrRecConstrID = -1;
+  //
+  return kTRUE;
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::CloseDataRecStorage()
+{
+  if (fTreeData) {
+    if (fDataRecFile->IsWritable()) {
+      fDataRecFile->cd();
+      fTreeData->Write();
+    }
+    delete fTreeData;  
+    fTreeData = 0;
+    fDataRecFile->Close();
+    delete fDataRecFile;
+    fDataRecFile = 0;
+  }
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::CloseConsRecStorage()
+{
+  if (fTreeConstr) {
+    if (fConsRecFile->IsWritable()) {
+      fConsRecFile->cd();
+      fTreeConstr->Write();
+    }
+    delete fTreeConstr;  
+    fTreeConstr = 0;
+    fConsRecFile->Close();
+    delete fConsRecFile;
+    fConsRecFile = 0;
+  }
+  //
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::ReadNextRecordData()
+{
+  // read next data record (if any)
+  if (!fTreeData || ++fCurrRecDataID >= fTreeData->GetEntries()) { fCurrRecDataID--; return kFALSE;}
+  fTreeData->GetEntry(fCurrRecDataID);
+  return kTRUE;
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::ReadNextRecordConstraint()
+{
+  // read next constraint record (if any)
+  if (!fTreeConstr || ++fCurrRecConstrID >= fTreeConstr->GetEntries()) { fCurrRecConstrID--; return kFALSE;}
+  fTreeConstr->GetEntry(fCurrRecConstrID);
+  return kTRUE;
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
+{
+  if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+  //
+  // write data of single measurement
+  if (lSigma<=0.0) { // If parameter is fixed, then no equation
+    for (int i=fNLocPar; i--;) derlc[i] = 0.0;
+    for (int i=fNGloPar; i--;) dergb[i] = 0.0;
+    return;
+  }
+  //
+  fRecord->AddResidual(lMeas);
+  //
+  // Retrieve local param interesting indices
+  for (int i=0;i<fNLocPar;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;}
+  //
+  fRecord->AddWeight( 1.0/lSigma/lSigma );
+  //
+  // Idem for global parameters
+  for (int i=0;i<fNGloPar;i++) if (dergb[i]!=0.0) {fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;}
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
+                                    double *derlc,int nlc,double lMeas,double lSigma)
+{      
+  // write data of single measurement
+  if (lSigma<=0.0) { // If parameter is fixed, then no equation
+    for (int i=nlc;i--;)  derlc[i] = 0.0;
+    for (int i=ngb;i--;)  dergb[i] = 0.0;
+    return;
+  }
+  //
+  if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+  //
+  fRecord->AddResidual(lMeas);
+  //
+  // Retrieve local param interesting indices
+  for (int i=0;i<nlc;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;}
+  //
+  fRecord->AddWeight( 1./lSigma/lSigma );
+  //
+  // Idem for global parameters
+  for (int i=0;i<ngb;i++) if (dergb[i]!=0.0) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} 
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetGlobalConstraint(double *dergb, double val)
+{      
+  // Define a constraint equation.
+  if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
+  //
+  fRecord->Reset();
+  fRecord->AddResidual(val);
+  fRecord->AddWeight(val);   // dummy
+  for (int i=0; i<fNGloPar; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(i,dergb[i]);
+  fNGloConstraints++;
+  SaveRecordConstraint();
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val)
+{      
+  // Define a constraint equation.
+  if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
+  fRecord->Reset();
+  fRecord->AddResidual(val);
+  fRecord->AddWeight(val);   // dummy
+  for (int i=0; i<ngb; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(indgb[i],dergb[i]);
+  fNGloConstraints++;
+  SaveRecordConstraint();
+  //
+}
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::LocalFit(double *localParams)
+{
+  /*
+    Perform local parameters fit once all the local equations have been set
+    -----------------------------------------------------------
+    localParams = (if !=0) will contain the fitted track parameters and
+    related errors
+  */
+  static int     nRefSize = 0;
+  static TArrayI RefLoc,RefGlo,nRefLoc,nRefGlo;
+  int nPoints = 0;
+  //
+  AliSymMatrix    &MatCLoc    = *fMatCLoc;
+  AliMatrixSq     &MatCGlo    = *fMatCGlo;
+  TMatrixD        &MatCGloLoc = *fMatCGloLoc;
+  //
+  memset(fVecBLoc,0,fNLocPar*sizeof(double));
+  MatCLoc.Reset();     
+  //
+  int cnt=0;
+  int recSz = fRecord->GetSize();
+  //
+  while(cnt<recSz) {  // Transfer the measurement records to matrices
+    //
+    // extract addresses of residual, weight and pointers on local and global derivatives for each point
+    if (nRefSize<=nPoints) {
+      nRefSize = 2*(nPoints+1); 
+      RefLoc.Set(nRefSize); 
+      RefGlo.Set(nRefSize);
+      nRefLoc.Set(nRefSize); 
+      nRefGlo.Set(nRefSize);
+    }
+    //
+    RefLoc[nPoints]  = ++cnt;
+    int nLoc = 0;
+    while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;}
+    nRefLoc[nPoints] = nLoc;
+    //
+    RefGlo[nPoints]  = ++cnt;
+    int nGlo = 0;
+    while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;} 
+    nRefGlo[nPoints] = nGlo;
+    //
+    nPoints++;
+  }
+  double vl;
+  //
+  for (int ip=nPoints;ip--;) {  // Transfer the measurement records to matrices
+    double  resid  = fRecord->GetValue( RefLoc[ip]-1 );
+    double  weight = fRecord->GetValue( RefGlo[ip]-1 );    
+    double *derLoc = fRecord->GetValue()+RefLoc[ip];
+    double *derGlo = fRecord->GetValue()+RefGlo[ip];
+    int    *indLoc = fRecord->GetIndex()+RefLoc[ip];
+    int    *indGlo = fRecord->GetIndex()+RefGlo[ip];
+    //
+    for (int i=nRefGlo[ip];i--;) {       // suppress the global part (only relevant with iterations)
+      int iID = indGlo[i];              // Global param indice
+      if ( fSigmaPar[iID] <= 0.) continue;                                    // fixed parameter RRRCheck
+      if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);  // linear parameter
+      else                resid -= derGlo[i]*fDeltaPar[iID];                  // nonlinear parameter
+    }
+    //
+    for (int i=nRefLoc[ip];i--;) {         // Fill local matrix and vector
+      int iID = indLoc[i];              // Local param indice (the matrix line) 
+      if ( (vl=weight*resid*derLoc[i])!=0 ) fVecBLoc[iID] += vl;
+      //                       
+      for (int j=i+1;j--;) {           // Symmetric matrix, don't bother j>i coeffs
+       int jID = indLoc[j];    
+       if ( (vl=weight*derLoc[i]*derLoc[j])!=0 ) MatCLoc(iID,jID) += vl;
+      }
+    }   
+    //
+  } // end of the transfer of the measurement record to matrices
+  //
+  // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
+  if (!MatCLoc.SolveChol(fVecBLoc,kTRUE)) {
+    AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
+    if (!MatCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { 
+      AliInfo("Failed to solve locals by Gaussian Elimination, skip...");
+      MatCLoc.Print("d");
+      return 0; // failed to solve
+    }
+  }
+  //
+  // If requested, store the track params and errors
+  if (localParams) for (int i=fNLocPar; i--;) {
+      localParams[2*i]   = fVecBLoc[i];
+      localParams[2*i+1] = TMath::Sqrt(TMath::Abs(MatCLoc.QuerryDiag(i)));
+    }
+  //
+  float lChi2 = 0;
+  int   nEq   = 0;
+  //
+  for (int ip=nPoints;ip--;) {  // Calculate residuals
+    double  resid  = fRecord->GetValue( RefLoc[ip]-1 );
+    double  weight = fRecord->GetValue( RefGlo[ip]-1 );    
+    double *derLoc = fRecord->GetValue()+RefLoc[ip];
+    double *derGlo = fRecord->GetValue()+RefGlo[ip];
+    int    *indLoc = fRecord->GetIndex()+RefLoc[ip];
+    int    *indGlo = fRecord->GetIndex()+RefGlo[ip];
+    //
+    // Suppress local and global contribution in residuals;
+    for (int i=nRefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ];     // local part 
+    //
+    for (int i=nRefGlo[ip];i--;) {                                             // global part
+      int iID = indGlo[i];
+      if ( fSigmaPar[iID] <= 0.) continue;                                     // fixed parameter RRRCheck
+      if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);   // linear parameter
+      else                resid -= derGlo[i]*fDeltaPar[iID];                   // nonlinear parameter
+    }
+    //
+    // reject the track if the residual is too large (outlier)
+    if ( (TMath::Abs(resid) >= fResCutInit && fIter ==1 ) ||
+        (TMath::Abs(resid) >= fResCut     && fIter > 1)) {
+      fNLocFitsRejected++;      
+      return 0;
+    }
+    // 
+    lChi2 += weight*resid*resid ; // total chi^2
+    nEq++;                        // number of equations                       
+  } // end of Calculate residuals
+  //
+  //
+  int nDoF = nEq-fNLocPar;
+  lChi2 = (nDoF>0) ? lChi2/nDoF : 0;  // Chi^2/dof  
+  //
+  if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
+    fNLocFitsRejected++;      
+    return 0;
+  }
+  //
+  fNLocFits++;
+  fNLocEquations += nEq;
+  //
+  //  local operations are finished, track is accepted 
+  //  We now update the global parameters (other matrices)
+  //
+  int nGloInFit = 0;
+  //
+  for (int ip=nPoints;ip--;) {  // Update matrices
+    double  resid  = fRecord->GetValue( RefLoc[ip]-1 );
+    double  weight = fRecord->GetValue( RefGlo[ip]-1 );    
+    double *derLoc = fRecord->GetValue()+RefLoc[ip];
+    double *derGlo = fRecord->GetValue()+RefGlo[ip];
+    int    *indLoc = fRecord->GetIndex()+RefLoc[ip];
+    int    *indGlo = fRecord->GetIndex()+RefGlo[ip];
+    //
+    for (int i=nRefGlo[ip];i--;) {    // suppress the global part
+      int iID = indGlo[i];           // Global param indice
+      if ( fSigmaPar[iID] <= 0.) continue;                                         // fixed parameter RRRCheck
+      if (fIsLinear[iID] == 0)  resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter
+      else                      resid -= derGlo[i]*fDeltaPar[iID];                 // nonlinear parameter
+    }
+    //
+    for (int ig=nRefGlo[ip];ig--;) {
+      int iIDg = indGlo[ig];   // Global param indice (the matrix line)          
+      if ( fSigmaPar[iIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
+      if ( (vl = weight*resid*derGlo[ig])!=0 ) fVecBGlo[ iIDg ] += vl;
+      //
+      // First of all, the global/global terms (exactly like local matrix)
+      for (int jg=ig+1;jg--;) {        // MatCGlo is symmetric by construction  
+       int jIDg = indGlo[jg];
+       if ( fSigmaPar[jIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
+       if ( (vl = weight*derGlo[ig]*derGlo[jg])!=0 ) MatCGlo(iIDg,jIDg) += vl; 
+      } 
+      //
+      // Now we have also rectangular matrices containing global/local terms.
+      int iCIDg = fGlo2CGlo[iIDg];  // compressed Index of index          
+      if (iCIDg == -1) {
+       for (int k=fNLocPar;k--;) MatCGloLoc(nGloInFit,k) = 0.0;  // reset the row
+       iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
+       fCGlo2Glo[nGloInFit++] = iIDg;
+      }
+      //
+      for (int il=nRefLoc[ip];il--;) {
+       int iIDl = indLoc[il];
+       if ( (vl = weight*derGlo[ig]*derLoc[il])!=0) MatCGloLoc(iCIDg,iIDl) += vl;
+      }
+      // update counter
+      fProcPnt[iIDg]++;
+      //
+    }
+  } // end of Update matrices
+  //
+  // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T
+  // and       fVecBGlo -= fMatCGloLoc * fVecBLoc
+  //
+  double vll;
+  for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) {
+    int iIDg = fCGlo2Glo[iCIDg];
+    //
+    vl = 0;
+    for (int kl=0;kl<fNLocPar;kl++) if ( (vll = MatCGloLoc(iCIDg,kl))!=0) vl += vll*fVecBLoc[kl];
+    if  (vl!=0) fVecBGlo[iIDg] -= vl;
+    //
+    for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) {  
+      int jIDg = fCGlo2Glo[jCIDg];
+      //
+      vl = 0;
+      for (int kl=0;kl<fNLocPar;kl++) {
+       // diag terms
+       if ( (vll=MatCGloLoc(iCIDg,kl)*MatCGloLoc(jCIDg,kl))!=0 ) vl += MatCLoc.QuerryDiag(kl)*vll;
+       //
+       // off-diag terms
+       for (int ll=0;ll<kl;ll++) {
+         if ( (vll=MatCGloLoc(iCIDg,kl)*MatCGloLoc(jCIDg,ll))!=0 ) vl += MatCLoc(kl,ll)*vll;
+         if ( (vll=MatCGloLoc(iCIDg,ll)*MatCGloLoc(jCIDg,kl))!=0 ) vl += MatCLoc(kl,ll)*vll;
+       }
+      }
+      if (vl!=0) MatCGlo(iIDg,jIDg) -= vl; 
+      //
+    }
+  }
+  //
+  // reset compressed index array
+  //
+  for (int i=nGloInFit;i--;) { fGlo2CGlo[ fCGlo2Glo[i] ] = -1; fCGlo2Glo[i] = -1;}
+  //
+  return 1;
+}
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull)
+{
+  // performs a requested number of global iterations
+  fIter = 1;
+  //
+  TStopwatch sw; sw.Start();
+  //
+  int res = 0;
+  AliInfo("Starting Global fit.");
+  while (fIter<=fMaxIter) {
+    //
+    res = GlobalFitIteration();
+    if (!res) break;
+    //
+    if (fChi2CutFactor != fChi2CutRef) {    
+      fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
+      if (fChi2CutFactor < 1.2*fChi2CutRef) {
+       fChi2CutFactor = fChi2CutRef;
+       fIter = fMaxIter - 1;     // Last iteration
+      }
+    }
+    fIter++;
+  }
+  //
+  sw.Stop();
+  AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime()));  
+  if (!res) return 0;
+  //
+  if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i]; 
+  //
+  if (fGloSolveStatus==gkInvert) { // errors on params are available
+    if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QuerryDiag(i))) : 0.;
+    if (pull)  for (int i=fNGloPar;i--;) pull[i]  = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QuerryDiag(i))>0. && fSigmaPar[i]>0 
+                                          ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QuerryDiag(i)) : 0;
+  }
+  //
+  return 1;
+}
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::GlobalFitIteration()
+{
+  // perform global parameters fit once all the local equations have been fitted
+  //
+  AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor));
+  //
+  if (!fNGloPar || !fTreeData) {
+    AliInfo("No data was stored, aborting iteration");
+    return 0;
+  }
+  TStopwatch sw; sw.Start();
+  //
+  if (!fConstrUsed) {
+    fConstrUsed = new Bool_t[fNGloConstraints];
+    memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t));
+  }
+  // Reset all info specific for this step
+  AliMatrixSq& MatCGlo = *fMatCGlo;
+  MatCGlo.Reset();
+  memset(fProcPnt,0,fNGloPar*sizeof(Int_t));
+  //
+  fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0;
+  //
+  // if needed, readjust the size of the global vector (for matrices this is done automatically)
+  if (!fVecBGlo || fNGloSize!=fNGloPar+fNGloConstraints) {
+    delete[] fVecBGlo;   // in case some constraint was added between the two manual iterations
+    fNGloSize = fNGloPar+fNGloConstraints;
+    fVecBGlo = new Double_t[fNGloSize];
+  }
+  memset(fVecBGlo,0,fNGloSize*sizeof(double));
+  //
+  fNLocFits         = 0;
+  fNLocFitsRejected = 0;
+  fNLocEquations      = 0;
+  //
+  //  Process data records and build the matrices
+  Long_t ndr = fTreeData->GetEntries();
+  AliInfo(Form("Building the Global matrix from %d data records",ndr));
+  for (Long_t i=0;i<ndr;i++) {
+    ReadRecordData(i);
+    LocalFit();
+  }
+  //
+  fNGloFix = 0;
+  for (int i=fNGloPar;i--;) fDiagCGlo[i] = MatCGlo.QuerryDiag(i); // save the diagonal elements  
+  //
+  //  Account for fixed parameters
+  for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
+    if ( fSigmaPar[i] <= 0. || fProcPnt[i]<fMinPntValid) {
+      fNGloFix++; 
+      fProcPnt[i] *= -1;
+      fVecBGlo[i] = 0.;
+      if (IsGlobalMatSparse()) {
+       AliVectorSparse& row = *((AliMatrixSparse*)fMatCGlo)->GetRow(i);
+       row.Clear();
+       row(i) = float(fNLocEquations);
+       for (int j=i+1;j<fNGloPar;j++) ((AliMatrixSparse*)fMatCGlo)->Zero(i,j);
+      }
+      else 
+       for (int j=fNGloPar;j--;) if (MatCGlo.Querry(i,j)) MatCGlo(i,j)=0;
+       MatCGlo.DiagElem(i) = float(fNLocEquations);
+    }
+    else MatCGlo.DiagElem(i) += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
+  }
+  //
+  // add constraint equations
+  int nVar = fNGloPar;                    // Current size of global matrix     
+  for (int i=0; i<fNGloConstraints; i++) {
+    ReadRecordConstraint(i);
+    double val   = fRecord->GetValue(0);  
+    int    *indV = fRecord->GetIndex()+2;
+    double *der  = fRecord->GetValue()+2;
+    int    csize = fRecord->GetSize()-2;
+    //
+    // check if after suppression of fixed variables this there are non-0 derivatives
+    int NSuppressed = 0;
+    for (int j=csize;j--;)  if (fProcPnt[indV[j]]<1) NSuppressed++;
+    //
+    if (NSuppressed==csize) {
+      AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize));
+      //
+      // was this constraint ever created ? 
+      if ( fConstrUsed[i] ) {
+       // to avoid empty row impose dummy constraint on "Lagrange multiplier"
+       MatCGlo.DiagElem(nVar) = 1.;
+       fVecBGlo[nVar++] = 0;
+      }
+      continue;
+    }
+    //
+    for (int j=csize; j--;) {
+      int jID = indV[j];
+      val -= der[j]*(fInitPar[jID]+fDeltaPar[jID]);
+      if (fProcPnt[jID]<1) continue;                      // this parameter was fixed, don't put it into constraint
+      MatCGlo(nVar,jID) = float(fNLocEquations)*der[j];   // fMatCGlo is symmetric, only lower triangle is filled  
+    }
+    //
+    if (MatCGlo.QuerryDiag(nVar)) MatCGlo.DiagElem(nVar) = 0.0;
+    fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? 
+    fConstrUsed[i] = kTRUE;
+  }
+  //
+  AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals",
+              fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix));
+
+  //
+  fGloSolveStatus = SolveGlobalMatEq();                     // obtain solution for this step
+  //
+  sw.Stop();
+  AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime()));
+  if (fGloSolveStatus==gkFailed) return 0;
+  //
+  for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
+  //
+  return 1;
+}
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::SolveGlobalMatEq()
+{
+  //
+  // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo
+  //
+  if (!fgIsMatGloSparse) {
+    //
+    if (fNGloConstraints==0) { // pos-def systems are faster to solve by Cholesky
+      if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, kTRUE) ) return gkInvert;
+      else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation");
+    }
+    //
+    if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert;
+    else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods");
+  }
+  // try to solve by minres
+  TVectorD SOL(fNGloSize);
+  //
+  AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo);
+  if (!slv) return gkFailed;
+  //
+  Bool_t res = kFALSE;
+  if      (fgIterSol == AliMinResSolve::kSolMinRes) 
+    res =  slv->SolveMinRes(SOL,fgMinResCondType,fgMinResMaxIter,fgMinResTol);
+  else if (fgIterSol == AliMinResSolve::kSolFGMRes) 
+    res =  slv->SolveFGMRES(SOL,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV);
+  else 
+    AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers));
+  //
+  if (!res) return gkFailed;
+  for (int i=fNGloSize;i--;) fVecBGlo[i] = SOL[i];
+  return gkMinRes;
+  //
+}
+
+//_____________________________________________________________________________________________
+Float_t AliMillePede2::Chi2DoFLim(int nSig, int nDoF) const
+{
+  /// return the limit in chi^2/nd for n sigmas stdev authorized
+  // Only n=1, 2, and 3 are expected in input
+  int lNSig;
+  float sn[3]        = {0.47523, 1.690140, 2.782170};
+  float table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
+                        1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
+                        1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
+                        1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
+                       {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
+                        1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
+                        1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
+                        1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
+                       {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
+                        2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
+                        2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
+                        1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
+  
+  if (nDoF < 1) {
+    return 0.0;
+  }
+  else {  
+    lNSig = TMath::Max(1,TMath::Min(nSig,3));
+    
+    if (nDoF <= 30) {    
+      return table[lNSig-1][nDoF-1];
+    }
+    else { // approximation
+      return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
+             (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
+    }
+  }
+}
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::SetIterations(double lChi2CutFac)
+{
+  // Number of iterations is calculated from lChi2CutFac 
+  fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
+  //
+  AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
+  fIter = 1; // Initializes the iteration process
+  return 1;
+}
+
+//_____________________________________________________________________________________________
+Double_t AliMillePede2::GetParError(int iPar) const
+{
+  // return error for parameter iPar
+  if (fGloSolveStatus==gkInvert) {
+    double res = fMatCGlo->QuerryDiag(iPar);
+    if (res>=0) return TMath::Sqrt(res);
+  } 
+  return -1.;
+}
+
+
+//_____________________________________________________________________________________________
+Int_t AliMillePede2::PrintGlobalParameters() const
+{
+  ///  Print the final results into the logfile
+  double lError = 0.;
+  double lGlobalCor =0.;
+       
+  AliInfo("");
+  AliInfo("   Result of fit for global parameters");
+  AliInfo("   ===================================");
+  AliInfo("    I       initial       final       differ        lastcor        error       gcor       Npnt");
+  AliInfo("----------------------------------------------------------------------------------------------");
+  //
+  for (int i=0; i<fNGloPar; i++) {
+    lError = GetParError(i);
+    lGlobalCor = 0.0;
+    //         
+    double dg;
+    if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QuerryDiag(i)) *fDiagCGlo[i]) > 0) {    
+      lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i])));
+      AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d",
+                  i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]));
+    }
+    else {    
+      AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF\t %6d",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],
+                  fDeltaPar[i],fVecBGlo[i],fProcPnt[i]));
+    }
+  }
+  return 1;
+}
diff --git a/STEER/AliMillePede2.h b/STEER/AliMillePede2.h
new file mode 100644 (file)
index 0000000..9c91225
--- /dev/null
@@ -0,0 +1,278 @@
+#ifndef ALIMILLEPEDE2_H
+#define ALIMILLEPEDE2_H
+
+#include <TObject.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TMatrixD.h>
+#include "AliSymMatrix.h"
+#include "AliMatrixSparse.h"
+#include "AliMinResSolve.h"
+class AliLog;
+class TStopwatch;
+
+/**********************************************************************************************/
+/* AliMillePedeRecords: class to store the data of single track processing                    */
+/* Format: for each measured point the data is stored consequtively                           */
+/* INDEX                                                      VALUE                           */
+/* -1                                                         residual                        */
+/* Local_param_id                                             dResidual/dLocal_param          */
+/* ...                                                        ...                             */
+/* -2                                                         weight of the measurement       */
+/* Global_param_od                                            dResidual/dGlobal_param         */
+/* ...                                                        ...                             */
+/*                                                                                            */
+/* The records for all processed tracks are stored in the temporary tree in orgder to be      */
+/* reused for multiple iterations of MillePede                                                */
+/*                                                                                            */
+/**********************************************************************************************/
+
+class AliMillePedeRecord : public TObject
+{
+ public:
+  AliMillePedeRecord();
+  AliMillePedeRecord(const AliMillePedeRecord& src);
+  AliMillePedeRecord& operator=(const AliMillePedeRecord& rhs);
+  //
+  virtual    ~AliMillePedeRecord();
+  void       Reset()                                             {fSize = 0;}
+  void       Print(const Option_t *opt="")                 const;
+  //
+  Int_t      GetSize()                                     const {return fSize;}
+  Int_t     *GetIndex()                                    const {return fIndex;}
+  Int_t      GetIndex(int i)                               const {return fIndex[i];}
+  //
+  void       GetIndexValue(int i,int &ind,double &val)     const {ind=fIndex[i]; val=fValue[i];}
+  void       AddIndexValue(int ind, double val);
+  void       AddResidual(double val)                             {AddIndexValue(-1,val);}
+  void       AddWeight(double val)                               {AddIndexValue(-2,val);}
+  Bool_t     IsResidual(int i)                             const {return fIndex[i]==-1;}
+  Bool_t     IsWeight(int i)                               const {return fIndex[i]==-2;}
+  //
+  Double_t  *GetValue()                                    const {return fValue;}
+  Double_t   GetValue(int i)                               const {return fValue[i];}
+  //
+ protected:
+  Int_t      GetBufferSize()                               const {return GetUniqueID();}
+  void       SetBufferSize(int sz)                               {SetUniqueID(sz);}
+  void       Expand(int bfsize);
+  //
+ protected:
+  Int_t      fSize;                             // size of the record
+  Int_t   *  fIndex;                            //[fSize] index of variables
+  Double_t*  fValue;                            //[fSize] array of values: derivs,residuals
+  //
+  ClassDef(AliMillePedeRecord,1)                // Record of track residuals and local/global deriavtives
+};
+
+inline void  AliMillePedeRecord::AddIndexValue(int ind, double val) 
+{
+  if (fSize>=GetBufferSize()) Expand(2*(fSize+1));
+  fIndex[fSize]=ind; 
+  fValue[fSize++]=val;
+}
+
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------
+
+class AliMillePede2: public TObject
+{
+ public:
+  //
+  enum {gkFailed,gkInvert,gkMinRes};    // used global matrix solution methods
+  //
+  AliMillePede2();
+  AliMillePede2(const AliMillePede2& src);
+  virtual ~AliMillePede2();
+  AliMillePede2& operator=(const AliMillePede2& )             {printf("Dummy\n"); return *this;}
+  //
+  Int_t                InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1.);
+  //
+  Int_t                GetNGloPar()                     const {return fNGloPar;}
+  Int_t                GetNLocPar()                     const {return fNLocPar;}
+  Long_t               GetNLocalEquations()             const {return fNLocEquations;}
+  Int_t                GetCurrentIteration()            const {return fIter;}
+  Int_t                GetNMaxIterations()              const {return fMaxIter;}
+  Int_t                GetNStdDev()                     const {return fNStdDev;} 
+  Int_t                GetNGlobalConstraints()          const {return fNGloConstraints;}
+  Long_t               GetNLocalFits()                  const {return fNLocFits;}
+  Long_t               GetNLocalFitsRejected()          const {return fNLocFitsRejected;}
+  Int_t                GetNGlobalsFixed()               const {return fNGloFix;}
+  Int_t                GetGlobalSolveStatus()           const {return fGloSolveStatus;}
+  Float_t              GetChi2CutFactor()               const {return fChi2CutFactor;}
+  Float_t              GetChi2CutRef()                  const {return fChi2CutRef;}
+  Float_t              GetResCurInit()                  const {return fResCutInit;}
+  Float_t              GetResCut()                      const {return fResCut;}
+  Int_t                GetMinPntValid()                 const {return fMinPntValid;}
+  Int_t                GetProcessedPoints(Int_t i)      const {return fProcPnt[i];}
+  Int_t*               GetProcessedPoints()             const {return fProcPnt;}
+  //
+  AliMatrixSq*         GetGlobalMatrix()                const {return fMatCGlo;}
+  Double_t*            GetGlobals()                     const {return fVecBGlo;}
+  Double_t*            GetDeltaPars()                   const {return fDeltaPar;}
+  Double_t*            GetInitPars()                    const {return fInitPar;}
+  Double_t*            GetSigmaPars()                   const {return fSigmaPar;}
+  Bool_t*              GetIsLinear()                    const {return fIsLinear;}        
+  Double_t             GetGlobal(Int_t i)               const {return fVecBGlo[i];}
+  Double_t             GetInitPar(Int_t i)              const {return fInitPar[i];}
+  Double_t             GetSigmaPar(Int_t i)             const {return fSigmaPar[i];}
+  Bool_t               GetIsLinear(Int_t i)             const {return fIsLinear[i];}
+  static Bool_t        IsGlobalMatSparse()                    {return fgIsMatGloSparse;}
+  //
+  void                 SetNGloPar(Int_t n)                    {fNGloPar = n;}
+  void                 SetNLocPar(Int_t n)                    {fNLocPar = n;}
+  void                 SetNMaxIterations(Int_t n=10)          {fMaxIter = n;}
+  void                 SetNStdDev(Int_t n)                    {fNStdDev = n;}
+  void                 SetChi2CutFactor(Float_t v)            {fChi2CutFactor = v;}
+  void                 SetChi2CutRef(Float_t v)               {fChi2CutRef = v;}
+  void                 SetResCurInit(Float_t v)               {fResCutInit = v;}
+  void                 SetResCut(Float_t v)                   {fResCut = v;}
+  void                 SetMinPntValid(Int_t n)                {fMinPntValid = n;}
+  static void          SetGlobalMatSparse(Bool_t v=kTRUE)     {fgIsMatGloSparse = v;}
+  //
+  void                 SetInitPars(const Double_t* par)       {memcpy(fInitPar,par,fNGloPar*sizeof(Double_t));}
+  void                 SetSigmaPars(const Double_t* par)      {memcpy(fSigmaPar,par,fNGloPar*sizeof(Double_t));}
+  void                 SetInitPar(Int_t i,Double_t par)       {fInitPar[i] = par;;}
+  void                 SetSigmaPar(Int_t i,Double_t par)      {fSigmaPar[i] = par;}
+  //
+  Int_t                GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
+  Int_t                GlobalFitIteration();
+  Int_t                SolveGlobalMatEq();
+  static void          SetMinResPrecondType(Int_t tp=0)       {fgMinResCondType = tp;}
+  static void          SetMinResTol(Double_t val=1e-12)       {fgMinResTol = val;}
+  static void          SetMinResMaxIter(Int_t val=2000)       {fgMinResMaxIter = val;}
+  static void          SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
+  static void          SetNKrylovV(Int_t val=60)              {fgNKrylovV = val;}
+  //
+  static Int_t         GetMinResPrecondType()                 {return fgMinResCondType;}
+  static Double_t      GetMinResTol()                         {return fgMinResTol;}
+  static Int_t         GetMinResMaxIter()                     {return fgMinResMaxIter;}
+  static Int_t         GetIterSolverType()                    {return fgIterSol;}
+  static Int_t         GetNKrylovV()                          {return fgNKrylovV;}
+  //
+  Double_t             GetParError(int iPar)           const;
+  Int_t                PrintGlobalParameters()         const;
+  //
+  //
+  Int_t                SetIterations(double lChi2CutFac);
+
+  //
+  // constraints
+  void                 SetGlobalConstraint(double *dergb, double val);
+  void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val);
+  //
+  // processing of the local measurement
+  void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
+  void                 SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc, 
+                                       double *derlc,int nlc,double lMeas,double lSigma);
+  //
+  // manipilation with processed data and costraints records and its buffer
+  void                 SetDataRecFName(const char* flname)   {fDataRecFName = flname;}
+  const Char_t*        GetDataRecFName()               const {return fDataRecFName.Data();}
+  void                 SetConsRecFName(const char* flname)   {fConstrRecFName = flname;}
+  const Char_t*        GetConsRecFName()               const {return fConstrRecFName.Data();}
+
+  Bool_t               InitDataRecStorage(Bool_t read=kFALSE);
+  Bool_t               InitConsRecStorage(Bool_t read=kFALSE);
+  Bool_t               ImposeDataRecFile(const char* fname);
+  Bool_t               ImposeConsRecFile(const char* fname);
+  void                 CloseDataRecStorage();
+  void                 CloseConsRecStorage();
+  void                 ReadRecordData(Long_t recID);
+  void                 ReadRecordConstraint(Long_t recID);
+  Bool_t               ReadNextRecordData();
+  Bool_t               ReadNextRecordConstraint();
+  void                 SaveRecordData();
+  void                 SaveRecordConstraint();
+  void                 ResetConstraints();
+  AliMillePedeRecord*  GetRecord()                      const {return fRecord;}
+  //
+  Float_t              Chi2DoFLim(int nSig, int nDoF)   const;
+  //
+  // aliases for compatibility with millipede1
+  void                 SetParSigma(Int_t i,Double_t par)      {SetSigmaPar(i,par);}
+  void                 SetGlobalParameters(Double_t *par)     {SetInitPars(par);}
+  //
+ protected:
+  //
+  Int_t                LocalFit(double *localParams=0);
+  //
+ protected:
+  //
+  Int_t                 fNLocPar;                        // number of local parameters
+  Int_t                 fNGloPar;                        // number of global parameters
+  Int_t                 fNGloSize;                       // final size of the global matrix (NGloPar+NConstraints)
+  //
+  Long_t                fNLocEquations;                  // Number of local equations 
+  Int_t                 fIter;                           // Current iteration
+  Int_t                 fMaxIter;                        // Maximum number of iterations
+  Int_t                 fNStdDev;                        // Number of standard deviations for chi2 cut 
+  Int_t                 fNGloConstraints;                // Number of constraint equations
+  Long_t                fNLocFits;                       // Number of local fits
+  Long_t                fNLocFitsRejected;               // Number of local fits rejected
+  Int_t                 fNGloFix;                        // Number of globals fixed by user
+  Int_t                 fGloSolveStatus;                 // Status of global solver at current step
+  //
+  Float_t               fChi2CutFactor;                  // Cut factor for chi2 cut to accept local fit 
+  Float_t               fChi2CutRef;                     // Reference cut for chi2 cut to accept local fit 
+  Float_t               fResCutInit;                     // Cut in residual for first iterartion
+  Float_t               fResCut;                         // Cut in residual for other iterartiona
+  Int_t                 fMinPntValid;                    // min number of points for global to vary
+  //
+  Int_t                *fProcPnt;                        // N of processed points per global variable
+  Double_t             *fVecBLoc;                        // Vector B local (parameters) 
+  Double_t             *fDiagCGlo;                       // Initial diagonal elements of C global matrix
+  Double_t             *fVecBGlo;                        // Vector B global (parameters)
+  //
+  Double_t             *fInitPar;                        // Initial global parameters
+  Double_t             *fDeltaPar;                       // Variation of global parameters
+  Double_t             *fSigmaPar;                       // Sigma of allowed variation of global parameter
+  //
+  Bool_t               *fIsLinear;                       // Flag for linear parameters
+  Bool_t               *fConstrUsed;                     // Flag for used constraints
+  //
+  Int_t                *fGlo2CGlo;                       // global ID to compressed ID buffer
+  Int_t                *fCGlo2Glo;                       // compressed ID to global ID buffer
+  //
+  // Matrices
+  AliSymMatrix         *fMatCLoc;                        // Matrix C local
+  AliMatrixSq          *fMatCGlo;                        // Matrix C global
+  TMatrixD             *fMatCGloLoc;                     // Rectangular matrix C g*l 
+  //
+  // processed data record bufferization   
+  TString               fDataRecFName;                   // Name of File for data records               
+  AliMillePedeRecord   *fRecord;                         // Buffer of measurements records
+  TFile                *fDataRecFile;                    // File of processed measurements records
+  TTree                *fTreeData;                       // Tree of processed measurements records
+  //
+  TString               fConstrRecFName;                 // Name of File for constraints records               
+  TTree                *fTreeConstr;                     // Tree of constraint records
+  TFile                *fConsRecFile;                    // File of processed constraints records
+  Long_t                fCurrRecDataID;                  // ID of the current data record
+  Long_t                fCurrRecConstrID;                // ID of the current constraint record
+  //
+  static Bool_t         fgIsMatGloSparse;                // Type of the global matrix (sparse ...)
+  static Int_t          fgMinResCondType;                // Type of the preconditioner for MinRes method 
+  static Double_t       fgMinResTol;                     // Tolerance for MinRes solution
+  static Int_t          fgMinResMaxIter;                 // Max number of iterations for the MinRes method
+  static Int_t          fgIterSol;                       // type of iterative solution: MinRes or FGMRES
+  static Int_t          fgNKrylovV;                      // size of Krylov vectors buffer in FGMRES
+  //
+  ClassDef(AliMillePede2,0)
+};
+
+//_____________________________________________________________________________________________
+inline void AliMillePede2::ReadRecordData(Long_t recID)       {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
+
+//_____________________________________________________________________________________________
+inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
+
+//_____________________________________________________________________________________________
+inline void AliMillePede2::SaveRecordData()                   {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
+
+//_____________________________________________________________________________________________
+inline void AliMillePede2::SaveRecordConstraint()             {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
+
+#endif
diff --git a/STEER/AliMinResSolve.cxx b/STEER/AliMinResSolve.cxx
new file mode 100644 (file)
index 0000000..bf821be
--- /dev/null
@@ -0,0 +1,1065 @@
+#include "AliMinResSolve.h"
+#include "AliLog.h"
+#include <TStopwatch.h>
+
+using namespace std;
+
+ClassImp(AliMinResSolve)
+
+//______________________________________________________
+AliMinResSolve::AliMinResSolve() : 
+fSize(0),fPrecon(0),fMatrix(0),fRHS(0),
+  fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
+  fPvv(0),fPvz(0),fPhh(0),
+  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
+{}
+
+//______________________________________________________
+AliMinResSolve::AliMinResSolve(const AliMinResSolve& src) : 
+  TObject(src),
+  fSize(src.fSize),fPrecon(src.fPrecon),fMatrix(src.fMatrix),fRHS(src.fRHS),
+  fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
+  fPvv(0),fPvz(0),fPhh(0),
+  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
+{}
+
+//______________________________________________________
+AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) :
+  fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs->GetMatrixArray()),
+  fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
+  fPvv(0),fPvz(0),fPhh(0),
+  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
+{}
+
+//______________________________________________________
+AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) :
+  fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs),
+  fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
+  fPvv(0),fPvz(0),fPhh(0),
+  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
+{}
+
+//______________________________________________________
+AliMinResSolve::~AliMinResSolve()
+{
+  ClearAux();
+}
+
+//______________________________________________________
+AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
+{
+  if (this != &src) {
+    fSize    = src.fSize;
+    fPrecon  = src.fPrecon;
+    fMatrix  = src.fMatrix;
+    fRHS     = src.fRHS;
+  }
+  return *this;
+}
+
+//_______________________________________________________________
+Int_t AliMinResSolve::BuildPrecon(Int_t prec)
+{
+  const Double_t kTiny = 1E-12;
+  fPrecon = prec;
+  //
+  if (fPrecon == kPreconDiag) return 0;     // nothing to do
+  //
+  else if (fPrecon == kPreconDILU) {
+    fPrecDiag = new Double_t[ fSize ];
+    fPrecAux  = new Double_t[ fSize ];
+    //
+    // copy inverse diagonal
+    for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QuerryDiag(i);
+    //
+    for (int i=0;i<fSize;i++) {
+      fPrecDiag[i] = TMath::Abs(fPrecDiag[i])>kTiny  ? 1./fPrecDiag[i] : 1./TMath::Sign(kTiny,fPrecDiag[i]);
+      for (int j=i+1;j<fSize;j++) {
+       double vl = fMatrix->Querry(j,i);
+       if (vl!=0) fPrecDiag[j] -= fPrecDiag[i]*vl*vl; 
+      }
+    }
+    return 0;
+  } // precon DILU
+  //
+  if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
+    if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(fPrecon-kPreconILU0);
+    else                                          return BuildPreconILUKDense(fPrecon-kPreconILU0);
+  }
+  //
+  return -1;
+}
+
+
+//________________________________ FGMRES METHODS ________________________________
+Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
+{
+  return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov);
+}
+
+//________________________________________________________________________________
+Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
+{
+  // Adapted from Y.Saad fgmrs.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
+  /*----------------------------------------------------------------------
+    |                 *** Preconditioned FGMRES ***                  
+    +-----------------------------------------------------------------------
+    | This is a simple version of the ARMS preconditioned FGMRES algorithm. 
+    +-----------------------------------------------------------------------
+    | Y. S. Dec. 2000. -- Apr. 2008  
+    +-----------------------------------------------------------------------
+    | VecSol  = real vector of length n containing an initial guess to the
+    | precon  = precondtioner id (0 = no precon)
+    | itnlim  = max n of iterations
+    | rtol     = tolerance for stopping iteration
+    | nkrylov = N of Krylov vectors to store
+    +---------------------------------------------------------------------*/
+  int l;
+  double status = kTRUE;
+  double t,beta,eps1=0;
+  const double epsmac    = 2.22E-16;
+  //
+  AliInfo(Form("Solution by FGMRes: Preconditioner#%d Max.iter.: %d Tol.: %e NKrylov: %d",
+              precon,itnlim,rtol,nkrylov));
+  //
+  int its = 0;
+  if (nkrylov<10) {
+    AliInfo(Form("Changing N Krylov vectors from %d to %d",nkrylov,10));
+    nkrylov = 10;
+  }
+  //
+  if (precon>0) {
+    if (precon>=kPreconsTot) {
+      AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon));
+    }
+    else {
+      if  (BuildPrecon(precon)<0) {
+       ClearAux();
+       AliInfo("FGMRES failed to build the preconditioner");
+       return kFALSE;
+      }
+    }
+  }
+  //
+  if (!InitAuxFGMRES(nkrylov)) return kFALSE;
+  //
+  for (l=fSize;l--;) VecSol[l] = 0;
+  //
+  //-------------------- outer loop starts here
+  TStopwatch timer; timer.Start();
+  while(1) {
+    //
+    //-------------------- compute initial residual vector
+    fMatrix->MultiplyByVec(VecSol,fPvv[0]);
+    for (l=fSize;l--;) fPvv[0][l] = fRHS[l] - fPvv[0][l];    //  fPvv[0]= initial residual
+    beta = 0;
+    for (l=fSize;l--;) beta += fPvv[0][l]*fPvv[0][l];
+    beta = TMath::Sqrt(beta);
+    //
+    if (beta == 0.0) break;                                      // success? 
+    t = 1.0 / beta;
+    //--------------------   normalize:  fPvv[0] = fPvv[0] / beta 
+    for (l=fSize;l--;) fPvv[0][l] *= t;
+    if (its == 0) eps1 = rtol*beta;
+    //
+    //    ** initialize 1-st term  of rhs of hessenberg system
+    fPVecV[0] = beta;
+    int i = -1;
+    do {
+      i++;
+      its++;
+      int i1 = i+1; 
+      //
+      //  (Right) Preconditioning Operation   z_{j} = M^{-1} v_{j}
+      //
+      if (precon>0) ApplyPrecon( fPvv[i], fPvz[i]);
+      else          for (l=fSize;l--;)  fPvz[i][l] = fPvv[i][l];
+      //
+      //-------------------- matvec operation w = A z_{j} = A M^{-1} v_{j}
+      fMatrix->MultiplyByVec(fPvz[i],fPvv[i1]);
+      //
+      // modified gram - schmidt...
+      // h_{i,j} = (w,v_{i})
+      // w  = w - h_{i,j} v_{i}
+      //
+      for (int j=0; j<=i; j++) {
+       for (t=0, l=fSize;l--;) t+= fPvv[j][l]*fPvv[i1][l];
+       fPhh[i][j] = t;
+       for (l=fSize;l--;) fPvv[i1][l] -= t*fPvv[j][l];
+      }
+      // -------------------- h_{j+1,j} = ||w||_{2}
+      for (t=0,l=fSize;l--;) t += fPvv[i1][l]*fPvv[i1][l]; t = TMath::Sqrt(t); 
+      fPhh[i][i1] = t;
+      if (t != 0.0) for (t=1./t, l=0; l<fSize; l++) fPvv[i1][l] *= t;  //  v_{j+1} = w / h_{j+1,j}
+      //
+      // done with modified gram schimdt and arnoldi step
+      // now  update factorization of fPhh
+      //
+      // perform previous transformations  on i-th column of h
+      //
+      for (l=1; l<=i; l++) {
+       int l1 = l-1;
+       t = fPhh[i][l1];
+       fPhh[i][l1] = fPVecR1[l1]*t + fPVecR2[l1]*fPhh[i][l];
+       fPhh[i][l] = -fPVecR2[l1]*t + fPVecR1[l1]*fPhh[i][l];
+      }
+      double gam = TMath::Sqrt( fPhh[i][i]*fPhh[i][i] + fPhh[i][i1]*fPhh[i][i1]);
+      //
+      // if gamma is zero then any small value will do...
+      // will affect only residual estimate
+      if (gam == 0.0) gam = epsmac;
+      //  get  next plane rotation
+      fPVecR1[i] = fPhh[i][i]/gam;
+      fPVecR2[i]  = fPhh[i][i1]/gam;
+      fPVecV[i1]  =-fPVecR2[i]*fPVecV[i];
+      fPVecV[i]  *= fPVecR1[i];
+      //
+      //  determine residual norm and test for convergence
+      fPhh[i][i] = fPVecR1[i]*fPhh[i][i] + fPVecR2[i]*fPhh[i][i1];
+      beta = TMath::Abs(fPVecV[i1]);
+      //
+    } while ( (i < nkrylov-1) && (beta > eps1) && (its < itnlim) );
+    //
+    // now compute solution. 1st, solve upper triangular system
+    fPVecV[i] = fPVecV[i]/fPhh[i][i];
+    for (int j=1; j<=i; j++) {
+      int k=i-j;
+      for (t=fPVecV[k],l=k+1; l<=i; l++) t -= fPhh[l][k]*fPVecV[l];
+      fPVecV[k] = t/fPhh[k][k];
+    }
+    // --------------------  linear combination of v[i]'s to get sol. 
+    for (int j=0; j<=i; j++) for (t=fPVecV[j],l=0; l<fSize; l++) VecSol[l] += t*fPvz[j][l];
+    //
+    // --------------------  restart outer loop if needed
+    //    
+    if (beta <= eps1) {
+      timer.Stop();
+      AliInfo(Form("FGMRES converged in %d iterations, CPU time: %.1f sec",its,timer.CpuTime()));
+      break;   // success
+    }
+    //
+    if (its >= itnlim) {
+      timer.Stop();
+      AliInfo(Form("%d iterations limit exceeded, CPU time: %.1f sec",itnlim,timer.CpuTime()));
+      status = kFALSE;
+      break;
+    }
+  }
+  //
+  ClearAux();
+  return status;
+  //
+}
+
+
+//________________________________ MINRES METHODS ________________________________
+Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
+{
+  return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
+}
+
+//________________________________________________________________________________
+Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double rtol)
+{
+  /*
+    Adapted from author's Fortran code:
+    Michael A. Saunders           na.msaunders@na-net.ornl.gov
+    
+    MINRES is an implementation of the algorithm described in the following reference:
+    C. C. Paige and M. A. Saunders (1975),
+    Solution of sparse indefinite systems of linear equations,
+    SIAM J. Numer. Anal. 12(4), pp. 617-629.  
+
+  */
+  if (!fMatrix->IsSymmetric()) {
+    AliInfo("MinRes cannot solve asymmetric matrices, use FGMRes instead");
+    return kFALSE;
+  }
+  //
+  ClearAux();
+  const double eps    = 2.22E-16;
+  double beta1;
+  //
+  if (precon>0) {
+    if (precon>=kPreconsTot) {
+      AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon));
+    }
+    else {
+      if  (BuildPrecon(precon)<0) {
+       ClearAux();
+       AliInfo("MinRes failed to build the preconditioner");
+       return kFALSE;
+      }
+    }
+  }
+  AliInfo(Form("Solution by MinRes: Preconditioner#%d Max.iter.: %d Tol.: %e",precon,itnlim,rtol));
+  //
+  // ------------------------ initialization  ---------------------->>>>
+  memset(VecSol,0,fSize*sizeof(double));
+  int status=0, itn=0;
+  double Anorm = 0;
+  double Acond = 0;
+  double ynorm = 0;
+  double rnorm = 0;
+  double gam,gmax=1,gmin=1,gbar,oldeps,epsa,epsx,epsr,diag, delta,phi,denom,z;
+  //
+  if (!InitAuxMinRes()) return kFALSE;
+  //
+  memset(VecSol,0,fSize*sizeof(double));
+  //
+  // ------------ init aux -------------------------<<<<  
+  //   Set up y and v for the first Lanczos vector v1.
+  //   y  =  beta1 P' v1,  where  P = C**(-1). v is really P' v1.
+  //
+  for (int i=fSize;i--;) fPVecY[i] = fPVecR1[i] = fRHS[i];
+  //
+  if ( precon>0 ) ApplyPrecon( fRHS, fPVecY);
+  beta1 = 0; for (int i=fSize;i--;) beta1 += fRHS[i]*fPVecY[i]; //
+  //
+  if (beta1 < 0) {
+    AliInfo(Form("Preconditioner is indefinite (init) (%e).",beta1));
+    ClearAux();
+    status = 7;
+    return kFALSE;
+  }
+  //
+  if (beta1 == 0) {
+    AliInfo(Form("RHS is zero or is the nullspace of the Preconditioner: Solution is {0}"));
+    ClearAux();
+    return kTRUE;
+  }  
+  //
+  beta1  = TMath::Sqrt( beta1 );       // Normalize y to get v1 later.
+  //
+  //      See if Msolve is symmetric. //RS: Skept
+  //      See if Aprod  is symmetric. //RS: Skept
+  //
+  double oldb   = 0;
+  double beta   = beta1;
+  double dbar   = 0;
+  double epsln  = 0;
+  double qrnorm = beta1;
+  double phibar = beta1;
+  double rhs1   = beta1;
+  double rhs2   = 0;
+  double tnorm2 = 0;
+  double ynorm2 = 0;
+  double cs     = -1;
+  double sn     = 0;
+  for (int i=fSize;i--;) fPVecR2[i] = fPVecR1[i];
+  //
+  TStopwatch timer; timer.Start();
+  while(status==0) { //-----------------  Main iteration loop ---------------------->>>>
+    //
+    itn++;
+    /*-----------------------------------------------------------------
+      Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
+      The general iteration is similar to the case k = 1 with v0 = 0:      
+      p1      = Operator * v1  -  beta1 * v0,
+      alpha1  = v1'p1,
+      q2      = p2  -  alpha1 * v1,
+      beta2^2 = q2'q2,
+      v2      = (1/beta2) q2.      
+      Again, y = betak P vk,  where  P = C**(-1).
+      .... more description needed.
+      -----------------------------------------------------------------*/
+    //
+    double s = 1./beta;                            // Normalize previous vector (in y).
+    for (int i=fSize;i--;) fPVecV[i] = s*fPVecY[i];    // v = vk if P = I
+    //
+    fMatrix->MultiplyByVec(fPVecV,fPVecY);   //      APROD (VecV, VecY);
+    //
+    if (itn>=2) {
+      double btrat = beta/oldb;
+      for (int i=fSize;i--;) fPVecY[i] -= btrat*fPVecR1[i];
+    }
+    double alfa = 0; for (int i=fSize;i--;) alfa += fPVecV[i]*fPVecY[i];    //      alphak
+    //
+    double alf2bt = alfa/beta;
+    for (int i=fSize;i--;) {
+      fPVecY[i] -= alf2bt*fPVecR2[i];
+      fPVecR1[i] = fPVecR2[i];
+      fPVecR2[i] = fPVecY[i];
+    }
+    //
+    if ( precon>0 ) ApplyPrecon(fPVecR2, fPVecY);
+    //
+    oldb  = beta;               //      oldb = betak
+    beta = 0; for (int i=fSize;i--;) beta += fPVecR2[i]*fPVecY[i];    // beta = betak+1^2
+    //
+    if (beta<0) {
+      AliInfo(Form("Preconditioner is indefinite (%e)",beta));
+      status = 7;
+      break;
+    }
+    // 
+    beta   = TMath::Sqrt(beta); //            beta = betak+1
+    tnorm2 += alfa*alfa + oldb*oldb + beta*beta;
+    //
+    if (itn == 1) {             //     Initialize a few things.
+      if (beta/beta1 <= 10.0*eps) {
+       status = 0; //-1   //?????  beta2 = 0 or ~ 0,  terminate later.
+       AliInfo("RHS is eigenvector");
+      }
+      //        !tnorm2 = alfa**2
+      gmax   = TMath::Abs(alfa);    //              alpha1
+      gmin   = gmax;                //              alpha1
+    }
+    //
+    /*
+      Apply previous rotation Qk-1 to get
+      [deltak epslnk+1] = [cs  sn][dbark    0   ]
+      [gbar k dbar k+1]   [sn -cs][alfak betak+1].
+    */
+    //
+    oldeps = epsln;
+    delta  = cs * dbar  +  sn * alfa;       //  delta1 = 0         deltak
+    gbar   = sn * dbar  -  cs * alfa;       //  gbar 1 = alfa1     gbar k
+    epsln  =               sn * beta;       //  epsln2 = 0         epslnk+1
+    dbar   =            -  cs * beta;       //  dbar 2 = beta2     dbar k+1
+    //
+    // Compute the next plane rotation Qk
+    //
+    gam    = TMath::Sqrt( gbar*gbar + beta*beta );  // gammak
+    cs     = gbar / gam;                            // ck
+    sn     = beta / gam;                            // sk
+    phi    = cs * phibar;                           // phik
+    phibar = sn * phibar;                           // phibark+1
+    //
+    // Update  x.
+    denom = 1./gam;
+    //
+    for (int i=fSize;i--;) {
+      fPVecW1[i] = fPVecW2[i];
+      fPVecW2[i] = fPVecW[i];
+      fPVecW[i]  = denom*(fPVecV[i]-oldeps*fPVecW1[i]-delta*fPVecW2[i]);
+      VecSol[i] += phi*fPVecW[i];
+    }
+    //
+    //  Go round again.
+    //
+    gmax = TMath::Max( gmax, gam );
+    gmin = TMath::Min( gmin, gam );
+    z    = rhs1 / gam;
+    ynorm2 += z*z;
+    rhs1   = rhs2  -  delta * z;
+    rhs2   =       -  epsln * z;
+    //
+    //   Estimate various norms and test for convergence.
+    Anorm  = TMath::Sqrt( tnorm2 );
+    ynorm  = TMath::Sqrt( ynorm2 );
+    epsa   = Anorm * eps;
+    epsx   = Anorm * ynorm * eps;
+    epsr   = Anorm * ynorm * rtol;
+    diag   = gbar;
+    if (diag == 0) diag = epsa;
+    //
+    qrnorm = phibar;
+    rnorm  = qrnorm;
+    /*
+      Estimate  cond(A).
+      In this version we look at the diagonals of  R  in the
+      factorization of the lower Hessenberg matrix,  Q * H = R,
+      where H is the tridiagonal matrix from Lanczos with one
+      extra row, beta(k+1) e_k^T.
+    */
+    Acond  = gmax / gmin;
+    //
+    // See if any of the stopping criteria are satisfied.
+    // In rare cases, istop is already -1 from above (Abar = const*I).
+    //
+    if (status == 0) {
+      if (itn    >= itnlim    ) {status = 5; AliInfo(Form("%d iterations limit exceeded",itnlim));}
+      if (Acond  >= 0.1/eps   ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",Acond,0.1/eps));}
+      if (epsx   >= beta1     ) {status = 3; AliInfo(Form("Approximate convergence"));}
+      if (qrnorm <= epsx      ) {status = 2; AliInfo(Form("Converged within machine precision"));}
+      if (qrnorm <= epsr      ) {status = 1; AliInfo(Form("Converged"));}
+    }
+    //
+  }  //-----------------  Main iteration loop ----------------------<<<
+  //
+  ClearAux();
+  //
+  timer.Stop();
+  AliInfo(Form("Exit from MinRes: CPU time: %.2f sec\n"
+              "Status    :  %2d\n"
+              "Iterations:  %4d\n"
+              "Norm      :  %+e\n"
+              "Condition :  %+e\n"
+              "Res.Norm  :  %+e\n"
+              "Sol.Norm  :  %+e",
+              timer.CpuTime(),status,itn,Anorm,Acond,rnorm,ynorm));
+  //
+  return status>=0 && status<=3;
+  //
+}
+
+//______________________________________________________________
+void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
+{
+  ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
+}
+
+//______________________________________________________________
+void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
+{
+  // Application of the preconditioner matrix:
+  // implicitly defines the matrix solving the M*VecOut = VecRHS 
+  const Double_t kTiny = 1E-12;
+  if (fPrecon==kPreconDiag) {
+    for (int i=fSize;i--;) {
+      double d = fMatrix->QuerryDiag(i);
+      vecOut[i] =  vecRHS[i]/(TMath::Abs(d)>kTiny ? d : kTiny);
+    }
+    //    return;
+  }
+  //
+  else if (fPrecon==kPreconDILU) {
+    for (int i=0;i<fSize;i++) {
+      double el = 0;
+      for (int j=i;j--;) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += fPrecAux[j]*vl;}
+      fPrecAux[i] = fPrecDiag[i]*(vecRHS[i]-el); 
+    }
+    for (int i=fSize;i--;) {
+      double el = 0;
+      for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += vl*vecOut[j];}
+      vecOut[i] = fPrecAux[i] - fPrecDiag[i]*el;
+    }
+    //    return;
+  }
+  //
+  else if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
+    //
+    for(int i=0; i<fSize; i++) {     // Block L solve
+      vecOut[i] = vecRHS[i];
+      AliVectorSparse &rowLi = *fMatL->GetRow(i);
+      int n = rowLi.GetNElems();
+      for(int j=0;j<n;j++) vecOut[i] -= vecOut[ rowLi.GetIndex(j) ] * rowLi.GetElem(j);
+
+    }
+    //
+    for(int i=fSize; i--; ) {    // Block -- U solve
+      AliVectorSparse &rowUi = *fMatU->GetRow(i);
+      int n = rowUi.GetNElems();
+      for(int j=0;j<n;j++ ) vecOut[i] -= vecOut[ rowUi.GetIndex(j) ] * rowUi.GetElem(j);
+      vecOut[i] *= fPrecDiag[i];
+    }
+    //
+  }
+  //
+}
+
+
+//___________________________________________________________
+Bool_t AliMinResSolve::InitAuxMinRes()
+{
+  try {
+    fPVecY   = new double[fSize];   
+    fPVecR1  = new double[fSize];   
+    fPVecR2  = new double[fSize];   
+    fPVecV   = new double[fSize];   
+    fPVecW   = new double[fSize];   
+    fPVecW1  = new double[fSize];   
+    fPVecW2  = new double[fSize];   
+  }
+  catch(bad_alloc&) {
+    AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize));
+    ClearAux();
+    return kFALSE;
+  }
+  //
+  for (int i=fSize;i--;) fPVecY[i]=fPVecR1[i]=fPVecR2[i]=fPVecV[i]=fPVecW[i]=fPVecW1[i]=fPVecW2[i]=0.0;
+  //
+  return kTRUE;
+}
+
+
+//___________________________________________________________
+Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
+{
+  try {
+    fPvv     = new double*[nkrylov+1];
+    fPvz     = new double*[nkrylov];
+    for (int i=0; i<=nkrylov; i++) fPvv[i] = new double[fSize];
+    fPhh     = new double*[nkrylov];
+    for (int i=0; i<nkrylov; i++) {
+      fPhh[i] = new double[i+2];
+      fPvz[i]  = new double[fSize];
+    }
+    //
+    fPVecR1  = new double[nkrylov];   
+    fPVecR2  = new double[nkrylov];   
+    fPVecV   = new double[nkrylov+1];
+  }
+  catch(bad_alloc&) {
+    AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize));
+    ClearAux();
+    return kFALSE;
+  }
+  //
+  return kTRUE;
+}
+
+
+//___________________________________________________________
+void AliMinResSolve::ClearAux()
+{
+  if (fPVecY)      delete[] fPVecY;  fPVecY   = 0;
+  if (fPVecR1)     delete[] fPVecR1; fPVecR1  = 0; 
+  if (fPVecR2)     delete[] fPVecR2; fPVecR2  = 0; 
+  if (fPVecV)      delete[] fPVecV;  fPVecV   = 0;
+  if (fPVecW)      delete[] fPVecW;  fPVecW   = 0;
+  if (fPVecW1)     delete[] fPVecW1; fPVecW1  = 0;
+  if (fPVecW2)     delete[] fPVecW2; fPVecW2  = 0;
+  if (fPrecDiag)   delete[] fPrecDiag; fPrecDiag = 0;
+  if (fPrecAux)    delete[] fPrecAux; fPrecAux = 0;
+  if (fMatL)       delete fMatL; fMatL = 0;
+  if (fMatU)       delete fMatU; fMatU = 0;
+}
+
+//___________________________________________________________
+Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
+{
+  /*----------------------------------------------------------------------------
+   * ILUK preconditioner
+   * incomplete LU factorization with level of fill dropping
+   * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
+   *----------------------------------------------------------------------------*/
+  //
+  AliInfo(Form("Building ILU%d preconditioner",lofM));
+  AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
+  //
+  TStopwatch sw; sw.Start();
+  fMatL = new AliMatrixSparse(fSize);
+  fMatU = new AliMatrixSparse(fSize);
+  fMatL->SetSymmetric(kFALSE);
+  fMatU->SetSymmetric(kFALSE);
+  fPrecDiag = new Double_t[fSize];
+  //
+  // symbolic factorization to calculate level of fill index arrays
+  if ( PreconILUKsymb(lofM)<0 ) {
+    ClearAux();
+    return -1;
+  }
+  //
+  Int_t *jw = new Int_t[fSize]; 
+  for(int j=fSize;j--;) jw[j] = -1;     // set indicator array jw to -1 
+  //
+  for(int i=0; i<fSize; i++ ) {            // beginning of main loop
+    /* setup array jw[], and initial i-th row */
+    AliVectorSparse& rowLi = *fMatL->GetRow(i);
+    AliVectorSparse& rowUi = *fMatU->GetRow(i);
+    AliVectorSparse& rowM  = *Matrix->GetRow(i);
+    //
+    for(int j=rowLi.GetNElems(); j--;) {  // initialize L part
+      int col = rowLi.GetIndex(j);
+      jw[col] = j;
+      rowLi.GetElem(j) = 0.;   // do we need this ?
+    }
+    jw[i] = i;
+    fPrecDiag[i] = 0;      // initialize diagonal
+    //
+    for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
+      int col = rowUi.GetIndex(j);
+      jw[col] = j;
+      rowUi.GetElem(j) = 0;
+    }
+    // copy row from csmat into L,U D
+    for(int j=rowM.GetNElems(); j--;) {  // L and D part 
+      if (rowM.GetElem(j)==0) continue;
+      int col = rowM.GetIndex(j);         // (the original matrix stores only lower triangle)
+      if( col < i )   rowLi.GetElem(jw[col]) = rowM.GetElem(j); 
+      else if(col==i) fPrecDiag[i] = rowM.GetElem(j);
+      else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
+    }
+    if (fMatrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) {      // part of the row I on the right of diagonal is stored as 
+       double vl = fMatrix->Querry(col,i);    // the lower part of the column I
+       if (vl==0) continue;                   
+       rowUi.GetElem(jw[col]) = vl;
+      }
+    //
+    // eliminate previous rows
+    for(int j=0; j<rowLi.GetNElems(); j++) {
+      int jrow = rowLi.GetIndex(j);
+      // get the multiplier for row to be eliminated (jrow)
+      rowLi.GetElem(j) *= fPrecDiag[jrow];
+      //
+      // combine current row and row jrow
+      AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
+      for(int k=0; k<rowUj.GetNElems(); k++ ) {
+       int col = rowUj.GetIndex(k);
+       int jpos = jw[col];
+       if( jpos == -1 ) continue;
+       if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
+      }
+    }
+    // reset double-pointer to -1 ( U-part) 
+    for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
+    jw[i] = -1;
+    for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
+    //
+    if( fPrecDiag[i] == 0 ) {
+      AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
+      delete[] jw;
+      return -1;
+    }
+    fPrecDiag[i] = 1.0 / fPrecDiag[i];
+  }
+  //
+  delete[] jw;
+  //
+  sw.Stop();
+  AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
+  AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
+  //
+  return 0;
+}
+
+//___________________________________________________________
+Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
+{
+  /*----------------------------------------------------------------------------
+   * ILUK preconditioner
+   * incomplete LU factorization with level of fill dropping
+   * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
+   *----------------------------------------------------------------------------*/
+  //
+  TStopwatch sw; sw.Start();
+  AliInfo(Form("Building ILU%d preconditioner for dense matrix",lofM));
+  //
+  fMatL = new AliMatrixSparse(fSize);
+  fMatU = new AliMatrixSparse(fSize);
+  fMatL->SetSymmetric(kFALSE);
+  fMatU->SetSymmetric(kFALSE);
+  fPrecDiag = new Double_t[fSize];
+  //
+  // symbolic factorization to calculate level of fill index arrays
+  if ( PreconILUKsymbDense(lofM)<0 ) {
+    ClearAux();
+    return -1;
+  }
+  //
+  Int_t *jw = new Int_t[fSize]; 
+  for(int j=fSize;j--;) jw[j] = -1;     // set indicator array jw to -1 
+  //
+  for(int i=0; i<fSize; i++ ) {            // beginning of main loop
+    /* setup array jw[], and initial i-th row */
+    AliVectorSparse& rowLi = *fMatL->GetRow(i);
+    AliVectorSparse& rowUi = *fMatU->GetRow(i);
+    //
+    for(int j=rowLi.GetNElems(); j--;) {  // initialize L part
+      int col = rowLi.GetIndex(j);
+      jw[col] = j;
+      rowLi.GetElem(j) = 0.;   // do we need this ?
+    }
+    jw[i] = i;
+    fPrecDiag[i] = 0;      // initialize diagonal
+    //
+    for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
+      int col = rowUi.GetIndex(j);
+      jw[col] = j;
+      rowUi.GetElem(j) = 0;
+    }
+    // copy row from csmat into L,U D
+    for(int j=fSize; j--;) {  // L and D part 
+      double vl = fMatrix->Querry(i,j);
+      if (vl==0) continue;
+      if( j < i )   rowLi.GetElem(jw[j]) = vl;
+      else if(j==i) fPrecDiag[i] = vl;
+      else rowUi.GetElem(jw[j]) = vl;
+    }
+    // eliminate previous rows
+    for(int j=0; j<rowLi.GetNElems(); j++) {
+      int jrow = rowLi.GetIndex(j);
+      // get the multiplier for row to be eliminated (jrow)
+      rowLi.GetElem(j) *= fPrecDiag[jrow];
+      //
+      // combine current row and row jrow
+      AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
+      for(int k=0; k<rowUj.GetNElems(); k++ ) {
+       int col = rowUj.GetIndex(k);
+       int jpos = jw[col];
+       if( jpos == -1 ) continue;
+       if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
+      }
+    }
+    // reset double-pointer to -1 ( U-part) 
+    for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
+    jw[i] = -1;
+    for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
+    //
+    if( fPrecDiag[i] == 0 ) {
+      AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
+      delete[] jw;
+      return -1;
+    }
+    fPrecDiag[i] = 1.0 / fPrecDiag[i];
+  }
+  //
+  delete[] jw;
+  //
+  sw.Stop();
+  AliInfo(Form("ILU%d dense preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
+  //  AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
+  //
+  return 0;
+}
+
+//___________________________________________________________
+Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
+{
+  /*----------------------------------------------------------------------------
+   * ILUK preconditioner
+   * incomplete LU factorization with level of fill dropping
+   * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
+   *----------------------------------------------------------------------------*/
+  //
+  AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
+  //
+  UChar_t **ulvl=0,*levls=0;
+  UShort_t *jbuf=0;
+  Int_t    *iw=0;
+  try {
+    ulvl = new UChar_t*[fSize];      // stores lev-fils for U part of ILU factorization
+    levls = new UChar_t[fSize];
+    jbuf = new UShort_t[fSize];
+    iw = new Int_t[fSize];
+  }
+  //
+  catch(bad_alloc&) {
+    AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb");
+    if (ulvl) delete[] ulvl;
+    if (levls) delete[] levls;
+    if (jbuf) delete[] jbuf;
+    if (iw) delete[] iw;
+    ClearAux();
+    return -1;
+  }
+  //
+  for(int j=fSize; j--;) iw[j] = -1;           // initialize iw 
+  for(int i=0; i<fSize; i++) {
+    int incl = 0;
+    int incu = i; 
+    AliVectorSparse& row = *Matrix->GetRow(i);
+    //
+    // assign lof = 0 for matrix elements
+    for(int j=0;j<row.GetNElems(); j++) {
+      int col = row.GetIndex(j);
+      if (row.GetElem(j)==0) continue;  // !!!! matrix is sparse but sometimes 0 appears 
+      if (col<i) {                      // L-part
+       jbuf[incl] = col;
+       levls[incl] = 0;
+       iw[col] = incl++;
+      }
+      else if (col>i) {                 // This works only for general matrix
+       jbuf[incu] = col;
+       levls[incu] = 0;
+       iw[col] = incu++;
+      }
+    }
+    if (Matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
+       if (fMatrix->Querry(col,i)==0) continue;    // Due to the symmetry  == matrix(i,col)
+       jbuf[incu] = col;
+       levls[incu] = 0;
+       iw[col] = incu++;
+      }
+    //
+    // symbolic k,i,j Gaussian elimination
+    int jpiv = -1; 
+    while (++jpiv < incl) {
+      int k = jbuf[jpiv] ;                        // select leftmost pivot
+      int kmin = k;
+      int jmin = jpiv; 
+      for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
+      //
+      // ------------------------------------  swap
+      if(jmin!=jpiv) {
+       jbuf[jpiv] = kmin; 
+       jbuf[jmin] = k; 
+       iw[kmin] = jpiv;
+       iw[k] = jmin; 
+       int tj = levls[jpiv] ;
+       levls[jpiv] = levls[jmin];
+       levls[jmin] = tj;
+       k = kmin; 
+      }
+      // ------------------------------------ symbolic linear combinaiton of rows
+      AliVectorSparse& rowU = *fMatU->GetRow(k);
+      for(int j=0; j<rowU.GetNElems(); j++ ) {
+       int col = rowU.GetIndex(j);
+       int it  = ulvl[k][j]+levls[jpiv]+1; 
+       if( it > lofM ) continue; 
+       int ip = iw[col];
+       if( ip == -1 ) {
+         if( col < i) {
+           jbuf[incl] = col;
+           levls[incl] = it;
+           iw[col] = incl++;
+         } 
+         else if( col > i ) {
+           jbuf[incu] = col;
+           levls[incu] = it;
+           iw[col] = incu++;
+         } 
+       }
+       else levls[ip] = TMath::Min(levls[ip], it); 
+      }
+      //
+    } // end - while loop
+    //
+    // reset iw
+    for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
+    for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
+    //
+    // copy L-part
+    AliVectorSparse& rowLi = *fMatL->GetRow(i);
+    rowLi.ReSize(incl);
+    if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
+    // copy U-part
+    int k = incu-i; 
+    AliVectorSparse& rowUi = *fMatU->GetRow(i);
+    rowUi.ReSize(k);
+    if( k > 0 ) {
+      memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
+      ulvl[i] = new UChar_t[k];   // update matrix of levels 
+      memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
+    }
+  }
+  //
+  // free temp space and leave
+  delete[] levls;
+  delete[] jbuf;
+  for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; 
+  delete[] ulvl; 
+  //
+  fMatL->SortIndices();
+  fMatU->SortIndices();
+  return 0;
+}
+
+
+//___________________________________________________________
+Int_t  AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
+{
+  /*----------------------------------------------------------------------------
+   * ILUK preconditioner
+   * incomplete LU factorization with level of fill dropping
+   * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
+   *----------------------------------------------------------------------------*/
+  //
+  UChar_t **ulvl=0,*levls=0;
+  UShort_t *jbuf=0;
+  Int_t    *iw=0;
+  try {
+    ulvl = new UChar_t*[fSize];      // stores lev-fils for U part of ILU factorization
+    levls = new UChar_t[fSize];
+    jbuf = new UShort_t[fSize];
+    iw = new Int_t[fSize];
+  }
+  //
+  catch(bad_alloc&) {
+    AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb");
+    if (ulvl) delete[] ulvl;
+    if (levls) delete[] levls;
+    if (jbuf) delete[] jbuf;
+    if (iw) delete[] iw;
+    ClearAux();
+    return -1;
+  }
+  //
+  for(int j=fSize; j--;) iw[j] = -1;           // initialize iw 
+  for(int i=0; i<fSize; i++) {
+    int incl = 0;
+    int incu = i; 
+    //
+    // assign lof = 0 for matrix elements
+    for(int j=0;j<fSize; j++) {
+      if (fMatrix->Querry(i,j)==0) continue;
+      if (j<i) {                      // L-part
+       jbuf[incl] = j;
+       levls[incl] = 0;
+       iw[j] = incl++;
+      }
+      else if (j>i) {                 // This works only for general matrix
+       jbuf[incu] = j;
+       levls[incu] = 0;
+       iw[j] = incu++;
+      }
+    }
+    //
+    // symbolic k,i,j Gaussian elimination
+    int jpiv = -1; 
+    while (++jpiv < incl) {
+      int k = jbuf[jpiv] ;                        // select leftmost pivot
+      int kmin = k;
+      int jmin = jpiv; 
+      for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
+      //
+      // ------------------------------------  swap
+      if(jmin!=jpiv) {
+       jbuf[jpiv] = kmin; 
+       jbuf[jmin] = k; 
+       iw[kmin] = jpiv;
+       iw[k] = jmin; 
+       int tj = levls[jpiv] ;
+       levls[jpiv] = levls[jmin];
+       levls[jmin] = tj;
+       k = kmin; 
+      }
+      // ------------------------------------ symbolic linear combinaiton of rows
+      AliVectorSparse& rowU = *fMatU->GetRow(k);
+      for(int j=0; j<rowU.GetNElems(); j++ ) {
+       int col = rowU.GetIndex(j);
+       int it  = ulvl[k][j]+levls[jpiv]+1; 
+       if( it > lofM ) continue; 
+       int ip = iw[col];
+       if( ip == -1 ) {
+         if( col < i) {
+           jbuf[incl] = col;
+           levls[incl] = it;
+           iw[col] = incl++;
+         } 
+         else if( col > i ) {
+           jbuf[incu] = col;
+           levls[incu] = it;
+           iw[col] = incu++;
+         } 
+       }
+       else levls[ip] = TMath::Min(levls[ip], it); 
+      }
+      //
+    } // end - while loop
+    //
+    // reset iw
+    for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
+    for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
+    //
+    // copy L-part
+    AliVectorSparse& rowLi = *fMatL->GetRow(i);
+    rowLi.ReSize(incl);
+    if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
+    // copy U-part
+    int k = incu-i; 
+    AliVectorSparse& rowUi = *fMatU->GetRow(i);
+    rowUi.ReSize(k);
+    if( k > 0 ) {
+      memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
+      ulvl[i] = new UChar_t[k];   // update matrix of levels 
+      memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
+    }
+  }
+  //
+  // free temp space and leave
+  delete[] levls;
+  delete[] jbuf;
+  for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; 
+  delete[] ulvl; 
+  //
+  fMatL->SortIndices();
+  fMatU->SortIndices();
+  return 0;
+}
diff --git a/STEER/AliMinResSolve.h b/STEER/AliMinResSolve.h
new file mode 100644 (file)
index 0000000..c9312e9
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef ALIMINRESSOLVE_H
+#define ALIMINRESSOLVE_H
+
+#include <TObject.h>
+#include <TVectorD.h>
+#include <TMath.h>
+#include "AliMatrixSq.h"
+#include "AliMatrixSparse.h"
+class AliLog;
+class TStopwatch;
+
+class AliMinResSolve : public TObject {
+  //
+ public:
+  enum {kPreconDiag=1,kPreconDILU=2,kPreconILU0=100,kPreconILU10=kPreconILU0+10,kPreconsTot};
+  enum {kSolMinRes,kSolFGMRes,kNSolvers};
+ public:
+  AliMinResSolve();
+  AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs);
+  AliMinResSolve(const AliMatrixSq *mat, const double  * rhs);
+  AliMinResSolve(const AliMinResSolve& src);
+  ~AliMinResSolve();
+  AliMinResSolve& operator=(const AliMinResSolve& rhs);
+  //
+  // ---------  MINRES method (for symmetric matrices)
+  Bool_t SolveMinRes(Double_t* VecSol,Int_t precon=0,int itnlim=2000,double rtol=1e-12);
+  Bool_t SolveMinRes(TVectorD &VecSol,Int_t precon=0,int itnlim=2000,double rtol=1e-12);
+  //
+  // ---------  FGMRES method (for general symmetric matrices)
+  Bool_t SolveFGMRES(Double_t* VecSol,Int_t precon=0,int itnlim=2000,double rtol=1e-12, int nkrylov=60);
+  Bool_t SolveFGMRES(TVectorD &VecSol,Int_t precon=0,int itnlim=2000,double rtol=1e-12, int nkrylov=60);  
+  //
+  Bool_t InitAuxMinRes();
+  Bool_t InitAuxFGMRES(int nkrylov);
+  void   ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut)     const;
+  void   ApplyPrecon(const double*   vecRHS, double*   vecOut)     const;
+  //
+  Int_t  BuildPrecon(Int_t val=0);
+  Int_t  GetPrecon()                                               const    {return fPrecon;} 
+  void   ClearAux();
+  //
+  Int_t  BuildPreconILUK(Int_t lofM);
+  Int_t  BuildPreconILUKDense(Int_t lofM);
+  Int_t  PreconILUKsymb(Int_t lofM);
+  Int_t  PreconILUKsymbDense(Int_t lofM);
+  //
+ protected:
+  //
+  Int_t               fSize;                             // dimension of the input matrix
+  Int_t               fPrecon;                           // preconditioner type
+  const AliMatrixSq*  fMatrix;                           // matrix defining the equations
+  const Double_t*     fRHS;                              // right hand side
+  //
+  Double_t            *fPVecY,*fPVecR1,*fPVecR2,*fPVecV,*fPVecW,*fPVecW1,*fPVecW2;// aux MinRes
+  Double_t            **fPvv,**fPvz,**fPhh;
+  Double_t            *fPrecDiag,*fPrecAux; // aux space
+  AliMatrixSparse     *fMatL,*fMatU;
+  //
+  ClassDef(AliMinResSolve,0)
+};
+
+#endif
diff --git a/STEER/AliSymMatrix.cxx b/STEER/AliSymMatrix.cxx
new file mode 100644 (file)
index 0000000..03020cc
--- /dev/null
@@ -0,0 +1,477 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream>
+//
+#include "TClass.h"
+#include "TMath.h"
+#include "AliSymMatrix.h"
+//
+
+using namespace std;
+
+ClassImp(AliSymMatrix)
+
+
+AliSymMatrix* AliSymMatrix::fgBuffer = 0; 
+Int_t         AliSymMatrix::fgCopyCnt = 0; 
+//___________________________________________________________
+AliSymMatrix::AliSymMatrix() 
+: fElems(0),fElemsAdd(0)
+{
+  fSymmetric = kTRUE;
+  fgCopyCnt++;
+}
+
+//___________________________________________________________
+AliSymMatrix::AliSymMatrix(Int_t size)
+  : AliMatrixSq(),fElems(0),fElemsAdd(0)
+{
+  //
+  fNrows = 0;
+  fNrowIndex = fNcols = size;
+  fElems     = new Double_t[fNcols*(fNcols+1)/2];
+  fSymmetric = kTRUE;
+  Reset();
+  fgCopyCnt++;
+  //
+}
+
+//___________________________________________________________
+AliSymMatrix::AliSymMatrix(const AliSymMatrix &src) 
+  : AliMatrixSq(src),fElems(0),fElemsAdd(0)
+{
+  fNrowIndex = fNcols = src.GetSize();
+  fNrows = 0;
+  if (fNcols) {
+    int nmainel = fNcols*(fNcols+1)/2;
+    fElems     = new Double_t[nmainel];
+    nmainel = src.fNcols*(src.fNcols+1)/2;
+    memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
+    if (src.fNrows) { // transfer extra rows to main matrix
+      Double_t *pnt = fElems + nmainel;
+      int ncl = src.fNcols + 1;
+      for (int ir=0;ir<src.fNrows;ir++) {
+       memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
+       pnt += ncl;
+       ncl++; 
+      }
+    }
+  }
+  else fElems = 0;
+  fElemsAdd = 0;
+  fgCopyCnt++;
+  //
+}
+
+//___________________________________________________________
+AliSymMatrix::~AliSymMatrix() 
+{
+  Clear();
+  if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;}
+}
+
+//___________________________________________________________
+AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
+{
+  //
+  if (this != &src) {
+    TObject::operator=(src);
+    if (fNcols!=src.fNcols && fNrows!=src.fNrows) {
+      // recreate the matrix
+      if (fElems) delete[] fElems;
+      for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
+      delete[] fElemsAdd;
+      //
+      fNrowIndex = fNcols = src.GetSize();
+      fNrows = 0;
+      fElems     = new Double_t[fNcols*(fNcols+1)/2];
+      int nmainel = src.fNcols*(src.fNcols+1);
+      memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
+      if (src.fNrows) { // transfer extra rows to main matrix
+       Double_t *pnt = fElems + nmainel*sizeof(Double_t);
+       int ncl = src.fNcols + 1;
+       for (int ir=0;ir<src.fNrows;ir++) {
+         ncl += ir; 
+         memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
+         pnt += ncl*sizeof(Double_t);
+       }
+      }
+      //
+    }
+    else {
+      memcpy(fElems,src.fElems,fNcols*(fNcols+1)/2*sizeof(Double_t));
+      int ncl = fNcols + 1;
+      for (int ir=0;ir<fNrows;ir++) { // dynamic rows
+       ncl += ir; 
+       memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
+      }
+    }
+  }
+  //
+  return *this;
+}
+
+//___________________________________________________________
+void AliSymMatrix::Clear(Option_t*)
+{
+  if (fElems) {delete[] fElems; fElems = 0;}
+  //  
+  if (fElemsAdd) {
+    for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
+    delete[] fElemsAdd;
+    fElemsAdd = 0;
+  }
+  fNrowIndex = fNcols = fNrows = 0;
+  //
+}
+
+//___________________________________________________________
+Float_t AliSymMatrix::GetDensity() const
+{
+  // get fraction of non-zero elements
+  Int_t nel = 0;
+  for (int i=GetSize();i--;) for (int j=i+1;j--;) if (GetEl(i,j)!=0) nel++;
+  return 2.*nel/( (GetSize()+1)*GetSize() );
+}
+
+//___________________________________________________________
+void AliSymMatrix::Print(Option_t* option) const
+{
+  printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
+  TString opt = option; opt.ToLower();
+  if (opt.IsNull()) return;
+  opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|";
+  for (Int_t i=0;i<fNrowIndex;i++) {
+    printf(opt,i);
+    for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j));
+    printf("\n");
+  }
+}
+
+//___________________________________________________________
+void AliSymMatrix::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
+{
+  // fill vecOut by matrix*vecIn
+  // vector should be of the same size as the matrix
+  for (int i=fNrowIndex;i--;) {
+    vecOut[i] = 0.0;
+    for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
+  }
+  //
+}
+
+//___________________________________________________________
+AliSymMatrix* AliSymMatrix::DecomposeChol() 
+{
+  // Return a matrix with Choleski decomposition
+  /*Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
+    consturcts Cholesky decomposition of SYMMETRIC and
+    POSITIVELY-DEFINED matrix a (a=L*Lt)
+    Only upper triangle of the matrix has to be filled.
+    In opposite to function from the book, the matrix is modified:
+    lower triangle and diagonal are refilled. */
+  //
+  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+    delete fgBuffer; 
+    try {
+      fgBuffer = new AliSymMatrix(*this);
+    }
+    catch(bad_alloc&) {
+      printf("Failed to allocate memory for Choleski decompostions\n");
+      return 0;
+    }
+  }
+  else (*fgBuffer) = *this;
+  //
+  AliSymMatrix& mchol = *fgBuffer;
+  //
+  for (int i=0;i<fNrowIndex;i++) {
+    for (int j=i;j<fNrowIndex;j++) {
+      double sum=(*this)(i,j);
+      for (int k=i-1;k>=0;k--) sum -= mchol(i,k)*mchol(j,k);
+      if (i == j) {
+       if (sum <= 0.0) { // not positive-definite
+         printf("The matrix is not positive definite [%e]\n"
+                "Choleski decomposition is not possible\n",sum);
+         return 0;
+       }
+       mchol(i,i) = TMath::Sqrt(sum);
+       //
+      } else mchol(j,i) = sum/mchol(i,i);
+    }
+  }
+  return fgBuffer;
+}
+
+//___________________________________________________________
+Bool_t AliSymMatrix::InvertChol() 
+{
+  // Invert matrix using Choleski decomposition
+  //
+  AliSymMatrix* mchol = DecomposeChol();
+  if (!mchol) {
+    printf("Failed to invert the matrix\n");
+    return kFALSE;
+  }
+  //
+  InvertChol(mchol);
+  return kTRUE;
+  //
+}
+//___________________________________________________________
+void AliSymMatrix::InvertChol(AliSymMatrix* pmchol) 
+{
+  // Invert matrix using Choleski decomposition, provided the Cholseki's L matrix
+  Int_t i,j,k;
+  Double_t sum;
+  AliSymMatrix& mchol = *pmchol;
+  //
+  // Invert decomposed triangular L matrix (Lower triangle is filled)
+  for (i=0;i<fNrowIndex;i++) { 
+    mchol(i,i) =  1.0/mchol(i,i);
+    for (j=i+1;j<fNrowIndex;j++) { 
+      sum = 0.0; 
+      for (k=i;k<j;k++) sum -= mchol(j,k)*mchol(k,i);
+      mchol(j,i) = sum/mchol(j,j);
+    } 
+  }
+  //
+  // take product of the inverted Choleski L matrix with its transposed
+  for (int i=fNrowIndex;i--;) {
+    for (int j=i+1;j--;) {
+      sum = 0;
+      for (int k=i;k<fNrowIndex;k++) sum += mchol(i,k)*mchol(j,k);
+      (*this)(j,i) = sum;
+    }
+  }
+  //
+}
+
+
+//___________________________________________________________
+Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert) 
+{
+  /* Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
+     Solves the set of n linear equations A x = b, 
+     where a is a positive-definite symmetric matrix. 
+     a[1..n][1..n] is the output of the routine CholDecomposw. 
+     Only the lower triangle of a is accessed. b[1..n] is input as the 
+     right-hand side vector. The solution vector is returned in b[1..n].*/
+  //
+  Int_t i,k;
+  Double_t sum;
+  //
+  AliSymMatrix *pmchol = DecomposeChol();
+  if (!pmchol) {
+    printf("SolveChol failed\n");
+    return kFALSE;
+  }
+  AliSymMatrix& mchol = *pmchol;
+  //
+  for (i=0;i<fNrowIndex;i++) {
+    for (sum=b[i],k=i-1;k>=0;k--) sum -= mchol(i,k)*b[k];
+    b[i]=sum/mchol(i,i);
+  }
+  for (i=fNrowIndex-1;i>=0;i--) {
+    for (sum=b[i],k=i+1;k<fNrowIndex;k++) sum -= mchol(k,i)*b[k];
+    b[i]=sum/mchol(i,i);
+  }
+  //
+  if (invert) InvertChol(pmchol);
+  return kTRUE;
+  //
+}
+
+//___________________________________________________________
+Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert) 
+{
+  return SolveChol((Double_t*)b.GetMatrixArray(),invert);
+}
+
+
+//___________________________________________________________
+Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert) 
+{
+  memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
+  return SolveChol(bsol,invert);  
+}
+
+//___________________________________________________________
+Bool_t AliSymMatrix::SolveChol(TVectorD &brhs, TVectorD &bsol,Bool_t invert) 
+{
+  bsol = brhs;
+  return SolveChol(bsol,invert);
+}
+
+//___________________________________________________________
+void AliSymMatrix::AddRows(int nrows)
+{
+  if (nrows<1) return;
+  Double_t **pnew = new Double_t*[nrows+fNrows];
+  for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows
+  for (int ir=0;ir<nrows;ir++) {
+    int ncl = GetSize()+1;
+    pnew[fNrows] = new Double_t[ncl];
+    memset(pnew[fNrows],0,ncl*sizeof(Double_t));
+    fNrows++;
+    fNrowIndex++;
+  }
+  delete[] fElemsAdd;
+  fElemsAdd = pnew;
+  //
+}
+
+//___________________________________________________________
+void AliSymMatrix::Reset()
+{
+  // if additional rows exist, regularize it
+  if (fElemsAdd) {
+    delete[] fElems;
+    for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
+    delete[] fElemsAdd; fElemsAdd = 0;
+    fNcols = fNrowIndex;
+    fElems = new Double_t[fNcols*(fNcols+1)/2];
+    fNrows = 0;
+  }
+  if (fElems) memset(fElems,0,fNcols*(fNcols+1)/2*sizeof(Double_t));
+  //
+}
+
+//___________________________________________________________
+int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
+{
+  //   Solution a la MP1: gaussian eliminations
+  ///  Obtain solution of a system of linear equations with symmetric matrix 
+  ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
+  //
+
+  Int_t nRank = 0;
+  int iPivot;
+  double vPivot = 0.;
+  double eps = 0.00000000000001;
+  int nGlo = GetSize();
+  bool   *bUnUsed = new bool[nGlo];
+  double *rowMax,*colMax=0;
+  rowMax  = new double[nGlo];
+  //
+  if (stabilize) {
+    colMax   = new double[nGlo];
+    for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0;
+    for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { 
+       double vl = TMath::Abs(Querry(i,j));
+       if (vl==0) continue;
+       if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i
+       if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j
+       if (i==j) continue;
+       if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j
+       if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i
+      }
+    //
+    for (Int_t i=nGlo; i--;) {
+      if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
+      if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
+    }
+    //
+  }
+  //
+  for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
+  //  
+  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+    delete fgBuffer; 
+    try {
+      fgBuffer = new AliSymMatrix(*this);
+    }
+    catch(bad_alloc&) {
+      printf("Failed to allocate memory for matrix inversion buffer\n");
+      return 0;
+    }
+  }
+  else (*fgBuffer) = *this;
+  //
+  if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) 
+      for (int j=0;j<=i; j++) {
+       double vl = Querry(i,j);
+       if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+      }
+      for (int j=i+1;j<nGlo;j++) {
+       double vl = Querry(j,i);
+       if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+      }
+    }
+  //
+  for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QuerryDiag(j)); // save diagonal elem absolute values 
+  //
+  for (Int_t i=0; i<nGlo; i++) {
+    vPivot = 0.0;
+    iPivot = -1;
+    //
+    for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element       
+      double vl;
+      if (bUnUsed[j] && (TMath::Abs(vl=QuerryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QuerryDiag(j)))) {    
+       vPivot = vl;
+       iPivot = j;
+      }
+    }
+    //
+    if (iPivot >= 0) {   // pivot found          
+      nRank++;
+      bUnUsed[iPivot] = false; // This value is used
+      vPivot = 1.0/vPivot;
+      DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse
+      //
+      for (Int_t j=0; j<nGlo; j++) {      
+       for (Int_t jj=0; jj<nGlo; jj++) {  
+         if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)         
+           double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j);
+           r -= vPivot* ( j>iPivot  ? Querry(j,iPivot)  : fgBuffer->Querry(iPivot,j) )
+             *          ( iPivot>jj ? Querry(iPivot,jj) : fgBuffer->Querry(jj,iPivot));
+         }
+       }
+      }
+      //
+      for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements 
+         (*this)(j,iPivot)     *= vPivot;
+         (*fgBuffer)(iPivot,j) *= vPivot;
+       }
+      //
+    }
+    else {  // No more pivot value (clear those elements)
+      for (Int_t j=0; j<nGlo; j++) {
+       if (bUnUsed[j]) {
+         vecB[j] = 0.0;
+         for (Int_t k=0; k<nGlo; k++) {
+           (*this)(j,k) = 0.;
+           if (j!=k) (*fgBuffer)(j,k) = 0;
+         }
+       }
+      }
+      break;  // No more pivots anyway, stop here
+    }
+  }
+  //
+  for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) {
+      double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V
+      if (i>=j) (*this)(i,j) *= vl;
+      else      (*fgBuffer)(j,i) *= vl;
+    }
+  //
+  for (Int_t j=0; j<nGlo; j++) {
+    rowMax[j] = 0.0;
+    for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
+      double vl;
+      if (j>=jj) vl = (*this)(j,jj)     = -Querry(j,jj);
+      else       vl = (*fgBuffer)(j,jj) = -fgBuffer->Querry(j,jj);
+      rowMax[j] += vl*vecB[jj];
+    }          
+  }
+  
+  for (Int_t j=0; j<nGlo; j++) {
+    vecB[j] = rowMax[j]; // The final result
+  }
+  //
+  delete [] bUnUsed;
+  delete [] rowMax;
+  if (stabilize) delete [] colMax;
+
+  return nRank;
+}
diff --git a/STEER/AliSymMatrix.h b/STEER/AliSymMatrix.h
new file mode 100644 (file)
index 0000000..b154e32
--- /dev/null
@@ -0,0 +1,115 @@
+#ifndef ALISYMMATRIX_H
+#define ALISYMMATRIX_H
+
+#include <string.h>
+#include <TObject.h>
+#include <TVectorD.h>
+#include "AliMatrixSq.h"
+
+
+
+class AliSymMatrix : public AliMatrixSq {
+  //
+ public:
+  AliSymMatrix();
+  AliSymMatrix(Int_t size);
+  AliSymMatrix(const AliSymMatrix &mat);
+  virtual ~AliSymMatrix();
+  //
+  void          Clear(Option_t* option="");
+  void          Reset();
+  //
+  Int_t         GetSize()                                        const {return fNrowIndex;}
+  Float_t       GetDensity()                                     const;
+  AliSymMatrix& operator=(const AliSymMatrix& src);
+  Double_t      operator()(Int_t rown, Int_t coln)               const;
+  Double_t&     operator()(Int_t rown, Int_t coln);
+  //
+  Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymMatrix*)this)(r,r);}
+  Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
+  //
+  void          Print(Option_t* option="")                       const;
+  void          AddRows(int nrows=1);
+  //
+  void          Scale(Double_t coeff);
+  void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
+  void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
+  //
+  // ---------------------------------- Dummy methods of MatrixBase
+  virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
+  virtual             Double_t   *GetMatrixArray  ()       {return (Double_t*)fElems;}
+  virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
+  virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
+  virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
+  virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
+  virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
+  virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
+  virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
+  virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
+  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
+  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
+  //
+  // ----------------------------- Choleski methods ----------------------------------------
+  AliSymMatrix*       DecomposeChol();                  //Obtain Cholesky decomposition L matrix 
+  void                InvertChol(AliSymMatrix* mchol);  //Invert using provided Choleski decomposition
+  Bool_t              InvertChol();                     //Invert 
+  Bool_t              SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
+  Bool_t              SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert=kFALSE);
+  Bool_t              SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
+  Bool_t              SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
+  //
+  int                 SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);
+
+ protected:
+  virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
+  Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
+  void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
+  //
+ protected:
+  Double_t*  fElems;     //   Elements booked by constructor
+  Double_t** fElemsAdd;  //   Elements (rows) added dynamicaly
+  //
+  static AliSymMatrix* fgBuffer;  // buffer for fast solution
+  static Int_t         fgCopyCnt;  // matrix copy counter
+  ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
+};
+
+
+//___________________________________________________________
+inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
+{
+  // lower triangle is actually filled
+  return ((row*(row+1))>>1) + col;
+}
+
+//___________________________________________________________
+inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
+{
+  //
+  if (row<col) Swap(row,col);
+  if (row>=fNrowIndex) return 0;
+  return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
+}
+
+//___________________________________________________________
+inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
+{
+  if (row<col) Swap(row,col);
+  if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
+  return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
+}
+
+//___________________________________________________________
+inline void AliSymMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
+{
+  MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
+}
+
+//___________________________________________________________
+inline void AliSymMatrix::Scale(Double_t coeff)
+{
+  for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
+}
+
+
+#endif