--- /dev/null
+#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);
+}
--- /dev/null
+#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
+
--- /dev/null
+#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);
+ //
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff