Logo ROOT  
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/** \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
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);
43
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
53
54 // Make eigenvectors and -values
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];
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) {
334 fEigenValues = source.fEigenValues;
335 }
336 return *this;
337}
ROOT::R::TRInterface & r
Definition: Object.C:4
#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 h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
double Double_t
Definition: RtypesCore.h:57
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
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.
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:1213
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:295
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:350
const Int_t n
Definition: legend1.C:16
return c2
Definition: legend2.C:14
return c3
Definition: legend3.C:15
static constexpr double s
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:57
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12