Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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/** \class TMatrixDSymEigen
13 \ingroup Matrix
14
15 TMatrixDSymEigen
16
17 Eigenvalues and eigenvectors of a real symmetric matrix.
18
19 If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
20 diagonal and the eigenvector matrix V is orthogonal. That is, the
21 diagonal values of D are the eigenvalues, and V*V' = I, where I is
22 the identity matrix. The columns of V represent the eigenvectors in
23 the sense that A*V = V*D.
24*/
25
26#include "TMatrixDSymEigen.h"
27#include "TMath.h"
28
29
30////////////////////////////////////////////////////////////////////////////////
31/// Constructor for eigen-problem of symmetric matrix A .
32
34{
35 R__ASSERT(a.IsValid());
36
37 const Int_t nRows = a.GetNrows();
38 const Int_t rowLwb = a.GetRowLwb();
39
42
44
47 if (nRows > kWorkMax) offDiag.ResizeTo(nRows);
48 else offDiag.Use(nRows,work);
49
50 // Tridiagonalize matrix
52
53 // Make eigenvectors and -values
55}
56
57////////////////////////////////////////////////////////////////////////////////
58/// Copy constructor
59
64
65////////////////////////////////////////////////////////////////////////////////
66/// This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and
67/// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
68/// Fortran subroutine in EISPACK.
69
71{
72 Double_t *pV = v.GetMatrixArray();
73 Double_t *pD = d.GetMatrixArray();
74 Double_t *pE = e.GetMatrixArray();
75
76 const Int_t n = v.GetNrows();
77
78 Int_t i,j,k;
79 Int_t off_n1 = (n-1)*n;
80 for (j = 0; j < n; j++)
81 pD[j] = pV[off_n1+j];
82
83 // Householder reduction to tridiagonal form.
84
85 for (i = n-1; i > 0; i--) {
86 const Int_t off_i1 = (i-1)*n;
87 const Int_t off_i = i*n;
88
89 // Scale to avoid under/overflow.
90
91 Double_t scale = 0.0;
92 Double_t h = 0.0;
93 for (k = 0; k < i; k++)
95 if (scale == 0.0) {
96 pE[i] = pD[i-1];
97 for (j = 0; j < i; j++) {
98 const Int_t off_j = j*n;
99 pD[j] = pV[off_i1+j];
100 pV[off_i+j] = 0.0;
101 pV[off_j+i] = 0.0;
102 }
103 } else {
104
105 // Generate Householder vector.
106
107 for (k = 0; k < i; k++) {
108 pD[k] /= scale;
109 h += pD[k]*pD[k];
110 }
111 Double_t f = pD[i-1];
113 if (f > 0)
114 g = -g;
115 pE[i] = scale*g;
116 h = h-f*g;
117 pD[i-1] = f-g;
118 for (j = 0; j < i; j++)
119 pE[j] = 0.0;
120
121 // Apply similarity transformation to remaining columns.
122
123 for (j = 0; j < i; j++) {
124 const Int_t off_j = j*n;
125 f = pD[j];
126 pV[off_j+i] = f;
127 g = pE[j]+pV[off_j+j]*f;
128 for (k = j+1; k <= i-1; k++) {
129 const Int_t off_k = k*n;
130 g += pV[off_k+j]*pD[k];
131 pE[k] += pV[off_k+j]*f;
132 }
133 pE[j] = g;
134 }
135 f = 0.0;
136 for (j = 0; j < i; j++) {
137 pE[j] /= h;
138 f += pE[j]*pD[j];
139 }
140 Double_t hh = f/(h+h);
141 for (j = 0; j < i; j++)
142 pE[j] -= hh*pD[j];
143 for (j = 0; j < i; j++) {
144 f = pD[j];
145 g = pE[j];
146 for (k = j; k <= i-1; k++) {
147 const Int_t off_k = k*n;
148 pV[off_k+j] -= (f*pE[k]+g*pD[k]);
149 }
150 pD[j] = pV[off_i1+j];
151 pV[off_i+j] = 0.0;
152 }
153 }
154 pD[i] = h;
155 }
156
157 // Accumulate transformations.
158
159 for (i = 0; i < n-1; i++) {
160 const Int_t off_i = i*n;
161 pV[off_n1+i] = pV[off_i+i];
162 pV[off_i+i] = 1.0;
163 Double_t h = pD[i+1];
164 if (h != 0.0) {
165 for (k = 0; k <= i; k++) {
166 const Int_t off_k = k*n;
167 pD[k] = pV[off_k+i+1]/h;
168 }
169 for (j = 0; j <= i; j++) {
170 Double_t g = 0.0;
171 for (k = 0; k <= i; k++) {
172 const Int_t off_k = k*n;
173 g += pV[off_k+i+1]*pV[off_k+j];
174 }
175 for (k = 0; k <= i; k++) {
176 const Int_t off_k = k*n;
177 pV[off_k+j] -= g*pD[k];
178 }
179 }
180 }
181 for (k = 0; k <= i; k++) {
182 const Int_t off_k = k*n;
183 pV[off_k+i+1] = 0.0;
184 }
185 }
186 for (j = 0; j < n; j++) {
187 pD[j] = pV[off_n1+j];
188 pV[off_n1+j] = 0.0;
189 }
190 pV[off_n1+n-1] = 1.0;
191 pE[0] = 0.0;
192}
193
194////////////////////////////////////////////////////////////////////////////////
195/// Symmetric tridiagonal QL algorithm.
196/// This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and
197/// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
198/// Fortran subroutine in EISPACK.
199
201{
202 Double_t *pV = v.GetMatrixArray();
203 Double_t *pD = d.GetMatrixArray();
204 Double_t *pE = e.GetMatrixArray();
205
206 const Int_t n = v.GetNrows();
207
208 Int_t i,j,k,l;
209 for (i = 1; i < n; i++)
210 pE[i-1] = pE[i];
211 pE[n-1] = 0.0;
212
213 Double_t f = 0.0;
214 Double_t tst1 = 0.0;
215 Double_t eps = TMath::Power(2.0,-52.0);
216 for (l = 0; l < n; l++) {
217
218 // Find small subdiagonal element
219
221 Int_t m = l;
222
223 // Original while-loop from Java code
224 while (m < n) {
225 if (TMath::Abs(pE[m]) <= eps*tst1) {
226 break;
227 }
228 m++;
229 }
230
231 // If m == l, pD[l] is an eigenvalue,
232 // otherwise, iterate.
233
234 if (m > l) {
235 Int_t iter = 0;
236 do {
237 if (iter++ == 30) { // (check iteration count here.)
238 Error("MakeEigenVectors","too many iterations");
239 break;
240 }
241
242 // Compute implicit shift
243
244 Double_t g = pD[l];
245 Double_t p = (pD[l+1]-g)/(2.0*pE[l]);
246 Double_t r = TMath::Hypot(p,1.0);
247 if (p < 0)
248 r = -r;
249 pD[l] = pE[l]/(p+r);
250 pD[l+1] = pE[l]*(p+r);
251 Double_t dl1 = pD[l+1];
252 Double_t h = g-pD[l];
253 for (i = l+2; i < n; i++)
254 pD[i] -= h;
255 f = f+h;
256
257 // Implicit QL transformation.
258
259 p = pD[m];
260 Double_t c = 1.0;
261 Double_t c2 = c;
262 Double_t c3 = c;
263 Double_t el1 = pE[l+1];
264 Double_t s = 0.0;
265 Double_t s2 = 0.0;
266 for (i = m-1; i >= l; i--) {
267 c3 = c2;
268 c2 = c;
269 s2 = s;
270 g = c*pE[i];
271 h = c*p;
272 r = TMath::Hypot(p,pE[i]);
273 pE[i+1] = s*r;
274 s = pE[i]/r;
275 c = p/r;
276 p = c*pD[i]-s*g;
277 pD[i+1] = h+s*(c*g+s*pD[i]);
278
279 // Accumulate transformation.
280
281 for (k = 0; k < n; k++) {
282 const Int_t off_k = k*n;
283 h = pV[off_k+i+1];
284 pV[off_k+i+1] = s*pV[off_k+i]+c*h;
285 pV[off_k+i] = c*pV[off_k+i]-s*h;
286 }
287 }
288 p = -s*s2*c3*el1*pE[l]/dl1;
289 pE[l] = s*p;
290 pD[l] = c*p;
291
292 // Check for convergence.
293
294 } while (TMath::Abs(pE[l]) > eps*tst1);
295 }
296 pD[l] = pD[l]+f;
297 pE[l] = 0.0;
298 }
299
300 // Sort eigenvalues and corresponding vectors.
301
302 for (i = 0; i < n-1; i++) {
303 k = i;
304 Double_t p = pD[i];
305 for (j = i+1; j < n; j++) {
306 if (pD[j] > p) {
307 k = j;
308 p = pD[j];
309 }
310 }
311 if (k != i) {
312 pD[k] = pD[i];
313 pD[i] = p;
314 for (j = 0; j < n; j++) {
315 const Int_t off_j = j*n;
316 p = pV[off_j+i];
317 pV[off_j+i] = pV[off_j+k];
318 pV[off_j+k] = p;
319 }
320 }
321 }
322}
323
324////////////////////////////////////////////////////////////////////////////////
325/// Assignment operator
326
328{
329 if (this != &source) {
330 fEigenVectors.ResizeTo(source.fEigenVectors);
331 fEigenValues.ResizeTo(source.fEigenValues);
332 fEigenVectors = source.fEigenVectors;
333 fEigenValues = source.fEigenValues;
334 }
335 return *this;
336}
#define d(i)
Definition RSha256.hxx:102
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
TMatrixDSymEigen.
static void MakeEigenVectors(TMatrixD &v, TVectorD &d, TVectorD &e)
Symmetric tridiagonal QL algorithm.
static void MakeTridiagonal(TMatrixD &v, TVectorD &d, TVectorD &e)
This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and Wilkinson,...
TMatrixDSymEigen & operator=(const TMatrixDSymEigen &source)
Assignment operator.
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:293
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14
return c3
Definition legend3.C:15
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4