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