+/**********************************************************************************************/
+/* Abstract class for matrix used for millepede2 operation. */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
+
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
//
#include "TClass.h"
#include "TMath.h"
-#include "TVectorD.h"
#include "AliMatrixSq.h"
//
#ifndef ALIMATRIXSQ_H
#define ALIMATRIXSQ_H
+/**********************************************************************************************/
+/* Abstract class for matrix used for millepede2 operation. */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
#include <TMatrixDBase.h>
#include <TVectorD.h>
#include <TStopwatch.h>
#include <TFile.h>
#include <TMath.h>
+#include <TVectorD.h>
#include "AliMatrixSq.h"
#include "AliSymMatrix.h"
#include "AliRectMatrix.h"
}
double vl;
//
+ Int_t maxLocUsed = 0;
+ //
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 );
// Symmetric matrix, don't bother j>i coeffs
for (int i=nrefLoc[ip];i--;) { // Fill local matrix and vector
fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i];
+ if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i];
for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j];
- }
+ }
//
} // end of the transfer of the measurement record to matrices
//
+ matCLoc.SetSizeUsed(++maxLocUsed); // data with B=0 may use less than declared nLocals
+ //
// 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 requested, store the track params and errors
// printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
- if (localParams) for (int i=fNLocPar; i--;) {
+ if (localParams) for (int i=maxLocUsed; i--;) {
localParams[2*i] = fVecBLoc[i];
localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
}
nEq++; // number of equations
} // end of Calculate residuals
//
- int nDoF = nEq-fNLocPar;
+ int nDoF = nEq-maxLocUsed;
lChi2 = (nDoF>0) ? lChi2/nDoF : 0; // Chi^2/dof
//
if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
int iCIDg = fGlo2CGlo[iIDg]; // compressed Index of index
if (iCIDg == -1) {
Double_t *rowGL = matCGloLoc(nGloInFit);
- for (int k=fNLocPar;k--;) rowGL[k] = 0.0; // reset the row
+ for (int k=maxLocUsed;k--;) rowGL[k] = 0.0; // reset the row
iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
fCGlo2Glo[nGloInFit++] = iIDg;
}
//
vl = 0;
Double_t *rowGLIDg = matCGloLoc(iCIDg);
- for (int kl=0;kl<fNLocPar;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
+ for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
if (vl!=0) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
//
int nfill = 0;
//
vl = 0;
Double_t *rowGLJDg = matCGloLoc(jCIDg);
- for (int kl=0;kl<fNLocPar;kl++) {
+ for (int kl=0;kl<maxLocUsed;kl++) {
// diag terms
if ( (vll=rowGLIDg[kl]*rowGLJDg[kl])!=0 ) vl += matCLoc.QueryDiag(kl)*vll;
//
for (Long_t i=0;i<ndr;i++) {
ReadRecordData(i);
LocalFit();
- if ( (i%int(0.1*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
+ if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
}
swt.Stop();
printf("%ld local fits done: ", ndr);
// add large number to diagonal of fixed params
//
for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
+ // printf("#%3d : Nproc : %5d grp: %d\n",i,fProcPnt[i],fParamGrID[i]);
if (fProcPnt[i]<1) {
fNGloFix++;
fVecBGlo[i] = 0.;
/**********************************************************************************************/
#include <TObject.h>
+#include <TString.h>
#include <TTree.h>
#include "AliMinResSolve.h"
#include "AliMillePedeRecord.h"
+/**********************************************************************************************/
+/* General class for solving large system of linear equations */
+/* Includes MINRES, FGMRES methods as well as a few precondiotiong methods */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
+
#include "AliMinResSolve.h"
#include "AliLog.h"
+#include "AliMatrixSq.h"
+#include "AliMatrixSparse.h"
+#include "AliSymBDMatrix.h"
#include <TStopwatch.h>
+#include <float.h>
+#include <TMath.h>
using namespace std;
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)
-{}
+ fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+ // default constructor
+}
//______________________________________________________
AliMinResSolve::AliMinResSolve(const AliMinResSolve& 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)
-{}
+ fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+ // copy constructor
+}
//______________________________________________________
AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) :
- fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs->GetMatrixArray()),
+ fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)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)
-{}
+ fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+ // copy accepting equation
+}
//______________________________________________________
AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) :
- fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs),
+ fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)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)
-{}
+ fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+ // copy accepting equation
+}
//______________________________________________________
AliMinResSolve::~AliMinResSolve()
{
+ // destructor
ClearAux();
}
//______________________________________________________
AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
{
+ // assignment op.
if (this != &src) {
fSize = src.fSize;
fPrecon = src.fPrecon;
//_______________________________________________________________
Int_t AliMinResSolve::BuildPrecon(Int_t prec)
{
- const Double_t kTiny = 1E-12;
+ // preconditioner building
+ // 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->QueryDiag(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->Query(j,i);
- if (vl!=0) fPrecDiag[j] -= fPrecDiag[i]*vl*vl;
- }
- }
- return 0;
- } // precon DILU
+ if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
+ return BuildPreconBD(fPrecon - kPreconBD); // with halfbandwidth + diagonal = fPrecon
+ }
//
if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(fPrecon-kPreconILU0);
return -1;
}
-
//________________________________ FGMRES METHODS ________________________________
Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
{
+ // solve by fgmres
return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov);
}
for (l=fSize;l--;) beta += fPvv[0][l]*fPvv[0][l];
beta = TMath::Sqrt(beta);
//
- if (beta == 0.0) break; // success?
+ if (beta < epsmac) break; // success?
t = 1.0 / beta;
//-------------------- normalize: fPvv[0] = fPvv[0] / beta
for (l=fSize;l--;) fPvv[0][l] *= t;
// -------------------- 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}
+ if (t > 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
//
// if gamma is zero then any small value will do...
// will affect only residual estimate
- if (gam == 0.0) gam = epsmac;
+ if (gam < epsmac) gam = epsmac;
// get next plane rotation
fPVecR1[i] = fPhh[i][i]/gam;
fPVecR2[i] = fPhh[i][i1]/gam;
//________________________________ MINRES METHODS ________________________________
Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
{
+ // solve by minres
return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
}
// ------------------------ initialization ---------------------->>>>
memset(VecSol,0,fSize*sizeof(double));
int status=0, itn=0;
- double Anorm = 0;
- double Acond = 0;
+ double normA = 0;
+ double condA = 0;
double ynorm = 0;
double rnorm = 0;
double gam,gmax=1,gmin=1,gbar,oldeps,epsa,epsx,epsr,diag, delta,phi,denom,z;
return kFALSE;
}
//
- if (beta1 == 0) {
+ if (beta1 < eps) {
AliInfo(Form("RHS is zero or is the nullspace of the Preconditioner: Solution is {0}"));
ClearAux();
return kTRUE;
rhs2 = - epsln * z;
//
// Estimate various norms and test for convergence.
- Anorm = TMath::Sqrt( tnorm2 );
+ normA = TMath::Sqrt( tnorm2 );
ynorm = TMath::Sqrt( ynorm2 );
- epsa = Anorm * eps;
- epsx = Anorm * ynorm * eps;
- epsr = Anorm * ynorm * rtol;
+ epsa = normA * eps;
+ epsx = normA * ynorm * eps;
+ epsr = normA * ynorm * rtol;
diag = gbar;
if (diag == 0) diag = epsa;
//
where H is the tridiagonal matrix from Lanczos with one
extra row, beta(k+1) e_k^T.
*/
- Acond = gmax / gmin;
+ condA = gmax / gmin;
//
// See if any of the stopping criteria are satisfied.
// In rare cases, istop is already -1 from above (Abar = const*I).
//
+ AliInfo(Form("#%5d |qnrm: %+.2e Anrm:%+.2e Cnd:%+.2e Rnrm:%+.2e Ynrm:%+.2e EpsR:%+.2e EpsX:%+.2e Beta1:%+.2e",
+ itn,qrnorm, normA,condA,rnorm,ynorm,epsr,epsx,beta1));
+
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 (condA >= 0.1/eps ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",condA,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"));}
"Condition : %+e\n"
"Res.Norm : %+e\n"
"Sol.Norm : %+e",
- timer.CpuTime(),status,itn,Anorm,Acond,rnorm,ynorm));
+ timer.CpuTime(),status,itn,normA,condA,rnorm,ynorm));
//
return status>=0 && status<=3;
//
//______________________________________________________________
void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
{
+ // apply precond.
ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
}
{
// 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->QueryDiag(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->Query(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->Query(i,j);if (vl!=0) el += vl*vecOut[j];}
- vecOut[i] = fPrecAux[i] - fPrecDiag[i]*el;
- }
+ // const Double_t kTiny = 1E-12;
+ if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
+ fMatBD->Solve(vecRHS,vecOut);
// return;
}
//
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];
+ vecOut[i] *= fDiagLU[i];
}
//
}
//___________________________________________________________
Bool_t AliMinResSolve::InitAuxMinRes()
{
+ // init auxiliary space for minres
try {
fPVecY = new double[fSize];
fPVecR1 = new double[fSize];
//___________________________________________________________
Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
{
+ // init auxiliary space for fgmres
try {
fPvv = new double*[nkrylov+1];
fPvz = new double*[nkrylov];
//___________________________________________________________
void AliMinResSolve::ClearAux()
{
+ // clear aux. space
if (fPVecY) delete[] fPVecY; fPVecY = 0;
if (fPVecR1) delete[] fPVecR1; fPVecR1 = 0;
if (fPVecR2) delete[] fPVecR2; fPVecR2 = 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 (fDiagLU) delete[] fDiagLU; fDiagLU = 0;
if (fMatL) delete fMatL; fMatL = 0;
if (fMatU) delete fMatU; fMatU = 0;
+ if (fMatBD) delete fMatBD; fMatBD = 0;
+}
+
+//___________________________________________________________
+Int_t AliMinResSolve::BuildPreconBD(Int_t hwidth)
+{
+ // build preconditioner
+ AliInfo(Form("Building Band-Diagonal preconditioner of half-width = %d",hwidth));
+ fMatBD = new AliSymBDMatrix( fMatrix->GetSize(), hwidth );
+ //
+ // fill the band-diagonal part of the matrix
+ if (fMatrix->InheritsFrom("AliMatrixSparse")) {
+ for (int ir=fMatrix->GetSize();ir--;) {
+ int jmin = TMath::Max(0,ir-hwidth);
+ AliVectorSparse& irow = *((AliMatrixSparse*)fMatrix)->GetRow(ir);
+ for(int j=irow.GetNElems();j--;) {
+ int jind = irow.GetIndex(j);
+ if (jind<jmin) break;
+ (*fMatBD)(ir,jind) = irow.GetElem(j);
+ }
+ }
+ }
+ else {
+ for (int ir=fMatrix->GetSize();ir--;) {
+ int jmin = TMath::Max(0,ir-hwidth);
+ for(int jr=jmin;jr<=ir;jr++) (*fMatBD)(ir,jr) = fMatrix->Query(ir,jr);
+ }
+ }
+ //
+ fMatBD->DecomposeLDLT();
+ //
+ return 0;
}
//___________________________________________________________
*----------------------------------------------------------------------------*/
//
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];
+ fDiagLU = new Double_t[fSize];
+ AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
//
// symbolic factorization to calculate level of fill index arrays
if ( PreconILUKsymb(lofM)<0 ) {
//
for(int i=0; i<fSize; i++ ) { // beginning of main loop
if ( (i%int(0.1*fSize)) == 0) {
- printf("BuildPrecon: row %d of %d\n",i,fSize);
+ AliInfo(Form("BuildPrecon: row %d of %d",i,fSize));
sw.Stop();
sw.Print();
sw.Start(kFALSE);
/* setup array jw[], and initial i-th row */
AliVectorSparse& rowLi = *fMatL->GetRow(i);
AliVectorSparse& rowUi = *fMatU->GetRow(i);
- AliVectorSparse& rowM = *Matrix->GetRow(i);
+ AliVectorSparse& rowM = *matrix->GetRow(i);
//
for(int j=rowLi.GetNElems(); j--;) { // initialize L part
int col = rowLi.GetIndex(j);
rowLi.GetElem(j) = 0.; // do we need this ?
}
jw[i] = i;
- fPrecDiag[i] = 0; // initialize diagonal
+ fDiagLU[i] = 0; // initialize diagonal
//
for(int j=rowUi.GetNElems();j--;) { // initialize U part
int col = rowUi.GetIndex(j);
}
// copy row from csmat into L,U D
for(int j=rowM.GetNElems(); j--;) { // L and D part
- if (rowM.GetElem(j)==0) continue;
+ if (TMath::Abs(rowM.GetElem(j))<DBL_MIN) 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 if(col==i) fDiagLU[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->Query(col,i); // the lower part of the column I
- if (vl==0) continue;
+ if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // part of the row I on the right of diagonal is stored as
+ double vl = matrix->Query(col,i); // the lower part of the column I
+ if (TMath::Abs(vl)<DBL_MIN) continue;
rowUi.GetElem(jw[col]) = vl;
}
//
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];
+ rowLi.GetElem(j) *= fDiagLU[jrow];
//
// combine current row and row jrow
AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
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 if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
}
}
jw[i] = -1;
for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
//
- if( fPrecDiag[i] == 0 ) {
+ if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
delete[] jw;
return -1;
}
- fPrecDiag[i] = 1.0 / fPrecDiag[i];
+ fDiagLU[i] = 1.0 / fDiagLU[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()));
+ AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
//
return 0;
}
fMatU = new AliMatrixSparse(fSize);
fMatL->SetSymmetric(kFALSE);
fMatU->SetSymmetric(kFALSE);
- fPrecDiag = new Double_t[fSize];
+ fDiagLU = new Double_t[fSize];
//
// symbolic factorization to calculate level of fill index arrays
if ( PreconILUKsymbDense(lofM)<0 ) {
rowLi.GetElem(j) = 0.; // do we need this ?
}
jw[i] = i;
- fPrecDiag[i] = 0; // initialize diagonal
+ fDiagLU[i] = 0; // initialize diagonal
//
for(int j=rowUi.GetNElems();j--;) { // initialize U part
int col = rowUi.GetIndex(j);
// copy row from csmat into L,U D
for(int j=fSize; j--;) { // L and D part
double vl = fMatrix->Query(i,j);
- if (vl==0) continue;
+ if (TMath::Abs(vl)<DBL_MIN) continue;
if( j < i ) rowLi.GetElem(jw[j]) = vl;
- else if(j==i) fPrecDiag[i] = vl;
+ else if(j==i) fDiagLU[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];
+ rowLi.GetElem(j) *= fDiagLU[jrow];
//
// combine current row and row jrow
AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
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 if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
}
}
jw[i] = -1;
for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
//
- if( fPrecDiag[i] == 0 ) {
+ if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
delete[] jw;
return -1;
}
- fPrecDiag[i] = 1.0 / fPrecDiag[i];
+ fDiagLU[i] = 1.0 / fDiagLU[i];
}
//
delete[] jw;
*----------------------------------------------------------------------------*/
//
TStopwatch sw;
- printf("PreconILUKsymb>>\n");
+ AliInfo("PreconILUKsymb>>");
+ AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
sw.Start();
- AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
//
UChar_t **ulvl=0,*levls=0;
UShort_t *jbuf=0;
for(int i=0; i<fSize; i++) {
int incl = 0;
int incu = i;
- AliVectorSparse& row = *Matrix->GetRow(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 (TMath::Abs(row.GetElem(j))<DBL_MIN) continue; // !!!! matrix is sparse but sometimes 0 appears
if (col<i) { // L-part
jbuf[incl] = col;
levls[incl] = 0;
iw[col] = incu++;
}
}
- if (Matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
- if (fMatrix->Query(col,i)==0) continue; // Due to the symmetry == matrix(i,col)
+ if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
+ if (TMath::Abs(matrix->Query(col,i))<DBL_MIN) continue; // Due to the symmetry == matrix(i,col)
jbuf[incu] = col;
levls[incu] = 0;
iw[col] = incu++;
fMatU->SortIndices();
sw.Stop();
sw.Print();
- printf("PreconILUKsymb<<\n");
+ AliInfo("PreconILUKsymb<<");
return 0;
}
//
// assign lof = 0 for matrix elements
for(int j=0;j<fSize; j++) {
- if (fMatrix->Query(i,j)==0) continue;
+ if (TMath::Abs(fMatrix->Query(i,j))<DBL_MIN) continue;
if (j<i) { // L-part
jbuf[incl] = j;
levls[incl] = 0;
#ifndef ALIMINRESSOLVE_H
#define ALIMINRESSOLVE_H
+/**********************************************************************************************/
+/* General class for solving large system of linear equations */
+/* Includes MINRES, FGMRES methods as well as a few precondiotiong methods */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
+
#include <TObject.h>
#include <TVectorD.h>
-#include <TMath.h>
-#include "AliMatrixSq.h"
-#include "AliMatrixSparse.h"
-class AliLog;
-class TStopwatch;
+class AliMatrixSq;
+class AliMatrixSparse;
+class AliSymBDMatrix;
+
class AliMinResSolve : public TObject {
//
public:
- enum {kPreconDiag=1,kPreconDILU=2,kPreconILU0=100,kPreconILU10=kPreconILU0+10,kPreconsTot};
+ enum {kPreconBD=1,kPreconILU0=100,kPreconILU10=kPreconILU0+10,kPreconsTot};
enum {kSolMinRes,kSolFGMRes,kNSolvers};
public:
AliMinResSolve();
Int_t GetPrecon() const {return fPrecon;}
void ClearAux();
//
+ Int_t BuildPreconBD(Int_t hwidth);
Int_t BuildPreconILUK(Int_t lofM);
Int_t BuildPreconILUKDense(Int_t lofM);
Int_t PreconILUKsymb(Int_t lofM);
//
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
+ AliMatrixSq* fMatrix; // matrix defining the equations
+ 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;
+ Double_t *fPVecY; // aux. space
+ Double_t *fPVecR1; // aux. space
+ Double_t *fPVecR2; // aux. space
+ Double_t *fPVecV; // aux. space
+ Double_t *fPVecW; // aux. space
+ Double_t *fPVecW1; // aux. space
+ Double_t *fPVecW2; // aux. space
+ Double_t **fPvv; // aux. space
+ Double_t **fPvz; // aux. space
+ Double_t **fPhh; // aux. space
+ Double_t *fDiagLU; // aux space
+ AliMatrixSparse *fMatL; // aux. space
+ AliMatrixSparse *fMatU; // aux. space
+ AliSymBDMatrix *fMatBD; // aux. space
//
ClassDef(AliMinResSolve,0)
};
--- /dev/null
+
+/*********************************************************************************/
+/* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
+/* Only lower triangle is stored in the "profile" format */
+/* */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/*********************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream>
+#include <float.h>
+//
+#include "TClass.h"
+#include "TMath.h"
+#include "AliSymBDMatrix.h"
+//
+
+using namespace std;
+
+ClassImp(AliSymBDMatrix)
+
+
+//___________________________________________________________
+AliSymBDMatrix::AliSymBDMatrix()
+: fElems(0)
+{
+ // def. c-tor
+ fSymmetric = kTRUE;
+}
+
+//___________________________________________________________
+AliSymBDMatrix::AliSymBDMatrix(Int_t size, Int_t w)
+ : AliMatrixSq(),fElems(0)
+{
+ // c-tor for given size
+ //
+ fNcols = size; // number of rows
+ if (w<0) w = 0;
+ if (w>=size) w = size-1;
+ fNrows = w;
+ fRowLwb = w+1;
+ fSymmetric = kTRUE;
+ //
+ // total number of stored elements
+ fNelems = size*(w+1) - w*(w+1)/2;
+ //
+ fElems = new Double_t[fNcols*fRowLwb];
+ memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
+ //
+}
+
+//___________________________________________________________
+AliSymBDMatrix::AliSymBDMatrix(const AliSymBDMatrix &src)
+ : AliMatrixSq(src),fElems(0)
+{
+ // copy c-tor
+ if (src.GetSize()<1) return;
+ fNcols = src.GetSize();
+ fNrows = src.GetBandHWidth();
+ fRowLwb = fNrows+1;
+ fNelems = src.GetNElemsStored();
+ fElems = new Double_t[fNcols*fRowLwb];
+ memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
+ //
+}
+
+//___________________________________________________________
+AliSymBDMatrix::~AliSymBDMatrix()
+{
+ // d-tor
+ Clear();
+}
+
+//___________________________________________________________
+AliSymBDMatrix& AliSymBDMatrix::operator=(const AliSymBDMatrix& src)
+{
+ // assignment
+ //
+ if (this != &src) {
+ TObject::operator=(src);
+ if (fNcols!=src.fNcols) {
+ // recreate the matrix
+ if (fElems) delete[] fElems;
+ fNcols = src.GetSize();
+ fNrows = src.GetBandHWidth();
+ fNelems = src.GetNElemsStored();
+ fRowLwb = src.fRowLwb;
+ fElems = new Double_t[fNcols*fRowLwb];
+ }
+ memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
+ fSymmetric = kTRUE;
+ }
+ //
+ return *this;
+ //
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::Clear(Option_t*)
+{
+ // clear dynamic part
+ if (fElems) {delete[] fElems; fElems = 0;}
+ //
+ fNelems = fNcols = fNrows = fRowLwb = 0;
+ //
+}
+
+//___________________________________________________________
+Float_t AliSymBDMatrix::GetDensity() const
+{
+ // get fraction of non-zero elements
+ if (!fNelems) return 0;
+ Int_t nel = 0;
+ for (int i=fNelems;i--;) if (TMath::Abs(fElems[i])>DBL_MIN) nel++;
+ return nel/fNelems;
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::Print(Option_t* option) const
+{
+ // print data
+ printf("Symmetric Band-Diagonal Matrix : Size = %d, half bandwidth = %d\n",
+ GetSize(),GetBandHWidth());
+ 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<GetSize();i++) {
+ printf(opt,i);
+ for (Int_t j=TMath::Max(0,i-GetBandHWidth());j<=i;j++) printf("%+.3e|",GetEl(i,j));
+ printf("\n");
+ }
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
+{
+ // fill vecOut by matrix*vecIn
+ // vector should be of the same size as the matrix
+ if (IsDecomposed()) {
+ for (int i=0;i<GetSize();i++) {
+ double sm = 0;
+ int jmax = TMath::Min(GetSize(),i+fRowLwb);
+ for (int j=i+1;j<jmax;j++) sm += vecIn[j]*Query(j,i);
+ vecOut[i] = QueryDiag(i)*(vecIn[i]+sm);
+ }
+ for (int i=GetSize();i--;) {
+ double sm = 0;
+ int jmin = TMath::Max(0,i - GetBandHWidth());
+ int jmax = i-1;
+ for (int j=jmin;j<jmax;j++) sm += vecOut[j]*Query(i,j);
+ vecOut[i] += sm;
+ }
+ }
+ else { // not decomposed
+ for (int i=GetSize();i--;) {
+ vecOut[i] = 0.0;
+ int jmin = TMath::Max(0,i - GetBandHWidth());
+ int jmax = TMath::Min(GetSize(),i + fRowLwb);
+ for (int j=jmin;j<jmax;j++) vecOut[i] += vecIn[j]*Query(i,j);
+ }
+ }
+ //
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::Reset()
+{
+ // set all elems to 0
+ if (fElems) memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
+ //
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
+{
+ // add list of elements to row r
+ for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::DecomposeLDLT()
+{
+ // decomposition to L Diag L^T
+ if (IsDecomposed()) return;
+ //
+ Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
+ //
+ Double_t dtmp,gamma=0.0,xi=0.0;
+ int iDiag;
+ //
+ // find max diag and number of non-0 diag.elements
+ for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
+ if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
+ if (gamma < dtmp) gamma = dtmp;
+ }
+ //
+ // find max. off-diag element
+ for (int ir=1;ir<iDiag;ir++) {
+ for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
+ if (ic<0) continue;
+ dtmp = TMath::Abs(Query(ir,ic));
+ if (xi<dtmp) xi = dtmp;
+ }
+ }
+ double delta = eps*TMath::Max(1.0,xi+gamma);
+ //
+ double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
+ double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn)));
+ //
+ for (int kr=1;kr<GetSize();kr++) {
+ int colKmin = TMath::Max(0,kr - GetBandHWidth());
+ double theta = 0.0;
+ //
+ for (int jr=colKmin;jr<=kr;jr++) {
+ int colJmin = TMath::Max(0,jr - GetBandHWidth());
+ //
+ dtmp = 0.0;
+ for (int i=TMath::Max(colKmin,colJmin);i<jr;i++)
+ dtmp += Query(kr,i)*QueryDiag(i)*Query(jr,i);
+ dtmp = (*this)(kr,jr) -= dtmp;
+ //
+ theta = TMath::Max(theta, TMath::Abs(dtmp));
+ //
+ if (jr!=kr) {
+ if (TMath::Abs(QueryDiag(jr))>DBL_MIN) (*this)(kr,jr) /= QueryDiag(jr);
+ else (*this)(kr,jr) = 0.0;
+ }
+ else if (kr<iDiag) {
+ dtmp = theta/beta;
+ dtmp *= dtmp;
+ dtmp = TMath::Max(dtmp, delta);
+ (*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp);
+ }
+ } // jr
+ } // kr
+ //
+ for (int i=0;i<GetSize();i++) {
+ dtmp = QueryDiag(i);
+ if (TMath::Abs(dtmp)>DBL_MIN) DiagElem(i) = 1./dtmp;
+ }
+ //
+ SetDecomposed();
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::Solve(Double_t *rhs)
+{
+ // solve matrix equation
+ //
+ if (!IsDecomposed()) DecomposeLDLT();
+ //
+ for (int kr=0;kr<GetSize();kr++)
+ for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
+ //
+ for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
+ //
+ for (int kr=GetSize();kr--;)
+ for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
+ //
+}
+
+//___________________________________________________________
+void AliSymBDMatrix::Solve(const Double_t *rhs,Double_t *sol)
+{
+ // solve matrix equation
+ memcpy(sol,rhs,GetSize()*sizeof(Double_t));
+ Solve(sol);
+ //
+}
--- /dev/null
+#ifndef ALISYMBDMATRIX_H
+#define ALISYMBDMATRIX_H
+
+/*********************************************************************************/
+/* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
+/* Only lower triangle is stored in the "profile" format */
+/* */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/*********************************************************************************/
+
+//#include <string.h>
+#include <TObject.h>
+#include <TVectorD.h>
+#include "AliMatrixSq.h"
+
+
+class AliSymBDMatrix : public AliMatrixSq {
+ //
+ public:
+ enum {kDecomposedBit = 0x1};
+ //
+ AliSymBDMatrix();
+ AliSymBDMatrix(Int_t size, Int_t w = 0);
+ AliSymBDMatrix(const AliSymBDMatrix &mat);
+ virtual ~AliSymBDMatrix();
+ //
+ Int_t GetBandHWidth() const {return fNrows;}
+ Int_t GetNElemsStored() const {return fNelems;}
+ void Clear(Option_t* option="");
+ void Reset();
+ //
+ Float_t GetDensity() const;
+ AliSymBDMatrix& operator=(const AliSymBDMatrix& src);
+ Double_t operator()(Int_t rown, Int_t coln) const;
+ Double_t& operator()(Int_t rown, Int_t coln);
+ Double_t operator()(Int_t rown) const;
+ Double_t& operator()(Int_t rown);
+ //
+ Double_t DiagElem(Int_t r) const {return (*(const AliSymBDMatrix*)this)(r,r);}
+ Double_t& DiagElem(Int_t r) {return (*this)(r,r);}
+ void DecomposeLDLT();
+ void Solve(Double_t *rhs);
+ void Solve(const Double_t *rhs,Double_t *sol);
+ void Solve(TVectorD &rhs) {Solve(rhs.GetMatrixArray());}
+ void Solve(const TVectorD &rhs,TVectorD &sol) {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());}
+ //
+ void Print(Option_t* option="") const;
+ void SetDecomposed() {SetBit(kDecomposedBit);}
+ Bool_t IsDecomposed() const {return TestBit(kDecomposedBit);}
+ //
+ void MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
+ void MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
+ void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
+ //
+ // protected:
+ virtual Int_t GetIndex(Int_t row,Int_t col) const;
+ virtual Int_t GetIndex(Int_t diagID) 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
+ //
+ ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
+};
+
+
+//___________________________________________________________
+inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
+{
+ // lower triangle band is actually filled
+ if (row<col) Swap(row,col);
+ col -= row;
+ if (col < -GetBandHWidth()) return -1;
+ return GetIndex(row) + col;
+}
+
+//___________________________________________________________
+inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
+{
+ // Get index of the diagonal element on row diagID
+ return (diagID+1)*fRowLwb-1;
+}
+
+//___________________________________________________________
+inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
+{
+ // query element
+ int idx = GetIndex(row,col);
+ return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
+}
+
+//___________________________________________________________
+inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
+{
+ // get element for assingment; assignment outside of the stored range has no effect
+ int idx = GetIndex(row,col);
+ if (idx>=0) return fElems[idx];
+ fTol = 0;
+ return fTol;
+}
+
+//___________________________________________________________
+inline Double_t AliSymBDMatrix::operator()(Int_t row) const
+{
+ // query diagonal
+ return (const Double_t&) fElems[GetIndex(row)];
+}
+
+//___________________________________________________________
+inline Double_t& AliSymBDMatrix::operator()(Int_t row)
+{
+ // get diagonal for assingment; assignment outside of the stored range has no effect
+ return fElems[GetIndex(row)];
+}
+
+//___________________________________________________________
+inline void AliSymBDMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
+{
+ MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
+}
+
+#endif
+/**********************************************************************************************/
+/* Fast symmetric matrix with dynamically expandable size. */
+/* Only part can be used for matrix operations. It is defined as: */
+/* fNCols: rows built by constructor (GetSizeBooked) */
+/* fNRows: number of rows added dynamically (automatically added on assignment to row) */
+/* GetNRowAdded */
+/* fNRowIndex: total size (fNCols+fNRows), GetSize */
+/* fRowLwb : actual size to used for given operation, by default = total size, GetSizeUsed */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
+#include <float.h>
//
-#include "TClass.h"
-#include "TMath.h"
+#include <TClass.h>
+#include <TMath.h>
#include "AliSymMatrix.h"
+#include "AliLog.h"
//
using namespace std;
{
//
fNrows = 0;
- fNrowIndex = fNcols = size;
+ fNrowIndex = fNcols = fRowLwb = size;
fElems = new Double_t[fNcols*(fNcols+1)/2];
fSymmetric = kTRUE;
Reset();
{
fNrowIndex = fNcols = src.GetSize();
fNrows = 0;
+ fRowLwb = src.GetSizeUsed();
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
+ if (src.GetSizeAdded()) { // transfer extra rows to main matrix
Double_t *pnt = fElems + nmainel;
- int ncl = src.fNcols + 1;
- for (int ir=0;ir<src.fNrows;ir++) {
+ int ncl = src.GetSizeBooked() + 1;
+ for (int ir=0;ir<src.GetSizeAdded();ir++) {
memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
pnt += ncl;
ncl++;
//
if (this != &src) {
TObject::operator=(src);
- if (fNcols!=src.fNcols && fNrows!=src.fNrows) {
+ if (GetSizeBooked()!=src.GetSizeBooked() && GetSizeAdded()!=src.GetSizeAdded()) {
// recreate the matrix
if (fElems) delete[] fElems;
- for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
+ for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i];
delete[] fElemsAdd;
//
fNrowIndex = src.GetSize();
fNcols = src.GetSize();
fNrows = 0;
- fElems = new Double_t[fNcols*(fNcols+1)/2];
- int nmainel = src.fNcols*(src.fNcols+1);
+ fRowLwb = src.GetSizeUsed();
+ fElems = new Double_t[GetSize()*(GetSize()+1)/2];
+ int nmainel = src.GetSizeBooked()*(src.GetSizeBooked()+1);
memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
- if (src.fNrows) { // transfer extra rows to main matrix
+ if (src.GetSizeAdded()) { // 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++) {
+ int ncl = src.GetSizeBooked() + 1;
+ for (int ir=0;ir<src.GetSizeAdded();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
+ memcpy(fElems,src.fElems,GetSizeBooked()*(GetSizeBooked()+1)/2*sizeof(Double_t));
+ int ncl = GetSizeBooked() + 1;
+ for (int ir=0;ir<GetSizeAdded();ir++) { // dynamic rows
ncl += ir;
memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
}
return *this;
}
+//___________________________________________________________
+AliSymMatrix& AliSymMatrix::operator+=(const AliSymMatrix& src)
+{
+ //
+ if (GetSizeUsed() != src.GetSizeUsed()) {
+ AliError("Matrix sizes are different");
+ return *this;
+ }
+ for (int i=0;i<GetSizeUsed();i++) for (int j=i;j<GetSizeUsed();j++) (*this)(j,i) += src(j,i);
+ 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];
+ for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i];
delete[] fElemsAdd;
fElemsAdd = 0;
}
- fNrowIndex = 0;
- fNcols = 0;
- fNrows = 0;
+ fNrowIndex = fNcols = fNrows = fRowLwb = 0;
//
}
{
// 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() );
+ for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (TMath::Abs(GetEl(i,j))>DBL_MIN) nel++;
+ return 2.*nel/( (GetSizeUsed()+1)*GetSizeUsed() );
}
//___________________________________________________________
void AliSymMatrix::Print(Option_t* option) const
{
- printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
+ printf("Symmetric Matrix: Size = %d (%d rows added dynamically), %d used\n",GetSize(),GetSizeAdded(),GetSizeUsed());
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++) {
+ for (Int_t i=0;i<GetSizeUsed();i++) {
printf(opt,i);
for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j));
printf("\n");
{
// fill vecOut by matrix*vecIn
// vector should be of the same size as the matrix
- for (int i=fNrowIndex;i--;) {
+ for (int i=GetSizeUsed();i--;) {
vecOut[i] = 0.0;
- for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
+ for (int j=GetSizeUsed();j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
}
//
}
// In opposite to function from the book, the matrix is modified:
// lower triangle and diagonal are refilled.
//
- if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+ if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) {
delete fgBuffer;
try {
fgBuffer = new AliSymMatrix(*this);
//
AliSymMatrix& mchol = *fgBuffer;
//
- for (int i=0;i<fNrowIndex;i++) {
+ for (int i=0;i<GetSizeUsed();i++) {
Double_t *rowi = mchol.GetRow(i);
- for (int j=i;j<fNrowIndex;j++) {
+ for (int j=i;j<GetSizeUsed();j++) {
Double_t *rowj = mchol.GetRow(j);
double sum = rowj[i];
for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k];
AliSymMatrix& mchol = *pmchol;
//
// Invert decomposed triangular L matrix (Lower triangle is filled)
- for (int i=0;i<fNrowIndex;i++) {
+ for (int i=0;i<GetSizeUsed();i++) {
mchol(i,i) = 1.0/mchol(i,i);
- for (int j=i+1;j<fNrowIndex;j++) {
+ for (int j=i+1;j<GetSizeUsed();j++) {
Double_t *rowj = mchol.GetRow(j);
sum = 0.0;
for (int k=i;k<j;k++) if (rowj[k]) {
}
//
// take product of the inverted Choleski L matrix with its transposed
- for (int i=fNrowIndex;i--;) {
+ for (int i=GetSizeUsed();i--;) {
for (int j=i+1;j--;) {
sum = 0;
- for (int k=i;k<fNrowIndex;k++) {
+ for (int k=i;k<GetSizeUsed();k++) {
double &mik = mchol(i,k);
if (mik) {
double &mjk = mchol(j,k);
}
AliSymMatrix& mchol = *pmchol;
//
- for (i=0;i<fNrowIndex;i++) {
+ for (i=0;i<GetSizeUsed();i++) {
Double_t *rowi = mchol.GetRow(i);
for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k];
b[i]=sum/rowi[i];
}
//
- for (i=fNrowIndex-1;i>=0;i--) {
- for (sum=b[i],k=i+1;k<fNrowIndex;k++) if (b[k]) {
+ for (i=GetSizeUsed()-1;i>=0;i--) {
+ for (sum=b[i],k=i+1;k<GetSizeUsed();k++) if (b[k]) {
double &mki=mchol(k,i); if (mki) sum -= mki*b[k];
}
b[i]=sum/mchol(i,i);
//___________________________________________________________
Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert)
{
- memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
+ memcpy(bsol,brhs,GetSizeUsed()*sizeof(Double_t));
return SolveChol(bsol,invert);
}
memset(pnew[fNrows],0,ncl*sizeof(Double_t));
fNrows++;
fNrowIndex++;
+ fRowLwb++;
}
delete[] fElemsAdd;
fElemsAdd = pnew;
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];
+ fNcols = fRowLwb = fNrowIndex;
+ fElems = new Double_t[GetSize()*(GetSize()+1)/2];
fNrows = 0;
}
- if (fElems) memset(fElems,0,fNcols*(fNcols+1)/2*sizeof(Double_t));
+ if (fElems) memset(fElems,0,GetSize()*(GetSize()+1)/2*sizeof(Double_t));
//
}
//___________________________________________________________
Double_t* AliSymMatrix::GetRow(Int_t r)
{
- if (r>=fNrowIndex) {
- int nn = fNrowIndex;
- AddRows(r-fNrowIndex+1);
+ if (r>=GetSize()) {
+ int nn = GetSize();
+ AddRows(r-GetSize()+1);
printf("create %d of %d\n",r, nn);
- return &((fElemsAdd[r-fNcols])[0]);
+ return &((fElemsAdd[r-GetSizeBooked()])[0]);
}
else return &fElems[GetIndex(r,0)];
}
Int_t nRank = 0;
int iPivot;
double vPivot = 0.;
- double eps = 0.00000000000001;
- int nGlo = GetSize();
+ double eps = 1e-14;
+ int nGlo = GetSizeUsed();
bool *bUnUsed = new bool[nGlo];
double *rowMax,*colMax=0;
rowMax = 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(Query(i,j));
- if (vl==0) continue;
+ if (vl<DBL_MIN) 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;
}
//
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
+ if (TMath::Abs(rowMax[i])>DBL_MIN) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
+ if (TMath::Abs(colMax[i])>DBL_MIN) colMax[i] = 1./colMax[i]; // Max elemt of column i
}
//
}
//
for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
//
- if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+ if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) {
delete fgBuffer;
try {
fgBuffer = new AliSymMatrix(*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 = Query(i,j);
- if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+ if (TMath::Abs(vl)>DBL_MIN) 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 = Query(j,i);
- if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+ if (TMath::Abs(vl)>DBL_MIN) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
}
}
//
#ifndef ALISYMMATRIX_H
#define ALISYMMATRIX_H
+/**********************************************************************************************/
+/* Fast symmetric matrix with dynamically expandable size. */
+/* Only part can be used for matrix operations. It is defined as: */
+/* fNCols: rows built by constructor (GetSizeBooked) */
+/* fNRows: number of rows added dynamically (automatically added on assignment to row) */
+/* GetNRowAdded */
+/* fNRowIndex: total size (fNCols+fNRows), GetSize */
+/* fRowLwb : actual size to used for given operation, by default = total size, GetSizeUsed */
+/* */
+/* Author: ruben.shahoyan@cern.ch */
+/* */
+/**********************************************************************************************/
#include <string.h>
#include <TObject.h>
void Reset();
//
Int_t GetSize() const {return fNrowIndex;}
+ Int_t GetSizeUsed() const {return fRowLwb;}
+ Int_t GetSizeBooked() const {return fNcols;}
+ Int_t GetSizeAdded() const {return fNrows;}
Float_t GetDensity() const;
AliSymMatrix& operator=(const AliSymMatrix& src);
+ AliSymMatrix& operator+=(const AliSymMatrix& src);
Double_t operator()(Int_t rown, Int_t coln) const;
Double_t& operator()(Int_t rown, Int_t coln);
//
//
void Print(Option_t* option="") const;
void AddRows(int nrows=1);
+ void SetSizeUsed(Int_t sz) {fRowLwb = sz;}
//
void Scale(Double_t coeff);
void MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;