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