2 /*********************************************************************************/
3 /* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
4 /* Only lower triangle is stored in the "profile" format */
7 /* Author: ruben.shahoyan@cern.ch */
9 /*********************************************************************************/
17 #include "AliSymBDMatrix.h"
22 ClassImp(AliSymBDMatrix)
25 //___________________________________________________________
26 AliSymBDMatrix::AliSymBDMatrix()
33 //___________________________________________________________
34 AliSymBDMatrix::AliSymBDMatrix(Int_t size, Int_t w)
35 : AliMatrixSq(),fElems(0)
37 // c-tor for given size
39 fNcols = size; // number of rows
41 if (w>=size) w = size-1;
46 // total number of stored elements
47 fNelems = size*(w+1) - w*(w+1)/2;
49 fElems = new Double_t[fNcols*fRowLwb];
50 memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
54 //___________________________________________________________
55 AliSymBDMatrix::AliSymBDMatrix(const AliSymBDMatrix &src)
56 : AliMatrixSq(src),fElems(0)
59 if (src.GetSize()<1) return;
60 fNcols = src.GetSize();
61 fNrows = src.GetBandHWidth();
63 fNelems = src.GetNElemsStored();
64 fElems = new Double_t[fNcols*fRowLwb];
65 memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
69 //___________________________________________________________
70 AliSymBDMatrix::~AliSymBDMatrix()
76 //___________________________________________________________
77 AliSymBDMatrix& AliSymBDMatrix::operator=(const AliSymBDMatrix& src)
82 TObject::operator=(src);
83 if (fNcols!=src.fNcols) {
84 // recreate the matrix
85 if (fElems) delete[] fElems;
86 fNcols = src.GetSize();
87 fNrows = src.GetBandHWidth();
88 fNelems = src.GetNElemsStored();
89 fRowLwb = src.fRowLwb;
90 fElems = new Double_t[fNcols*fRowLwb];
92 memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
100 //___________________________________________________________
101 void AliSymBDMatrix::Clear(Option_t*)
103 // clear dynamic part
104 if (fElems) {delete[] fElems; fElems = 0;}
106 fNelems = fNcols = fNrows = fRowLwb = 0;
110 //___________________________________________________________
111 Float_t AliSymBDMatrix::GetDensity() const
113 // get fraction of non-zero elements
114 if (!fNelems) return 0;
116 for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
120 //___________________________________________________________
121 void AliSymBDMatrix::Print(Option_t* option) const
124 printf("Symmetric Band-Diagonal Matrix : Size = %d, half bandwidth = %d\n",
125 GetSize(),GetBandHWidth());
126 TString opt = option; opt.ToLower();
127 if (opt.IsNull()) return;
128 opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|";
129 for (Int_t i=0;i<GetSize();i++) {
131 for (Int_t j=TMath::Max(0,i-GetBandHWidth());j<=i;j++) printf("%+.3e|",GetEl(i,j));
136 //___________________________________________________________
137 void AliSymBDMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const
139 // fill vecOut by matrix*vecIn
140 // vector should be of the same size as the matrix
141 if (IsDecomposed()) {
142 for (int i=0;i<GetSize();i++) {
144 int jmax = TMath::Min(GetSize(),i+fRowLwb);
145 for (int j=i+1;j<jmax;j++) sm += vecIn[j]*Query(j,i);
146 vecOut[i] = QueryDiag(i)*(vecIn[i]+sm);
148 for (int i=GetSize();i--;) {
150 int jmin = TMath::Max(0,i - GetBandHWidth());
152 for (int j=jmin;j<jmax;j++) sm += vecOut[j]*Query(i,j);
156 else { // not decomposed
157 for (int i=GetSize();i--;) {
159 int jmin = TMath::Max(0,i - GetBandHWidth());
160 int jmax = TMath::Min(GetSize(),i + fRowLwb);
161 for (int j=jmin;j<jmax;j++) vecOut[i] += vecIn[j]*Query(i,j);
167 //___________________________________________________________
168 void AliSymBDMatrix::Reset()
170 // set all elems to 0
171 if (fElems) memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
172 SetDecomposed(kFALSE);
176 //___________________________________________________________
177 void AliSymBDMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
179 // add list of elements to row r
180 for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
183 //___________________________________________________________
184 void AliSymBDMatrix::DecomposeLDLT()
186 // decomposition to L Diag L^T
187 if (IsDecomposed()) return;
189 Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
191 Double_t dtmp,gamma=0.0,xi=0.0;
194 // find max diag and number of non-0 diag.elements
195 for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
196 if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
197 if (gamma < dtmp) gamma = dtmp;
200 // find max. off-diag element
201 for (int ir=1;ir<iDiag;ir++) {
202 for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
204 dtmp = TMath::Abs(Query(ir,ic));
205 if (xi<dtmp) xi = dtmp;
208 double delta = eps*TMath::Max(1.0,xi+gamma);
210 double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
211 double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn)));
213 for (int kr=1;kr<GetSize();kr++) {
214 int colKmin = TMath::Max(0,kr - GetBandHWidth());
217 for (int jr=colKmin;jr<=kr;jr++) {
218 int colJmin = TMath::Max(0,jr - GetBandHWidth());
221 for (int i=TMath::Max(colKmin,colJmin);i<jr;i++)
222 dtmp += Query(kr,i)*QueryDiag(i)*Query(jr,i);
223 dtmp = (*this)(kr,jr) -= dtmp;
225 theta = TMath::Max(theta, TMath::Abs(dtmp));
228 if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
229 else (*this)(kr,jr) = 0.0;
234 dtmp = TMath::Max(dtmp, delta);
235 (*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp);
240 for (int i=0;i<GetSize();i++) {
242 if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
248 //___________________________________________________________
249 void AliSymBDMatrix::Solve(Double_t *rhs)
251 // solve matrix equation
253 if (!IsDecomposed()) DecomposeLDLT();
255 for (int kr=0;kr<GetSize();kr++)
256 for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
258 for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
260 for (int kr=GetSize();kr--;)
261 for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
265 //___________________________________________________________
266 void AliSymBDMatrix::Solve(const Double_t *rhs,Double_t *sol)
268 // solve matrix equation
269 memcpy(sol,rhs,GetSize()*sizeof(Double_t));