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