ROOT  6.06/09
Reference Guide
TMatrixDSymEigen.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Dec 2003
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 //////////////////////////////////////////////////////////////////////////
13 // //
14 // TMatrixDSymEigen //
15 // //
16 // Eigenvalues and eigenvectors of a real symmetric matrix. //
17 // //
18 // If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is //
19 // diagonal and the eigenvector matrix V is orthogonal. That is, the //
20 // diagonal values of D are the eigenvalues, and V*V' = I, where I is //
21 // the identity matrix. The columns of V represent the eigenvectors in //
22 // the sense that A*V = V*D. //
23 // //
24 //////////////////////////////////////////////////////////////////////////
25 
26 #include "TMatrixDSymEigen.h"
27 #include "TMath.h"
28 
30 
31 ////////////////////////////////////////////////////////////////////////////////
32 /// Constructor for eigen-problem of symmetric matrix A .
33 
35 {
36  R__ASSERT(a.IsValid());
37 
38  const Int_t nRows = a.GetNrows();
39  const Int_t rowLwb = a.GetRowLwb();
40 
41  fEigenValues.ResizeTo(rowLwb,rowLwb+nRows-1);
42  fEigenVectors.ResizeTo(a);
43 
44  fEigenVectors = a;
45 
46  TVectorD offDiag;
47  Double_t work[kWorkMax];
48  if (nRows > kWorkMax) offDiag.ResizeTo(nRows);
49  else offDiag.Use(nRows,work);
50 
51  // Tridiagonalize matrix
52  MakeTridiagonal(fEigenVectors,fEigenValues,offDiag);
53 
54  // Make eigenvectors and -values
55  MakeEigenVectors(fEigenVectors,fEigenValues,offDiag);
56 }
57 
58 ////////////////////////////////////////////////////////////////////////////////
59 /// Copy constructor
60 
62 {
63  *this = another;
64 }
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and
68 /// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
69 /// Fortran subroutine in EISPACK.
70 
72 {
73  Double_t *pV = v.GetMatrixArray();
74  Double_t *pD = d.GetMatrixArray();
75  Double_t *pE = e.GetMatrixArray();
76 
77  const Int_t n = v.GetNrows();
78 
79  Int_t i,j,k;
80  Int_t off_n1 = (n-1)*n;
81  for (j = 0; j < n; j++)
82  pD[j] = pV[off_n1+j];
83 
84  // Householder reduction to tridiagonal form.
85 
86  for (i = n-1; i > 0; i--) {
87  const Int_t off_i1 = (i-1)*n;
88  const Int_t off_i = i*n;
89 
90  // Scale to avoid under/overflow.
91 
92  Double_t scale = 0.0;
93  Double_t h = 0.0;
94  for (k = 0; k < i; k++)
95  scale = scale+TMath::Abs(pD[k]);
96  if (scale == 0.0) {
97  pE[i] = pD[i-1];
98  for (j = 0; j < i; j++) {
99  const Int_t off_j = j*n;
100  pD[j] = pV[off_i1+j];
101  pV[off_i+j] = 0.0;
102  pV[off_j+i] = 0.0;
103  }
104  } else {
105 
106  // Generate Householder vector.
107 
108  for (k = 0; k < i; k++) {
109  pD[k] /= scale;
110  h += pD[k]*pD[k];
111  }
112  Double_t f = pD[i-1];
113  Double_t g = TMath::Sqrt(h);
114  if (f > 0)
115  g = -g;
116  pE[i] = scale*g;
117  h = h-f*g;
118  pD[i-1] = f-g;
119  for (j = 0; j < i; j++)
120  pE[j] = 0.0;
121 
122  // Apply similarity transformation to remaining columns.
123 
124  for (j = 0; j < i; j++) {
125  const Int_t off_j = j*n;
126  f = pD[j];
127  pV[off_j+i] = f;
128  g = pE[j]+pV[off_j+j]*f;
129  for (k = j+1; k <= i-1; k++) {
130  const Int_t off_k = k*n;
131  g += pV[off_k+j]*pD[k];
132  pE[k] += pV[off_k+j]*f;
133  }
134  pE[j] = g;
135  }
136  f = 0.0;
137  for (j = 0; j < i; j++) {
138  pE[j] /= h;
139  f += pE[j]*pD[j];
140  }
141  Double_t hh = f/(h+h);
142  for (j = 0; j < i; j++)
143  pE[j] -= hh*pD[j];
144  for (j = 0; j < i; j++) {
145  f = pD[j];
146  g = pE[j];
147  for (k = j; k <= i-1; k++) {
148  const Int_t off_k = k*n;
149  pV[off_k+j] -= (f*pE[k]+g*pD[k]);
150  }
151  pD[j] = pV[off_i1+j];
152  pV[off_i+j] = 0.0;
153  }
154  }
155  pD[i] = h;
156  }
157 
158  // Accumulate transformations.
159 
160  for (i = 0; i < n-1; i++) {
161  const Int_t off_i = i*n;
162  pV[off_n1+i] = pV[off_i+i];
163  pV[off_i+i] = 1.0;
164  Double_t h = pD[i+1];
165  if (h != 0.0) {
166  for (k = 0; k <= i; k++) {
167  const Int_t off_k = k*n;
168  pD[k] = pV[off_k+i+1]/h;
169  }
170  for (j = 0; j <= i; j++) {
171  Double_t g = 0.0;
172  for (k = 0; k <= i; k++) {
173  const Int_t off_k = k*n;
174  g += pV[off_k+i+1]*pV[off_k+j];
175  }
176  for (k = 0; k <= i; k++) {
177  const Int_t off_k = k*n;
178  pV[off_k+j] -= g*pD[k];
179  }
180  }
181  }
182  for (k = 0; k <= i; k++) {
183  const Int_t off_k = k*n;
184  pV[off_k+i+1] = 0.0;
185  }
186  }
187  for (j = 0; j < n; j++) {
188  pD[j] = pV[off_n1+j];
189  pV[off_n1+j] = 0.0;
190  }
191  pV[off_n1+n-1] = 1.0;
192  pE[0] = 0.0;
193 }
194 
195 ////////////////////////////////////////////////////////////////////////////////
196 /// Symmetric tridiagonal QL algorithm.
197 /// This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and
198 /// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
199 /// Fortran subroutine in EISPACK.
200 
202 {
203  Double_t *pV = v.GetMatrixArray();
204  Double_t *pD = d.GetMatrixArray();
205  Double_t *pE = e.GetMatrixArray();
206 
207  const Int_t n = v.GetNrows();
208 
209  Int_t i,j,k,l;
210  for (i = 1; i < n; i++)
211  pE[i-1] = pE[i];
212  pE[n-1] = 0.0;
213 
214  Double_t f = 0.0;
215  Double_t tst1 = 0.0;
216  Double_t eps = TMath::Power(2.0,-52.0);
217  for (l = 0; l < n; l++) {
218 
219  // Find small subdiagonal element
220 
221  tst1 = TMath::Max(tst1,TMath::Abs(pD[l])+TMath::Abs(pE[l]));
222  Int_t m = l;
223 
224  // Original while-loop from Java code
225  while (m < n) {
226  if (TMath::Abs(pE[m]) <= eps*tst1) {
227  break;
228  }
229  m++;
230  }
231 
232  // If m == l, pD[l] is an eigenvalue,
233  // otherwise, iterate.
234 
235  if (m > l) {
236  Int_t iter = 0;
237  do {
238  if (iter++ == 30) { // (check iteration count here.)
239  Error("MakeEigenVectors","too many iterations");
240  break;
241  }
242 
243  // Compute implicit shift
244 
245  Double_t g = pD[l];
246  Double_t p = (pD[l+1]-g)/(2.0*pE[l]);
247  Double_t r = TMath::Hypot(p,1.0);
248  if (p < 0)
249  r = -r;
250  pD[l] = pE[l]/(p+r);
251  pD[l+1] = pE[l]*(p+r);
252  Double_t dl1 = pD[l+1];
253  Double_t h = g-pD[l];
254  for (i = l+2; i < n; i++)
255  pD[i] -= h;
256  f = f+h;
257 
258  // Implicit QL transformation.
259 
260  p = pD[m];
261  Double_t c = 1.0;
262  Double_t c2 = c;
263  Double_t c3 = c;
264  Double_t el1 = pE[l+1];
265  Double_t s = 0.0;
266  Double_t s2 = 0.0;
267  for (i = m-1; i >= l; i--) {
268  c3 = c2;
269  c2 = c;
270  s2 = s;
271  g = c*pE[i];
272  h = c*p;
273  r = TMath::Hypot(p,pE[i]);
274  pE[i+1] = s*r;
275  s = pE[i]/r;
276  c = p/r;
277  p = c*pD[i]-s*g;
278  pD[i+1] = h+s*(c*g+s*pD[i]);
279 
280  // Accumulate transformation.
281 
282  for (k = 0; k < n; k++) {
283  const Int_t off_k = k*n;
284  h = pV[off_k+i+1];
285  pV[off_k+i+1] = s*pV[off_k+i]+c*h;
286  pV[off_k+i] = c*pV[off_k+i]-s*h;
287  }
288  }
289  p = -s*s2*c3*el1*pE[l]/dl1;
290  pE[l] = s*p;
291  pD[l] = c*p;
292 
293  // Check for convergence.
294 
295  } while (TMath::Abs(pE[l]) > eps*tst1);
296  }
297  pD[l] = pD[l]+f;
298  pE[l] = 0.0;
299  }
300 
301  // Sort eigenvalues and corresponding vectors.
302 
303  for (i = 0; i < n-1; i++) {
304  k = i;
305  Double_t p = pD[i];
306  for (j = i+1; j < n; j++) {
307  if (pD[j] > p) {
308  k = j;
309  p = pD[j];
310  }
311  }
312  if (k != i) {
313  pD[k] = pD[i];
314  pD[i] = p;
315  for (j = 0; j < n; j++) {
316  const Int_t off_j = j*n;
317  p = pV[off_j+i];
318  pV[off_j+i] = pV[off_j+k];
319  pV[off_j+k] = p;
320  }
321  }
322  }
323 }
324 
325 ////////////////////////////////////////////////////////////////////////////////
326 /// Assignment operator
327 
329 {
330  if (this != &source) {
333  }
334  return *this;
335 }
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:291
static void MakeEigenVectors(TMatrixD &v, TVectorD &d, TVectorD &e)
Symmetric tridiagonal QL algorithm.
TH1 * h
Definition: legend2.C:5
#define R__ASSERT(e)
Definition: TError.h:98
TMatrixDSymEigen & operator=(const TMatrixDSymEigen &source)
Assignment operator.
int Int_t
Definition: RtypesCore.h:41
TArc * a
Definition: textangle.C:12
ClassImp(TMatrixDSymEigen) TMatrixDSymEigen
Constructor for eigen-problem of symmetric matrix A .
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1201
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:346
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
void Error(const char *location, const char *msgfmt,...)
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
Element * GetMatrixArray()
Definition: TVectorT.h:84
ROOT::R::TRInterface & r
Definition: Object.C:4
SVector< double, 2 > v
Definition: Dict.h:5
TMarker * m
Definition: textangle.C:8
TLine * l
Definition: textangle.C:4
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:211
static void MakeTridiagonal(TMatrixD &v, TVectorD &d, TVectorD &e)
This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and Wilkinson, Handbook for Auto.
return c2
Definition: legend2.C:14
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:60
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
return c3
Definition: legend3.C:15
const Int_t n
Definition: legend1.C:16