Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TMatrixDEigen.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 TMatrixDEigen
13  \ingroup Matrix
14 
15  TMatrixDEigen
16 
17  Eigenvalues and eigenvectors of a real matrix.
18 
19  If A is not symmetric, then the eigenvalue matrix D is block
20  diagonal with the real eigenvalues in 1-by-1 blocks and any complex
21  eigenvalues, a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if
22  the complex eigenvalues look like
23 
24 ~~~
25  u + iv . . . . .
26  . u - iv . . . .
27  . . a + ib . . .
28  . . . a - ib . .
29  . . . . x .
30  . . . . . y
31 ~~~
32 
33  then D looks like
34 
35 ~~~
36  u v . . . .
37  -v u . . . .
38  . . a b . .
39  . . -b a . .
40  . . . . x .
41  . . . . . y
42 ~~~
43 
44  This keeps V a real matrix in both symmetric and non-symmetric
45  cases, and A*V = V*D.
46 
47 */
48 
49 #include "TMatrixDEigen.h"
50 #include "TMath.h"
51 
52 ClassImp(TMatrixDEigen);
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 /// Constructor for eigen-problem of matrix A .
56 
57 TMatrixDEigen::TMatrixDEigen(const TMatrixD &a)
58 {
59  R__ASSERT(a.IsValid());
60 
61  const Int_t nRows = a.GetNrows();
62  const Int_t nCols = a.GetNcols();
63  const Int_t rowLwb = a.GetRowLwb();
64  const Int_t colLwb = a.GetColLwb();
65 
66  if (nRows != nCols || rowLwb != colLwb)
67  {
68  Error("TMatrixDEigen(TMatrixD &)","matrix should be square");
69  return;
70  }
71 
72  const Int_t rowUpb = rowLwb+nRows-1;
73  fEigenVectors.ResizeTo(rowLwb,rowUpb,rowLwb,rowUpb);
74  fEigenValuesRe.ResizeTo(rowLwb,rowUpb);
75  fEigenValuesIm.ResizeTo(rowLwb,rowUpb);
76 
77  TVectorD ortho;
78  Double_t work[kWorkMax];
79  if (nRows > kWorkMax) ortho.ResizeTo(nRows);
80  else ortho.Use(nRows,work);
81 
82  TMatrixD mH = a;
83 
84  // Reduce to Hessenberg form.
85  MakeHessenBerg(fEigenVectors,ortho,mH);
86 
87  // Reduce Hessenberg to real Schur form.
88  MakeSchurr(fEigenVectors,fEigenValuesRe,fEigenValuesIm,mH);
89 
90  // Sort eigenvalues and corresponding vectors in descending order of Re^2+Im^2
91  // of the complex eigenvalues .
92  Sort(fEigenVectors,fEigenValuesRe,fEigenValuesIm);
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Copy constructor
97 
98 TMatrixDEigen::TMatrixDEigen(const TMatrixDEigen &another)
99 {
100  *this = another;
101 }
102 
103 ////////////////////////////////////////////////////////////////////////////////
104 /// Nonsymmetric reduction to Hessenberg form.
105 /// This is derived from the Algol procedures orthes and ortran, by Martin and Wilkinson,
106 /// Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
107 /// Fortran subroutines in EISPACK.
108 
109 void TMatrixDEigen::MakeHessenBerg(TMatrixD &v,TVectorD &ortho,TMatrixD &H)
110 {
111  Double_t *pV = v.GetMatrixArray();
112  Double_t *pO = ortho.GetMatrixArray();
113  Double_t *pH = H.GetMatrixArray();
114 
115  const UInt_t n = v.GetNrows();
116 
117  const UInt_t low = 0;
118  const UInt_t high = n-1;
119 
120  UInt_t i,j,m;
121  for (m = low+1; m <= high-1; m++) {
122  const UInt_t off_m = m*n;
123 
124  // Scale column.
125 
126  Double_t scale = 0.0;
127  for (i = m; i <= high; i++) {
128  const UInt_t off_i = i*n;
129  scale = scale + TMath::Abs(pH[off_i+m-1]);
130  }
131  if (scale != 0.0) {
132 
133  // Compute Householder transformation.
134 
135  Double_t h = 0.0;
136  for (i = high; i >= m; i--) {
137  const UInt_t off_i = i*n;
138  pO[i] = pH[off_i+m-1]/scale;
139  h += pO[i]*pO[i];
140  }
141  Double_t g = TMath::Sqrt(h);
142  if (pO[m] > 0)
143  g = -g;
144  h = h-pO[m]*g;
145  pO[m] = pO[m]-g;
146 
147  // Apply Householder similarity transformation
148  // H = (I-u*u'/h)*H*(I-u*u')/h)
149 
150  for (j = m; j < n; j++) {
151  Double_t f = 0.0;
152  for (i = high; i >= m; i--) {
153  const UInt_t off_i = i*n;
154  f += pO[i]*pH[off_i+j];
155  }
156  f = f/h;
157  for (i = m; i <= high; i++) {
158  const UInt_t off_i = i*n;
159  pH[off_i+j] -= f*pO[i];
160  }
161  }
162 
163  for (i = 0; i <= high; i++) {
164  const UInt_t off_i = i*n;
165  Double_t f = 0.0;
166  for (j = high; j >= m; j--)
167  f += pO[j]*pH[off_i+j];
168  f = f/h;
169  for (j = m; j <= high; j++)
170  pH[off_i+j] -= f*pO[j];
171  }
172  pO[m] = scale*pO[m];
173  pH[off_m+m-1] = scale*g;
174  }
175  }
176 
177  // Accumulate transformations (Algol's ortran).
178 
179  for (i = 0; i < n; i++) {
180  const UInt_t off_i = i*n;
181  for (j = 0; j < n; j++)
182  pV[off_i+j] = (i == j ? 1.0 : 0.0);
183  }
184 
185  for (m = high-1; m >= low+1; m--) {
186  const UInt_t off_m = m*n;
187  if (pH[off_m+m-1] != 0.0) {
188  for (i = m+1; i <= high; i++) {
189  const UInt_t off_i = i*n;
190  pO[i] = pH[off_i+m-1];
191  }
192  for (j = m; j <= high; j++) {
193  Double_t g = 0.0;
194  for (i = m; i <= high; i++) {
195  const UInt_t off_i = i*n;
196  g += pO[i]*pV[off_i+j];
197  }
198  // Double division avoids possible underflow
199  g = (g/pO[m])/pH[off_m+m-1];
200  for (i = m; i <= high; i++) {
201  const UInt_t off_i = i*n;
202  pV[off_i+j] += g*pO[i];
203  }
204  }
205  }
206  }
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////////
210 /// Complex scalar division.
211 
212 static Double_t gCdivr, gCdivi;
213 static void cdiv(Double_t xr,Double_t xi,Double_t yr,Double_t yi) {
214  Double_t r,d;
215  if (TMath::Abs(yr) > TMath::Abs(yi)) {
216  r = yi/yr;
217  d = yr+r*yi;
218  gCdivr = (xr+r*xi)/d;
219  gCdivi = (xi-r*xr)/d;
220  } else {
221  r = yr/yi;
222  d = yi+r*yr;
223  gCdivr = (r*xr+xi)/d;
224  gCdivi = (r*xi-xr)/d;
225  }
226 }
227 
228 ////////////////////////////////////////////////////////////////////////////////
229 /// Nonsymmetric reduction from Hessenberg to real Schur form.
230 /// This is derived from the Algol procedure hqr2, by Martin and Wilkinson,
231 /// Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
232 /// Fortran subroutine in EISPACK.
233 
234 void TMatrixDEigen::MakeSchurr(TMatrixD &v,TVectorD &d,TVectorD &e,TMatrixD &H)
235 {
236  // Initialize
237 
238  const Int_t nn = v.GetNrows();
239  Int_t n = nn-1;
240  const Int_t low = 0;
241  const Int_t high = nn-1;
242  const Double_t eps = TMath::Power(2.0,-52.0);
243  Double_t exshift = 0.0;
244  Double_t p=0,q=0,r=0,s=0,z=0,t,w,x,y;
245 
246  Double_t *pV = v.GetMatrixArray();
247  Double_t *pD = d.GetMatrixArray();
248  Double_t *pE = e.GetMatrixArray();
249  Double_t *pH = H.GetMatrixArray();
250 
251  // Store roots isolated by balanc and compute matrix norm
252 
253  Double_t norm = 0.0;
254  Int_t i,j,k;
255  for (i = 0; i < nn; i++) {
256  const Int_t off_i = i*nn;
257  if ((i < low) || (i > high)) {
258  pD[i] = pH[off_i+i];
259  pE[i] = 0.0;
260  }
261  for (j = TMath::Max(i-1,0); j < nn; j++)
262  norm += TMath::Abs(pH[off_i+j]);
263  }
264 
265  // Outer loop over eigenvalue index
266 
267  Int_t iter = 0;
268  while (n >= low) {
269  const Int_t off_n = n*nn;
270  const Int_t off_n1 = (n-1)*nn;
271 
272  // Look for single small sub-diagonal element
273 
274  Int_t l = n;
275  while (l > low) {
276  const Int_t off_l1 = (l-1)*nn;
277  const Int_t off_l = l*nn;
278  s = TMath::Abs(pH[off_l1+l-1])+TMath::Abs(pH[off_l+l]);
279  if (s == 0.0)
280  s = norm;
281  if (TMath::Abs(pH[off_l+l-1]) < eps*s)
282  break;
283  l--;
284  }
285 
286  // Check for convergence
287  // One root found
288 
289  if (l == n) {
290  pH[off_n+n] = pH[off_n+n]+exshift;
291  pD[n] = pH[off_n+n];
292  pE[n] = 0.0;
293  n--;
294  iter = 0;
295 
296  // Two roots found
297 
298  } else if (l == n-1) {
299  w = pH[off_n+n-1]*pH[off_n1+n];
300  p = (pH[off_n1+n-1]-pH[off_n+n])/2.0;
301  q = p*p+w;
302  z = TMath::Sqrt(TMath::Abs(q));
303  pH[off_n+n] = pH[off_n+n]+exshift;
304  pH[off_n1+n-1] = pH[off_n1+n-1]+exshift;
305  x = pH[off_n+n];
306 
307  // Double_t pair
308 
309  if (q >= 0) {
310  if (p >= 0)
311  z = p+z;
312  else
313  z = p-z;
314  pD[n-1] = x+z;
315  pD[n] = pD[n-1];
316  if (z != 0.0)
317  pD[n] = x-w/z;
318  pE[n-1] = 0.0;
319  pE[n] = 0.0;
320  x = pH[off_n+n-1];
321  s = TMath::Abs(x)+TMath::Abs(z);
322  p = x/s;
323  q = z/s;
324  r = TMath::Sqrt((p*p)+(q*q));
325  p = p/r;
326  q = q/r;
327 
328  // Row modification
329 
330  for (j = n-1; j < nn; j++) {
331  z = pH[off_n1+j];
332  pH[off_n1+j] = q*z+p*pH[off_n+j];
333  pH[off_n+j] = q*pH[off_n+j]-p*z;
334  }
335 
336  // Column modification
337 
338  for (i = 0; i <= n; i++) {
339  const Int_t off_i = i*nn;
340  z = pH[off_i+n-1];
341  pH[off_i+n-1] = q*z+p*pH[off_i+n];
342  pH[off_i+n] = q*pH[off_i+n]-p*z;
343  }
344 
345  // Accumulate transformations
346 
347  for (i = low; i <= high; i++) {
348  const Int_t off_i = i*nn;
349  z = pV[off_i+n-1];
350  pV[off_i+n-1] = q*z+p*pV[off_i+n];
351  pV[off_i+n] = q*pV[off_i+n]-p*z;
352  }
353 
354  // Complex pair
355 
356  } else {
357  pD[n-1] = x+p;
358  pD[n] = x+p;
359  pE[n-1] = z;
360  pE[n] = -z;
361  }
362  n = n-2;
363  iter = 0;
364 
365  // No convergence yet
366 
367  } else {
368 
369  // Form shift
370 
371  x = pH[off_n+n];
372  y = 0.0;
373  w = 0.0;
374  if (l < n) {
375  y = pH[off_n1+n-1];
376  w = pH[off_n+n-1]*pH[off_n1+n];
377  }
378 
379  // Wilkinson's original ad hoc shift
380 
381  if (iter == 10) {
382  exshift += x;
383  for (i = low; i <= n; i++) {
384  const Int_t off_i = i*nn;
385  pH[off_i+i] -= x;
386  }
387  s = TMath::Abs(pH[off_n+n-1])+TMath::Abs(pH[off_n1+n-2]);
388  x = y = 0.75*s;
389  w = -0.4375*s*s;
390  }
391 
392  // MATLAB's new ad hoc shift
393 
394  if (iter == 30) {
395  s = (y-x)/2.0;
396  s = s*s+w;
397  if (s > 0) {
398  s = TMath::Sqrt(s);
399  if (y<x)
400  s = -s;
401  s = x-w/((y-x)/2.0+s);
402  for (i = low; i <= n; i++) {
403  const Int_t off_i = i*nn;
404  pH[off_i+i] -= s;
405  }
406  exshift += s;
407  x = y = w = 0.964;
408  }
409  }
410 
411  if (iter++ == 50) { // (check iteration count here.)
412  Error("MakeSchurr","too many iterations");
413  break;
414  }
415 
416  // Look for two consecutive small sub-diagonal elements
417 
418  Int_t m = n-2;
419  while (m >= l) {
420  const Int_t off_m = m*nn;
421  const Int_t off_m_1 = (m-1)*nn;
422  const Int_t off_m1 = (m+1)*nn;
423  const Int_t off_m2 = (m+2)*nn;
424  z = pH[off_m+m];
425  r = x-z;
426  s = y-z;
427  p = (r*s-w)/pH[off_m1+m]+pH[off_m+m+1];
428  q = pH[off_m1+m+1]-z-r-s;
429  r = pH[off_m2+m+1];
430  s = TMath::Abs(p)+TMath::Abs(q)+TMath::Abs(r);
431  p = p /s;
432  q = q/s;
433  r = r/s;
434  if (m == l)
435  break;
436  if (TMath::Abs(pH[off_m+m-1])*(TMath::Abs(q)+TMath::Abs(r)) <
437  eps*(TMath::Abs(p)*(TMath::Abs(pH[off_m_1+m-1])+TMath::Abs(z)+
438  TMath::Abs(pH[off_m1+m+1]))))
439  break;
440  m--;
441  }
442 
443  for (i = m+2; i <= n; i++) {
444  const Int_t off_i = i*nn;
445  pH[off_i+i-2] = 0.0;
446  if (i > m+2)
447  pH[off_i+i-3] = 0.0;
448  }
449 
450  // Double QR step involving rows l:n and columns m:n
451 
452  for (k = m; k <= n-1; k++) {
453  const Int_t off_k = k*nn;
454  const Int_t off_k1 = (k+1)*nn;
455  const Int_t off_k2 = (k+2)*nn;
456  const Int_t notlast = (k != n-1);
457  if (k != m) {
458  p = pH[off_k+k-1];
459  q = pH[off_k1+k-1];
460  r = (notlast ? pH[off_k2+k-1] : 0.0);
461  x = TMath::Abs(p)+TMath::Abs(q)+TMath::Abs(r);
462  if (x != 0.0) {
463  p = p/x;
464  q = q/x;
465  r = r/x;
466  }
467  }
468  if (x == 0.0)
469  break;
470  s = TMath::Sqrt(p*p+q*q+r*r);
471  if (p < 0) {
472  s = -s;
473  }
474  if (s != 0) {
475  if (k != m)
476  pH[off_k+k-1] = -s*x;
477  else if (l != m)
478  pH[off_k+k-1] = -pH[off_k+k-1];
479  p = p+s;
480  x = p/s;
481  y = q/s;
482  z = r/s;
483  q = q/p;
484  r = r/p;
485 
486  // Row modification
487 
488  for (j = k; j < nn; j++) {
489  p = pH[off_k+j]+q*pH[off_k1+j];
490  if (notlast) {
491  p = p+r*pH[off_k2+j];
492  pH[off_k2+j] = pH[off_k2+j]-p*z;
493  }
494  pH[off_k+j] = pH[off_k+j]-p*x;
495  pH[off_k1+j] = pH[off_k1+j]-p*y;
496  }
497 
498  // Column modification
499 
500  for (i = 0; i <= TMath::Min(n,k+3); i++) {
501  const Int_t off_i = i*nn;
502  p = x*pH[off_i+k]+y*pH[off_i+k+1];
503  if (notlast) {
504  p = p+z*pH[off_i+k+2];
505  pH[off_i+k+2] = pH[off_i+k+2]-p*r;
506  }
507  pH[off_i+k] = pH[off_i+k]-p;
508  pH[off_i+k+1] = pH[off_i+k+1]-p*q;
509  }
510 
511  // Accumulate transformations
512 
513  for (i = low; i <= high; i++) {
514  const Int_t off_i = i*nn;
515  p = x*pV[off_i+k]+y*pV[off_i+k+1];
516  if (notlast) {
517  p = p+z*pV[off_i+k+2];
518  pV[off_i+k+2] = pV[off_i+k+2]-p*r;
519  }
520  pV[off_i+k] = pV[off_i+k]-p;
521  pV[off_i+k+1] = pV[off_i+k+1]-p*q;
522  }
523  } // (s != 0)
524  } // k loop
525  } // check convergence
526  } // while (n >= low)
527 
528  // Backsubstitute to find vectors of upper triangular form
529 
530  if (norm == 0.0)
531  return;
532 
533  for (n = nn-1; n >= 0; n--) {
534  p = pD[n];
535  q = pE[n];
536 
537  // Double_t vector
538 
539  const Int_t off_n = n*nn;
540  if (q == 0) {
541  Int_t l = n;
542  pH[off_n+n] = 1.0;
543  for (i = n-1; i >= 0; i--) {
544  const Int_t off_i = i*nn;
545  const Int_t off_i1 = (i+1)*nn;
546  w = pH[off_i+i]-p;
547  r = 0.0;
548  for (j = l; j <= n; j++) {
549  const Int_t off_j = j*nn;
550  r = r+pH[off_i+j]*pH[off_j+n];
551  }
552  if (pE[i] < 0.0) {
553  z = w;
554  s = r;
555  } else {
556  l = i;
557  if (pE[i] == 0.0) {
558  if (w != 0.0)
559  pH[off_i+n] = -r/w;
560  else
561  pH[off_i+n] = -r/(eps*norm);
562 
563  // Solve real equations
564 
565  } else {
566  x = pH[off_i+i+1];
567  y = pH[off_i1+i];
568  q = (pD[i]-p)*(pD[i]-p)+pE[i]*pE[i];
569  t = (x*s-z*r)/q;
570  pH[off_i+n] = t;
571  if (TMath::Abs(x) > TMath::Abs(z))
572  pH[i+1+n] = (-r-w*t)/x;
573  else
574  pH[i+1+n] = (-s-y*t)/z;
575  }
576 
577  // Overflow control
578 
579  t = TMath::Abs(pH[off_i+n]);
580  if ((eps*t)*t > 1) {
581  for (j = i; j <= n; j++) {
582  const Int_t off_j = j*nn;
583  pH[off_j+n] = pH[off_j+n]/t;
584  }
585  }
586  }
587  }
588 
589  // Complex vector
590 
591  } else if (q < 0) {
592  Int_t l = n-1;
593  const Int_t off_n1 = (n-1)*nn;
594 
595  // Last vector component imaginary so matrix is triangular
596 
597  if (TMath::Abs(pH[off_n+n-1]) > TMath::Abs(pH[off_n1+n])) {
598  pH[off_n1+n-1] = q/pH[off_n+n-1];
599  pH[off_n1+n] = -(pH[off_n+n]-p)/pH[off_n+n-1];
600  } else {
601  cdiv(0.0,-pH[off_n1+n],pH[off_n1+n-1]-p,q);
602  pH[off_n1+n-1] = gCdivr;
603  pH[off_n1+n] = gCdivi;
604  }
605  pH[off_n+n-1] = 0.0;
606  pH[off_n+n] = 1.0;
607  for (i = n-2; i >= 0; i--) {
608  const Int_t off_i = i*nn;
609  const Int_t off_i1 = (i+1)*nn;
610  Double_t ra = 0.0;
611  Double_t sa = 0.0;
612  for (j = l; j <= n; j++) {
613  const Int_t off_j = j*nn;
614  ra += pH[off_i+j]*pH[off_j+n-1];
615  sa += pH[off_i+j]*pH[off_j+n];
616  }
617  w = pH[off_i+i]-p;
618 
619  if (pE[i] < 0.0) {
620  z = w;
621  r = ra;
622  s = sa;
623  } else {
624  l = i;
625  if (pE[i] == 0) {
626  cdiv(-ra,-sa,w,q);
627  pH[off_i+n-1] = gCdivr;
628  pH[off_i+n] = gCdivi;
629  } else {
630 
631  // Solve complex equations
632 
633  x = pH[off_i+i+1];
634  y = pH[off_i1+i];
635  Double_t vr = (pD[i]-p)*(pD[i]-p)+pE[i]*pE[i]-q*q;
636  Double_t vi = (pD[i]-p)*2.0*q;
637  if ((vr == 0.0) && (vi == 0.0)) {
638  vr = eps*norm*(TMath::Abs(w)+TMath::Abs(q)+
639  TMath::Abs(x)+TMath::Abs(y)+TMath::Abs(z));
640  }
641  cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
642  pH[off_i+n-1] = gCdivr;
643  pH[off_i+n] = gCdivi;
644  if (TMath::Abs(x) > (TMath::Abs(z)+TMath::Abs(q))) {
645  pH[off_i1+n-1] = (-ra-w*pH[off_i+n-1]+q*pH[off_i+n])/x;
646  pH[off_i1+n] = (-sa-w*pH[off_i+n]-q*pH[off_i+n-1])/x;
647  } else {
648  cdiv(-r-y*pH[off_i+n-1],-s-y*pH[off_i+n],z,q);
649  pH[off_i1+n-1] = gCdivr;
650  pH[off_i1+n] = gCdivi;
651  }
652  }
653 
654  // Overflow control
655 
656  t = TMath::Max(TMath::Abs(pH[off_i+n-1]),TMath::Abs(pH[off_i+n]));
657  if ((eps*t)*t > 1) {
658  for (j = i; j <= n; j++) {
659  const Int_t off_j = j*nn;
660  pH[off_j+n-1] = pH[off_j+n-1]/t;
661  pH[off_j+n] = pH[off_j+n]/t;
662  }
663  }
664  }
665  }
666  }
667  }
668 
669  // Vectors of isolated roots
670 
671  for (i = 0; i < nn; i++) {
672  if (i < low || i > high) {
673  const Int_t off_i = i*nn;
674  for (j = i; j < nn; j++)
675  pV[off_i+j] = pH[off_i+j];
676  }
677  }
678 
679  // Back transformation to get eigenvectors of original matrix
680 
681  for (j = nn-1; j >= low; j--) {
682  for (i = low; i <= high; i++) {
683  const Int_t off_i = i*nn;
684  z = 0.0;
685  for (k = low; k <= TMath::Min(j,high); k++) {
686  const Int_t off_k = k*nn;
687  z = z+pV[off_i+k]*pH[off_k+j];
688  }
689  pV[off_i+j] = z;
690  }
691  }
692 
693 }
694 
695 ////////////////////////////////////////////////////////////////////////////////
696 /// Sort eigenvalues and corresponding vectors in descending order of Re^2+Im^2
697 /// of the complex eigenvalues .
698 
699 void TMatrixDEigen::Sort(TMatrixD &v,TVectorD &d,TVectorD &e)
700 {
701  // Sort eigenvalues and corresponding vectors.
702  Double_t *pV = v.GetMatrixArray();
703  Double_t *pD = d.GetMatrixArray();
704  Double_t *pE = e.GetMatrixArray();
705 
706  const Int_t n = v.GetNrows();
707 
708  for (Int_t i = 0; i < n-1; i++) {
709  Int_t k = i;
710  Double_t norm = pD[i]*pD[i]+pE[i]*pE[i];
711  Int_t j;
712  for (j = i+1; j < n; j++) {
713  const Double_t norm_new = pD[j]*pD[j]+pE[j]*pE[j];
714  if (norm_new > norm) {
715  k = j;
716  norm = norm_new;
717  }
718  }
719  if (k != i) {
720  Double_t tmp;
721  tmp = pD[k];
722  pD[k] = pD[i];
723  pD[i] = tmp;
724  tmp = pE[k];
725  pE[k] = pE[i];
726  pE[i] = tmp;
727  for (j = 0; j < n; j++) {
728  const Int_t off_j = j*n;
729  tmp = pV[off_j+i];
730  pV[off_j+i] = pV[off_j+k];
731  pV[off_j+k] = tmp;
732  }
733  }
734  }
735 }
736 
737 ////////////////////////////////////////////////////////////////////////////////
738 /// Assignment operator
739 
740 TMatrixDEigen &TMatrixDEigen::operator=(const TMatrixDEigen &source)
741 {
742  if (this != &source) {
743  fEigenVectors.ResizeTo(source.fEigenVectors);
744  fEigenValuesRe.ResizeTo(source.fEigenValuesRe);
745  fEigenValuesIm.ResizeTo(source.fEigenValuesIm);
746  fEigenVectors = source.fEigenVectors;
747  fEigenValuesRe = source.fEigenValuesRe;
748  fEigenValuesIm = source.fEigenValuesIm;
749  }
750  return *this;
751 }
752 
753 ////////////////////////////////////////////////////////////////////////////////
754 /// Computes the block diagonal eigenvalue matrix.
755 /// If the original matrix A is not symmetric, then the eigenvalue
756 /// matrix D is block diagonal with the real eigenvalues in 1-by-1
757 /// blocks and any complex eigenvalues,
758 /// a + i*b, in 2-by-2 blocks, [a, b; -b, a].
759 /// That is, if the complex eigenvalues look like
760 ///
761 /// ~~~
762 /// u + iv . . . . .
763 /// . u - iv . . . .
764 /// . . a + ib . . .
765 /// . . . a - ib . .
766 /// . . . . x .
767 /// . . . . . y
768 /// ~~~
769 ///
770 /// then D looks like
771 ///
772 /// ~~~
773 /// u v . . . .
774 /// -v u . . . .
775 /// . . a b . .
776 /// . . -b a . .
777 /// . . . . x .
778 /// . . . . . y
779 /// ~~~
780 ///
781 /// This keeps V a real matrix in both symmetric and non-symmetric
782 /// cases, and A*V = V*D.
783 ///
784 /// Indexing:
785 /// If matrix A has the index/shape (rowLwb,rowUpb,rowLwb,rowUpb)
786 /// each eigen-vector must have the shape (rowLwb,rowUpb) .
787 /// For convenience, the column index of the eigen-vector matrix
788 /// also runs from rowLwb to rowUpb so that the returned matrix
789 /// has also index/shape (rowLwb,rowUpb,rowLwb,rowUpb) .
790 ///
791 
792 const TMatrixD TMatrixDEigen::GetEigenValues() const
793 {
794  const Int_t nrows = fEigenVectors.GetNrows();
795  const Int_t rowLwb = fEigenVectors.GetRowLwb();
796  const Int_t rowUpb = rowLwb+nrows-1;
797 
798  TMatrixD mD(rowLwb,rowUpb,rowLwb,rowUpb);
799 
800  Double_t *pD = mD.GetMatrixArray();
801  const Double_t * const pd = fEigenValuesRe.GetMatrixArray();
802  const Double_t * const pe = fEigenValuesIm.GetMatrixArray();
803 
804  for (Int_t i = 0; i < nrows; i++) {
805  const Int_t off_i = i*nrows;
806  for (Int_t j = 0; j < nrows; j++)
807  pD[off_i+j] = 0.0;
808  pD[off_i+i] = pd[i];
809  if (pe[i] > 0) {
810  pD[off_i+i+1] = pe[i];
811  } else if (pe[i] < 0) {
812  pD[off_i+i-1] = pe[i];
813  }
814  }
815 
816  return mD;
817 }