```// @(#)root/matrix:\$Id\$
// Authors: Fons Rademakers, Eddy Offermann  Dec 2003

/*************************************************************************
*                                                                       *
* For the licensing terms see \$ROOTSYS/LICENSE.                         *
* For the list of contributors see \$ROOTSYS/README/CREDITS.             *
*************************************************************************/

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TMatrixDSymEigen                                                     //
//                                                                      //
// Eigenvalues and eigenvectors of a real symmetric matrix.             //
//                                                                      //
// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is  //
// diagonal and the eigenvector matrix V is orthogonal. That is, the    //
// diagonal values of D are the eigenvalues, and V*V' = I, where I is   //
// the identity matrix.  The columns of V represent the eigenvectors in //
// the sense that A*V = V*D.                                            //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "TMatrixDSymEigen.h"
#include "TMath.h"

ClassImp(TMatrixDSymEigen)

//______________________________________________________________________________
TMatrixDSymEigen::TMatrixDSymEigen(const TMatrixDSym &a)
{
// Constructor for eigen-problem of symmetric matrix A .

R__ASSERT(a.IsValid());

const Int_t nRows  = a.GetNrows();
const Int_t rowLwb = a.GetRowLwb();

fEigenValues.ResizeTo(rowLwb,rowLwb+nRows-1);
fEigenVectors.ResizeTo(a);

fEigenVectors = a;

TVectorD offDiag;
Double_t work[kWorkMax];
if (nRows > kWorkMax) offDiag.ResizeTo(nRows);
else                  offDiag.Use(nRows,work);

// Tridiagonalize matrix
MakeTridiagonal(fEigenVectors,fEigenValues,offDiag);

// Make eigenvectors and -values
MakeEigenVectors(fEigenVectors,fEigenValues,offDiag);
}

//______________________________________________________________________________
TMatrixDSymEigen::TMatrixDSymEigen(const TMatrixDSymEigen &another)
{
// Copy constructor

*this = another;
}

//______________________________________________________________________________
void TMatrixDSymEigen::MakeTridiagonal(TMatrixD &v,TVectorD &d,TVectorD &e)
{
// This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and
// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.

Double_t *pV = v.GetMatrixArray();
Double_t *pD = d.GetMatrixArray();
Double_t *pE = e.GetMatrixArray();

const Int_t n = v.GetNrows();

Int_t i,j,k;
Int_t off_n1 = (n-1)*n;
for (j = 0; j < n; j++)
pD[j] = pV[off_n1+j];

// Householder reduction to tridiagonal form.

for (i = n-1; i > 0; i--) {
const Int_t off_i1 = (i-1)*n;
const Int_t off_i  = i*n;

// Scale to avoid under/overflow.

Double_t scale = 0.0;
Double_t h = 0.0;
for (k = 0; k < i; k++)
scale = scale+TMath::Abs(pD[k]);
if (scale == 0.0) {
pE[i] = pD[i-1];
for (j = 0; j < i; j++) {
const Int_t off_j = j*n;
pD[j] = pV[off_i1+j];
pV[off_i+j] = 0.0;
pV[off_j+i] = 0.0;
}
} else {

// Generate Householder vector.

for (k = 0; k < i; k++) {
pD[k] /= scale;
h += pD[k]*pD[k];
}
Double_t f = pD[i-1];
Double_t g = TMath::Sqrt(h);
if (f > 0)
g = -g;
pE[i]   = scale*g;
h       = h-f*g;
pD[i-1] = f-g;
for (j = 0; j < i; j++)
pE[j] = 0.0;

// Apply similarity transformation to remaining columns.

for (j = 0; j < i; j++) {
const Int_t off_j = j*n;
f = pD[j];
pV[off_j+i] = f;
g = pE[j]+pV[off_j+j]*f;
for (k = j+1; k <= i-1; k++) {
const Int_t off_k = k*n;
g += pV[off_k+j]*pD[k];
pE[k] += pV[off_k+j]*f;
}
pE[j] = g;
}
f = 0.0;
for (j = 0; j < i; j++) {
pE[j] /= h;
f += pE[j]*pD[j];
}
Double_t hh = f/(h+h);
for (j = 0; j < i; j++)
pE[j] -= hh*pD[j];
for (j = 0; j < i; j++) {
f = pD[j];
g = pE[j];
for (k = j; k <= i-1; k++) {
const Int_t off_k = k*n;
pV[off_k+j] -= (f*pE[k]+g*pD[k]);
}
pD[j] = pV[off_i1+j];
pV[off_i+j] = 0.0;
}
}
pD[i] = h;
}

// Accumulate transformations.

for (i = 0; i < n-1; i++) {
const Int_t off_i  = i*n;
pV[off_n1+i] = pV[off_i+i];
pV[off_i+i] = 1.0;
Double_t h = pD[i+1];
if (h != 0.0) {
for (k = 0; k <= i; k++) {
const Int_t off_k = k*n;
pD[k] = pV[off_k+i+1]/h;
}
for (j = 0; j <= i; j++) {
Double_t g = 0.0;
for (k = 0; k <= i; k++) {
const Int_t off_k = k*n;
g += pV[off_k+i+1]*pV[off_k+j];
}
for (k = 0; k <= i; k++) {
const Int_t off_k = k*n;
pV[off_k+j] -= g*pD[k];
}
}
}
for (k = 0; k <= i; k++) {
const Int_t off_k = k*n;
pV[off_k+i+1] = 0.0;
}
}
for (j = 0; j < n; j++) {
pD[j] = pV[off_n1+j];
pV[off_n1+j] = 0.0;
}
pV[off_n1+n-1] = 1.0;
pE[0] = 0.0;
}

//______________________________________________________________________________
void TMatrixDSymEigen::MakeEigenVectors(TMatrixD &v,TVectorD &d,TVectorD &e)
{
// Symmetric tridiagonal QL algorithm.
// This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and
// Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.

Double_t *pV = v.GetMatrixArray();
Double_t *pD = d.GetMatrixArray();
Double_t *pE = e.GetMatrixArray();

const Int_t n = v.GetNrows();

Int_t i,j,k,l;
for (i = 1; i < n; i++)
pE[i-1] = pE[i];
pE[n-1] = 0.0;

Double_t f = 0.0;
Double_t tst1 = 0.0;
Double_t eps = TMath::Power(2.0,-52.0);
for (l = 0; l < n; l++) {

// Find small subdiagonal element

tst1 = TMath::Max(tst1,TMath::Abs(pD[l])+TMath::Abs(pE[l]));
Int_t m = l;

// Original while-loop from Java code
while (m < n) {
if (TMath::Abs(pE[m]) <= eps*tst1) {
break;
}
m++;
}

// If m == l, pD[l] is an eigenvalue,
// otherwise, iterate.

if (m > l) {
Int_t iter = 0;
do {
if (iter++ == 30) {  // (check iteration count here.)
Error("MakeEigenVectors","too many iterations");
break;
}

// Compute implicit shift

Double_t g = pD[l];
Double_t p = (pD[l+1]-g)/(2.0*pE[l]);
Double_t r = TMath::Hypot(p,1.0);
if (p < 0)
r = -r;
pD[l] = pE[l]/(p+r);
pD[l+1] = pE[l]*(p+r);
Double_t dl1 = pD[l+1];
Double_t h = g-pD[l];
for (i = l+2; i < n; i++)
pD[i] -= h;
f = f+h;

// Implicit QL transformation.

p = pD[m];
Double_t c = 1.0;
Double_t c2 = c;
Double_t c3 = c;
Double_t el1 = pE[l+1];
Double_t s = 0.0;
Double_t s2 = 0.0;
for (i = m-1; i >= l; i--) {
c3 = c2;
c2 = c;
s2 = s;
g = c*pE[i];
h = c*p;
r = TMath::Hypot(p,pE[i]);
pE[i+1] = s*r;
s = pE[i]/r;
c = p/r;
p = c*pD[i]-s*g;
pD[i+1] = h+s*(c*g+s*pD[i]);

// Accumulate transformation.

for (k = 0; k < n; k++) {
const Int_t off_k = k*n;
h = pV[off_k+i+1];
pV[off_k+i+1] = s*pV[off_k+i]+c*h;
pV[off_k+i]   = c*pV[off_k+i]-s*h;
}
}
p = -s*s2*c3*el1*pE[l]/dl1;
pE[l] = s*p;
pD[l] = c*p;

// Check for convergence.

} while (TMath::Abs(pE[l]) > eps*tst1);
}
pD[l] = pD[l]+f;
pE[l] = 0.0;
}

// Sort eigenvalues and corresponding vectors.

for (i = 0; i < n-1; i++) {
k = i;
Double_t p = pD[i];
for (j = i+1; j < n; j++) {
if (pD[j] > p) {
k = j;
p = pD[j];
}
}
if (k != i) {
pD[k] = pD[i];
pD[i] = p;
for (j = 0; j < n; j++) {
const Int_t off_j = j*n;
p = pV[off_j+i];
pV[off_j+i] = pV[off_j+k];
pV[off_j+k] = p;
}
}
}
}

//______________________________________________________________________________
TMatrixDSymEigen &TMatrixDSymEigen::operator=(const TMatrixDSymEigen &source)
{
// Assignment operator

if (this != &source) {
fEigenVectors.ResizeTo(source.fEigenVectors);
fEigenValues.ResizeTo(source.fEigenValues);
}
return *this;
}
```
TMatrixDSymEigen.cxx:1
TMatrixDSymEigen.cxx:2
TMatrixDSymEigen.cxx:3
TMatrixDSymEigen.cxx:4
TMatrixDSymEigen.cxx:5
TMatrixDSymEigen.cxx:6
TMatrixDSymEigen.cxx:7
TMatrixDSymEigen.cxx:8
TMatrixDSymEigen.cxx:9
TMatrixDSymEigen.cxx:10
TMatrixDSymEigen.cxx:11
TMatrixDSymEigen.cxx:12
TMatrixDSymEigen.cxx:13
TMatrixDSymEigen.cxx:14
TMatrixDSymEigen.cxx:15
TMatrixDSymEigen.cxx:16
TMatrixDSymEigen.cxx:17
TMatrixDSymEigen.cxx:18
TMatrixDSymEigen.cxx:19
TMatrixDSymEigen.cxx:20
TMatrixDSymEigen.cxx:21
TMatrixDSymEigen.cxx:22
TMatrixDSymEigen.cxx:23
TMatrixDSymEigen.cxx:24
TMatrixDSymEigen.cxx:25
TMatrixDSymEigen.cxx:26
TMatrixDSymEigen.cxx:27
TMatrixDSymEigen.cxx:28
TMatrixDSymEigen.cxx:29
TMatrixDSymEigen.cxx:30
TMatrixDSymEigen.cxx:31
TMatrixDSymEigen.cxx:32
TMatrixDSymEigen.cxx:33
TMatrixDSymEigen.cxx:34
TMatrixDSymEigen.cxx:35
TMatrixDSymEigen.cxx:36
TMatrixDSymEigen.cxx:37
TMatrixDSymEigen.cxx:38
TMatrixDSymEigen.cxx:39
TMatrixDSymEigen.cxx:40
TMatrixDSymEigen.cxx:41
TMatrixDSymEigen.cxx:42
TMatrixDSymEigen.cxx:43
TMatrixDSymEigen.cxx:44
TMatrixDSymEigen.cxx:45
TMatrixDSymEigen.cxx:46
TMatrixDSymEigen.cxx:47
TMatrixDSymEigen.cxx:48
TMatrixDSymEigen.cxx:49
TMatrixDSymEigen.cxx:50
TMatrixDSymEigen.cxx:51
TMatrixDSymEigen.cxx:52
TMatrixDSymEigen.cxx:53
TMatrixDSymEigen.cxx:54
TMatrixDSymEigen.cxx:55
TMatrixDSymEigen.cxx:56
TMatrixDSymEigen.cxx:57
TMatrixDSymEigen.cxx:58
TMatrixDSymEigen.cxx:59
TMatrixDSymEigen.cxx:60
TMatrixDSymEigen.cxx:61
TMatrixDSymEigen.cxx:62
TMatrixDSymEigen.cxx:63
TMatrixDSymEigen.cxx:64
TMatrixDSymEigen.cxx:65
TMatrixDSymEigen.cxx:66
TMatrixDSymEigen.cxx:67
TMatrixDSymEigen.cxx:68
TMatrixDSymEigen.cxx:69
TMatrixDSymEigen.cxx:70
TMatrixDSymEigen.cxx:71
TMatrixDSymEigen.cxx:72
TMatrixDSymEigen.cxx:73
TMatrixDSymEigen.cxx:74
TMatrixDSymEigen.cxx:75
TMatrixDSymEigen.cxx:76
TMatrixDSymEigen.cxx:77
TMatrixDSymEigen.cxx:78
TMatrixDSymEigen.cxx:79
TMatrixDSymEigen.cxx:80
TMatrixDSymEigen.cxx:81
TMatrixDSymEigen.cxx:82
TMatrixDSymEigen.cxx:83
TMatrixDSymEigen.cxx:84
TMatrixDSymEigen.cxx:85
TMatrixDSymEigen.cxx:86
TMatrixDSymEigen.cxx:87
TMatrixDSymEigen.cxx:88
TMatrixDSymEigen.cxx:89
TMatrixDSymEigen.cxx:90
TMatrixDSymEigen.cxx:91
TMatrixDSymEigen.cxx:92
TMatrixDSymEigen.cxx:93
TMatrixDSymEigen.cxx:94
TMatrixDSymEigen.cxx:95
TMatrixDSymEigen.cxx:96
TMatrixDSymEigen.cxx:97
TMatrixDSymEigen.cxx:98
TMatrixDSymEigen.cxx:99
TMatrixDSymEigen.cxx:100
TMatrixDSymEigen.cxx:101
TMatrixDSymEigen.cxx:102
TMatrixDSymEigen.cxx:103
TMatrixDSymEigen.cxx:104
TMatrixDSymEigen.cxx:105
TMatrixDSymEigen.cxx:106
TMatrixDSymEigen.cxx:107
TMatrixDSymEigen.cxx:108
TMatrixDSymEigen.cxx:109
TMatrixDSymEigen.cxx:110
TMatrixDSymEigen.cxx:111
TMatrixDSymEigen.cxx:112
TMatrixDSymEigen.cxx:113
TMatrixDSymEigen.cxx:114
TMatrixDSymEigen.cxx:115
TMatrixDSymEigen.cxx:116
TMatrixDSymEigen.cxx:117
TMatrixDSymEigen.cxx:118
TMatrixDSymEigen.cxx:119
TMatrixDSymEigen.cxx:120
TMatrixDSymEigen.cxx:121
TMatrixDSymEigen.cxx:122
TMatrixDSymEigen.cxx:123
TMatrixDSymEigen.cxx:124
TMatrixDSymEigen.cxx:125
TMatrixDSymEigen.cxx:126
TMatrixDSymEigen.cxx:127
TMatrixDSymEigen.cxx:128
TMatrixDSymEigen.cxx:129
TMatrixDSymEigen.cxx:130
TMatrixDSymEigen.cxx:131
TMatrixDSymEigen.cxx:132
TMatrixDSymEigen.cxx:133
TMatrixDSymEigen.cxx:134
TMatrixDSymEigen.cxx:135
TMatrixDSymEigen.cxx:136
TMatrixDSymEigen.cxx:137
TMatrixDSymEigen.cxx:138
TMatrixDSymEigen.cxx:139
TMatrixDSymEigen.cxx:140
TMatrixDSymEigen.cxx:141
TMatrixDSymEigen.cxx:142
TMatrixDSymEigen.cxx:143
TMatrixDSymEigen.cxx:144
TMatrixDSymEigen.cxx:145
TMatrixDSymEigen.cxx:146
TMatrixDSymEigen.cxx:147
TMatrixDSymEigen.cxx:148
TMatrixDSymEigen.cxx:149
TMatrixDSymEigen.cxx:150
TMatrixDSymEigen.cxx:151
TMatrixDSymEigen.cxx:152
TMatrixDSymEigen.cxx:153
TMatrixDSymEigen.cxx:154
TMatrixDSymEigen.cxx:155
TMatrixDSymEigen.cxx:156
TMatrixDSymEigen.cxx:157
TMatrixDSymEigen.cxx:158
TMatrixDSymEigen.cxx:159
TMatrixDSymEigen.cxx:160
TMatrixDSymEigen.cxx:161
TMatrixDSymEigen.cxx:162
TMatrixDSymEigen.cxx:163
TMatrixDSymEigen.cxx:164
TMatrixDSymEigen.cxx:165
TMatrixDSymEigen.cxx:166
TMatrixDSymEigen.cxx:167
TMatrixDSymEigen.cxx:168
TMatrixDSymEigen.cxx:169
TMatrixDSymEigen.cxx:170
TMatrixDSymEigen.cxx:171
TMatrixDSymEigen.cxx:172
TMatrixDSymEigen.cxx:173
TMatrixDSymEigen.cxx:174
TMatrixDSymEigen.cxx:175
TMatrixDSymEigen.cxx:176
TMatrixDSymEigen.cxx:177
TMatrixDSymEigen.cxx:178
TMatrixDSymEigen.cxx:179
TMatrixDSymEigen.cxx:180
TMatrixDSymEigen.cxx:181
TMatrixDSymEigen.cxx:182
TMatrixDSymEigen.cxx:183
TMatrixDSymEigen.cxx:184
TMatrixDSymEigen.cxx:185
TMatrixDSymEigen.cxx:186
TMatrixDSymEigen.cxx:187
TMatrixDSymEigen.cxx:188
TMatrixDSymEigen.cxx:189
TMatrixDSymEigen.cxx:190
TMatrixDSymEigen.cxx:191
TMatrixDSymEigen.cxx:192
TMatrixDSymEigen.cxx:193
TMatrixDSymEigen.cxx:194
TMatrixDSymEigen.cxx:195
TMatrixDSymEigen.cxx:196
TMatrixDSymEigen.cxx:197
TMatrixDSymEigen.cxx:198
TMatrixDSymEigen.cxx:199
TMatrixDSymEigen.cxx:200
TMatrixDSymEigen.cxx:201
TMatrixDSymEigen.cxx:202
TMatrixDSymEigen.cxx:203
TMatrixDSymEigen.cxx:204
TMatrixDSymEigen.cxx:205
TMatrixDSymEigen.cxx:206
TMatrixDSymEigen.cxx:207
TMatrixDSymEigen.cxx:208
TMatrixDSymEigen.cxx:209
TMatrixDSymEigen.cxx:210
TMatrixDSymEigen.cxx:211
TMatrixDSymEigen.cxx:212
TMatrixDSymEigen.cxx:213
TMatrixDSymEigen.cxx:214
TMatrixDSymEigen.cxx:215
TMatrixDSymEigen.cxx:216
TMatrixDSymEigen.cxx:217
TMatrixDSymEigen.cxx:218
TMatrixDSymEigen.cxx:219
TMatrixDSymEigen.cxx:220
TMatrixDSymEigen.cxx:221
TMatrixDSymEigen.cxx:222
TMatrixDSymEigen.cxx:223
TMatrixDSymEigen.cxx:224
TMatrixDSymEigen.cxx:225
TMatrixDSymEigen.cxx:226
TMatrixDSymEigen.cxx:227
TMatrixDSymEigen.cxx:228
TMatrixDSymEigen.cxx:229
TMatrixDSymEigen.cxx:230
TMatrixDSymEigen.cxx:231
TMatrixDSymEigen.cxx:232
TMatrixDSymEigen.cxx:233
TMatrixDSymEigen.cxx:234
TMatrixDSymEigen.cxx:235
TMatrixDSymEigen.cxx:236
TMatrixDSymEigen.cxx:237
TMatrixDSymEigen.cxx:238
TMatrixDSymEigen.cxx:239
TMatrixDSymEigen.cxx:240
TMatrixDSymEigen.cxx:241
TMatrixDSymEigen.cxx:242
TMatrixDSymEigen.cxx:243
TMatrixDSymEigen.cxx:244
TMatrixDSymEigen.cxx:245
TMatrixDSymEigen.cxx:246
TMatrixDSymEigen.cxx:247
TMatrixDSymEigen.cxx:248
TMatrixDSymEigen.cxx:249
TMatrixDSymEigen.cxx:250
TMatrixDSymEigen.cxx:251
TMatrixDSymEigen.cxx:252
TMatrixDSymEigen.cxx:253
TMatrixDSymEigen.cxx:254
TMatrixDSymEigen.cxx:255
TMatrixDSymEigen.cxx:256
TMatrixDSymEigen.cxx:257
TMatrixDSymEigen.cxx:258
TMatrixDSymEigen.cxx:259
TMatrixDSymEigen.cxx:260
TMatrixDSymEigen.cxx:261
TMatrixDSymEigen.cxx:262
TMatrixDSymEigen.cxx:263
TMatrixDSymEigen.cxx:264
TMatrixDSymEigen.cxx:265
TMatrixDSymEigen.cxx:266
TMatrixDSymEigen.cxx:267
TMatrixDSymEigen.cxx:268
TMatrixDSymEigen.cxx:269
TMatrixDSymEigen.cxx:270
TMatrixDSymEigen.cxx:271
TMatrixDSymEigen.cxx:272
TMatrixDSymEigen.cxx:273
TMatrixDSymEigen.cxx:274
TMatrixDSymEigen.cxx:275
TMatrixDSymEigen.cxx:276
TMatrixDSymEigen.cxx:277
TMatrixDSymEigen.cxx:278
TMatrixDSymEigen.cxx:279
TMatrixDSymEigen.cxx:280
TMatrixDSymEigen.cxx:281
TMatrixDSymEigen.cxx:282
TMatrixDSymEigen.cxx:283
TMatrixDSymEigen.cxx:284
TMatrixDSymEigen.cxx:285
TMatrixDSymEigen.cxx:286
TMatrixDSymEigen.cxx:287
TMatrixDSymEigen.cxx:288
TMatrixDSymEigen.cxx:289
TMatrixDSymEigen.cxx:290
TMatrixDSymEigen.cxx:291
TMatrixDSymEigen.cxx:292
TMatrixDSymEigen.cxx:293
TMatrixDSymEigen.cxx:294
TMatrixDSymEigen.cxx:295
TMatrixDSymEigen.cxx:296
TMatrixDSymEigen.cxx:297
TMatrixDSymEigen.cxx:298
TMatrixDSymEigen.cxx:299
TMatrixDSymEigen.cxx:300
TMatrixDSymEigen.cxx:301
TMatrixDSymEigen.cxx:302
TMatrixDSymEigen.cxx:303
TMatrixDSymEigen.cxx:304
TMatrixDSymEigen.cxx:305
TMatrixDSymEigen.cxx:306
TMatrixDSymEigen.cxx:307
TMatrixDSymEigen.cxx:308
TMatrixDSymEigen.cxx:309
TMatrixDSymEigen.cxx:310
TMatrixDSymEigen.cxx:311
TMatrixDSymEigen.cxx:312
TMatrixDSymEigen.cxx:313
TMatrixDSymEigen.cxx:314
TMatrixDSymEigen.cxx:315
TMatrixDSymEigen.cxx:316
TMatrixDSymEigen.cxx:317
TMatrixDSymEigen.cxx:318
TMatrixDSymEigen.cxx:319
TMatrixDSymEigen.cxx:320
TMatrixDSymEigen.cxx:321
TMatrixDSymEigen.cxx:322
TMatrixDSymEigen.cxx:323
TMatrixDSymEigen.cxx:324
TMatrixDSymEigen.cxx:325
TMatrixDSymEigen.cxx:326
TMatrixDSymEigen.cxx:327
TMatrixDSymEigen.cxx:328
TMatrixDSymEigen.cxx:329
TMatrixDSymEigen.cxx:330
TMatrixDSymEigen.cxx:331
TMatrixDSymEigen.cxx:332
TMatrixDSymEigen.cxx:333
TMatrixDSymEigen.cxx:334
TMatrixDSymEigen.cxx:335