Logo ROOT   6.16/01
Reference Guide
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
53
54////////////////////////////////////////////////////////////////////////////////
55/// Constructor for eigen-problem of matrix A .
56
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.
86
87 // Reduce Hessenberg to real Schur form.
89
90 // Sort eigenvalues and corresponding vectors in descending order of Re^2+Im^2
91 // of the complex eigenvalues .
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Copy constructor
97
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
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 }
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
213static 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
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;
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];
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);
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)+
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
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
741{
742 if (this != &source) {
746 }
747 return *this;
748}
749
750////////////////////////////////////////////////////////////////////////////////
751/// Computes the block diagonal eigenvalue matrix.
752/// If the original matrix A is not symmetric, then the eigenvalue
753/// matrix D is block diagonal with the real eigenvalues in 1-by-1
754/// blocks and any complex eigenvalues,
755/// a + i*b, in 2-by-2 blocks, [a, b; -b, a].
756/// That is, if the complex eigenvalues look like
757///
758/// ~~~
759/// u + iv . . . . .
760/// . u - iv . . . .
761/// . . a + ib . . .
762/// . . . a - ib . .
763/// . . . . x .
764/// . . . . . y
765/// ~~~
766///
767/// then D looks like
768///
769/// ~~~
770/// u v . . . .
771/// -v u . . . .
772/// . . a b . .
773/// . . -b a . .
774/// . . . . x .
775/// . . . . . y
776/// ~~~
777///
778/// This keeps V a real matrix in both symmetric and non-symmetric
779/// cases, and A*V = V*D.
780///
781/// Indexing:
782/// If matrix A has the index/shape (rowLwb,rowUpb,rowLwb,rowUpb)
783/// each eigen-vector must have the shape (rowLwb,rowUpb) .
784/// For convenience, the column index of the eigen-vector matrix
785/// also runs from rowLwb to rowUpb so that the returned matrix
786/// has also index/shape (rowLwb,rowUpb,rowLwb,rowUpb) .
787///
788
790{
791 const Int_t nrows = fEigenVectors.GetNrows();
792 const Int_t rowLwb = fEigenVectors.GetRowLwb();
793 const Int_t rowUpb = rowLwb+nrows-1;
794
795 TMatrixD mD(rowLwb,rowUpb,rowLwb,rowUpb);
796
797 Double_t *pD = mD.GetMatrixArray();
798 const Double_t * const pd = fEigenValuesRe.GetMatrixArray();
799 const Double_t * const pe = fEigenValuesIm.GetMatrixArray();
800
801 for (Int_t i = 0; i < nrows; i++) {
802 const Int_t off_i = i*nrows;
803 for (Int_t j = 0; j < nrows; j++)
804 pD[off_i+j] = 0.0;
805 pD[off_i+i] = pd[i];
806 if (pe[i] > 0) {
807 pD[off_i+i+1] = pe[i];
808 } else if (pe[i] < 0) {
809 pD[off_i+i-1] = pe[i];
810 }
811 }
812
813 return mD;
814}
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
#define g(i)
Definition: RSha256.hxx:105
#define h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:363
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
float * q
Definition: THbookFile.cxx:87
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.
Definition: TMatrixDEigen.h:27
TVectorD fEigenValuesIm
Definition: TMatrixDEigen.h:36
static void MakeHessenBerg(TMatrixD &v, TVectorD &ortho, TMatrixD &H)
Nonsymmetric reduction to Hessenberg form.
TVectorD fEigenValuesRe
Definition: TMatrixDEigen.h:35
TMatrixD fEigenVectors
Definition: TMatrixDEigen.h:34
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
Definition: TMatrixTBase.h:124
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
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:1210
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:347
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)
static constexpr double s
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12