EigenMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Code created 2014-2018 by Alberto Passalacqua
9  Contributed 2018-07-31 to the OpenFOAM Foundation
10  Copyright (C) 2018 OpenFOAM Foundation
11  Copyright (C) 2019-2020 Alberto Passalacqua
12  Copyright (C) 2020-2022 OpenCFD Ltd.
13 -------------------------------------------------------------------------------
14 License
15  This file is part of OpenFOAM.
16 
17  OpenFOAM is free software: you can redistribute it and/or modify it
18  under the terms of the GNU General Public License as published by
19  the Free Software Foundation, either version 3 of the License, or
20  (at your option) any later version.
21 
22  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
23  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
24  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25  for more details.
26 
27  You should have received a copy of the GNU General Public License
28  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
29 
30 \*---------------------------------------------------------------------------*/
31 
32 #include "EigenMatrix.H"
33 
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35 
36 template<class cmptType>
37 void Foam::EigenMatrix<cmptType>::tridiagonaliseSymmMatrix()
38 {
39  for (label j = 0; j < n_; ++j)
40  {
41  EValsRe_[j] = EVecs_[n_ - 1][j];
42  }
43 
44  // Householder reduction to tridiagonal form
45  for (label i = n_ - 1; i > 0; --i)
46  {
47  // Scale to avoid under/overflow.
48  cmptType scale = Zero;
49  cmptType h = Zero;
50 
51  for (label k = 0; k < i; ++k)
52  {
53  scale = scale + mag(EValsRe_[k]);
54  }
55 
56  if (scale == 0.0)
57  {
58  EValsIm_[i] = EValsRe_[i - 1];
59 
60  for (label j = 0; j < i; ++j)
61  {
62  EValsRe_[j] = EVecs_[i - 1][j];
63  EVecs_[i][j] = Zero;
64  EVecs_[j][i] = Zero;
65  }
66  }
67  else
68  {
69  // Generate Householder vector
70  for (label k = 0; k < i; ++k)
71  {
72  EValsRe_[k] /= scale;
73  h += EValsRe_[k]*EValsRe_[k];
74  }
75 
76  cmptType f = EValsRe_[i - 1];
77  cmptType g = Foam::sqrt(h);
78 
79  if (f > 0)
80  {
81  g = -g;
82  }
83 
84  EValsIm_[i] = scale*g;
85  h -= f*g;
86  EValsRe_[i - 1] = f - g;
87 
88  for (label j = 0; j < i; ++j)
89  {
90  EValsIm_[j] = Zero;
91  }
92 
93  // Apply similarity transformation to remaining columns
94  for (label j = 0; j < i; ++j)
95  {
96  f = EValsRe_[j];
97  EVecs_[j][i] = f;
98  g = EValsIm_[j] + EVecs_[j][j]*f;
99 
100  for (label k = j + 1; k <= i - 1; ++k)
101  {
102  g += EVecs_[k][j]*EValsRe_[k];
103  EValsIm_[k] += EVecs_[k][j]*f;
104  }
105 
106  EValsIm_[j] = g;
107  }
108 
109  f = Zero;
110 
111  for (label j = 0; j < i; ++j)
112  {
113  EValsIm_[j] /= h;
114  f += EValsIm_[j]*EValsRe_[j];
115  }
116 
117  cmptType hh = f/(2.0*h);
118 
119  for (label j = 0; j < i; ++j)
120  {
121  EValsIm_[j] -= hh*EValsRe_[j];
122  }
123 
124  for (label j = 0; j < i; ++j)
125  {
126  f = EValsRe_[j];
127  g = EValsIm_[j];
128 
129  for (label k = j; k <= i - 1; ++k)
130  {
131  EVecs_[k][j] -=
132  (f*EValsIm_[k] + g*EValsRe_[k]);
133  }
134 
135  EValsRe_[j] = EVecs_[i - 1][j];
136  EVecs_[i][j] = Zero;
137  }
138  }
139 
140  EValsRe_[i] = h;
141  }
142 
143  // Accumulate transformations
144  for (label i = 0; i < n_ - 1; ++i)
145  {
146  EVecs_[n_ - 1][i] = EVecs_[i][i];
147  EVecs_[i][i] = cmptType(1);
148  cmptType h = EValsRe_[i + 1];
149 
150  if (h != 0.0)
151  {
152  for (label k = 0; k <= i; ++k)
153  {
154  EValsRe_[k] = EVecs_[k][i + 1]/h;
155  }
156 
157  for (label j = 0; j <= i; ++j)
158  {
159  cmptType g = Zero;
160 
161  for (label k = 0; k <= i; ++k)
162  {
163  g += EVecs_[k][i + 1]*EVecs_[k][j];
164  }
165 
166  for (label k = 0; k <= i; ++k)
167  {
168  EVecs_[k][j] -= g*EValsRe_[k];
169  }
170  }
171  }
172 
173  for (label k = 0; k <= i; ++k)
174  {
175  EVecs_[k][i + 1] = Zero;
176  }
177  }
178 
179  for (label j = 0; j < n_; ++j)
180  {
181  EValsRe_[j] = EVecs_[n_ - 1][j];
182  EVecs_[n_ - 1][j] = Zero;
183  }
184 
185  EVecs_[n_ - 1][n_ - 1] = cmptType(1);
186  EValsIm_[0] = Zero;
187 }
188 
189 
190 template<class cmptType>
191 void Foam::EigenMatrix<cmptType>::symmTridiagonalQL()
192 {
193  for (label i = 1; i < n_; ++i)
194  {
195  EValsIm_[i - 1] = EValsIm_[i];
196  }
197 
198  EValsIm_[n_-1] = Zero;
199 
200  cmptType f = Zero;
201  cmptType tst1 = Zero;
202  cmptType eps = pow(2.0, -52.0);
203 
204  for (label l = 0; l < n_; l++)
205  {
206  // Find SMALL subdiagonal element
207  tst1 = max(tst1, mag(EValsRe_[l]) + mag(EValsIm_[l]));
208  label m = l;
209 
210  // Original while-loop from Java code
211  while (m < n_)
212  {
213  if (mag(EValsIm_[m]) <= eps*tst1)
214  {
215  break;
216  }
217 
218  ++m;
219  }
220 
221  // If m == l, EValsRe_[l] is an eigenvalue, otherwise, iterate
222  if (m > l)
223  {
224  do
225  {
226  // Compute implicit shift
227  cmptType g = EValsRe_[l];
228  cmptType p = (EValsRe_[l + 1] - g)/(2.0*EValsIm_[l]);
229  cmptType r = Foam::hypot(p, cmptType(1));
230 
231  if (p < 0)
232  {
233  r = -r;
234  }
235 
236  EValsRe_[l] = EValsIm_[l]/(p + r);
237  EValsRe_[l + 1] = EValsIm_[l]*(p + r);
238  cmptType dl1 = EValsRe_[l + 1];
239  cmptType h = g - EValsRe_[l];
240 
241  for (label i = l + 2; i < n_; ++i)
242  {
243  EValsRe_[i] -= h;
244  }
245 
246  f += h;
247 
248  // Implicit QL transformation.
249  p = EValsRe_[m];
250  cmptType c = cmptType(1);
251  cmptType c2 = c;
252  cmptType c3 = c;
253  cmptType el1 = EValsIm_[l + 1];
254  cmptType s = Zero;
255  cmptType s2 = Zero;
256 
257  for (label i = m - 1; i >= l; --i)
258  {
259  c3 = c2;
260  c2 = c;
261  s2 = s;
262  g = c*EValsIm_[i];
263  h = c*p;
264  r = std::hypot(p, EValsIm_[i]);
265  EValsIm_[i + 1] = s*r;
266  s = EValsIm_[i]/r;
267  c = p/r;
268  p = c*EValsRe_[i] - s*g;
269  EValsRe_[i + 1] = h + s*(c*g + s*EValsRe_[i]);
270 
271  // Accumulate transformation
272  for (label k = 0; k < n_; ++k)
273  {
274  h = EVecs_[k][i + 1];
275  EVecs_[k][i + 1] = s*EVecs_[k][i] + c*h;
276  EVecs_[k][i] = c*EVecs_[k][i] - s*h;
277  }
278  }
279 
280  p = -s*s2*c3*el1*EValsIm_[l]/dl1;
281  EValsIm_[l] = s*p;
282  EValsRe_[l] = c*p;
283  }
284  while (mag(EValsIm_[l]) > eps*tst1); // Convergence check
285  }
286 
287  EValsRe_[l] = EValsRe_[l] + f;
288  EValsIm_[l] = Zero;
289  }
290 
291  // Sorting eigenvalues and corresponding vectors
292  for (label i = 0; i < n_ - 1; ++i)
293  {
294  label k = i;
295  cmptType p = EValsRe_[i];
296 
297  for (label j = i + 1; j < n_; ++j)
298  {
299  if (EValsRe_[j] < p)
300  {
301  k = j;
302  p = EValsRe_[j];
303  }
304  }
305 
306  if (k != i)
307  {
308  EValsRe_[k] = EValsRe_[i];
309  EValsRe_[i] = p;
310 
311  for (label j = 0; j < n_; ++j)
312  {
313  p = EVecs_[j][i];
314  EVecs_[j][i] = EVecs_[j][k];
315  EVecs_[j][k] = p;
316  }
317  }
318  }
319 }
320 
321 
322 template<class cmptType>
323 void Foam::EigenMatrix<cmptType>::Hessenberg()
324 {
325  DiagonalMatrix<cmptType> orth_(n_);
326 
327  label low = 0;
328  label high = n_ - 1;
329 
330  for (label m = low + 1; m <= high - 1; ++m)
331  {
332  // Scale column
333  cmptType scale = Zero;
334 
335  for (label i = m; i <= high; ++i)
336  {
337  scale = scale + mag(H_[i][m - 1]);
338  }
339 
340  if (scale != 0.0)
341  {
342  // Compute Householder transformation
343  cmptType h = Zero;
344 
345  for (label i = high; i >= m; --i)
346  {
347  orth_[i] = H_[i][m - 1]/scale;
348  h += orth_[i]*orth_[i];
349  }
350 
351  cmptType g = Foam::sqrt(h);
352 
353  if (orth_[m] > 0)
354  {
355  g = -g;
356  }
357 
358  h -= orth_[m]*g;
359  orth_[m] -= g;
360 
361  // Apply Householder similarity transformation
362  // H = (I - u*u'/h)*H*(I - u*u')/h)
363  for (label j = m; j < n_; ++j)
364  {
365  cmptType f = Zero;
366 
367  for (label i = high; i >= m; --i)
368  {
369  f += orth_[i]*H_[i][j];
370  }
371 
372  f /= h;
373 
374  for (label i = m; i <= high; ++i)
375  {
376  H_[i][j] -= f*orth_[i];
377  }
378  }
379 
380  for (label i = 0; i <= high; ++i)
381  {
382  cmptType f = Zero;
383 
384  for (label j = high; j >= m; --j)
385  {
386  f += orth_[j]*H_[i][j];
387  }
388 
389  f /= h;
390 
391  for (label j = m; j <= high; ++j)
392  {
393  H_[i][j] -= f*orth_[j];
394  }
395  }
396 
397  orth_[m] = scale*orth_[m];
398  H_[m][m-1] = scale*g;
399  }
400  }
401 
402  // Accumulate transformations
403  for (label i = 0; i < n_; ++i)
404  {
405  for (label j = 0; j < n_; ++j)
406  {
407  EVecs_[i][j] = (i == j ? 1.0 : 0.0);
408  }
409  }
410 
411  for (label m = high - 1; m >= low + 1; --m)
412  {
413  if (H_[m][m - 1] != 0.0)
414  {
415  for (label i = m + 1; i <= high; ++i)
416  {
417  orth_[i] = H_[i][m - 1];
418  }
419 
420  for (label j = m; j <= high; ++j)
421  {
422  cmptType g = Zero;
423 
424  for (label i = m; i <= high; ++i)
425  {
426  g += orth_[i]*EVecs_[i][j];
427  }
428 
429  // Double division avoids possible underflow
430  g = (g/orth_[m])/H_[m][m - 1];
431 
432  for (label i = m; i <= high; ++i)
433  {
434  EVecs_[i][j] += g*orth_[i];
435  }
436  }
437  }
438  }
439 }
440 
441 
442 template<class cmptType>
443 void Foam::EigenMatrix<cmptType>::realSchur()
444 {
445  // Initialize
446  label nn = n_;
447  label n = nn - 1;
448  label low = 0;
449  label high = nn - 1;
450  cmptType eps = pow(2.0, -52.0);
451  cmptType exshift = Zero;
452  cmptType p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
453 
454  // Store roots isolated by balance and compute matrix norm
455  cmptType norm = Zero;
456 
457  for (label i = 0; i < nn; ++i)
458  {
459  if ((i < low) || (i > high))
460  {
461  EValsRe_[i] = H_[i][i];
462  EValsIm_[i] = Zero;
463  }
464 
465  for (label j = max(i - 1, 0); j < nn; ++j)
466  {
467  norm += mag(H_[i][j]);
468  }
469  }
470 
471  label iter = 0;
472 
473  while (n >= low)
474  {
475  // Look for single SMALL sub-diagonal element
476  label l = n;
477 
478  while (l > low)
479  {
480  s = mag(H_[l - 1][l - 1]) + mag(H_[l][l]);
481 
482  if (s == 0.0)
483  {
484  s = norm;
485  }
486 
487  if (mag(H_[l][l - 1]) < eps*s)
488  {
489  break;
490  }
491 
492  l--;
493  }
494 
495  // Check for convergence
496  if (l == n) // One root found
497  {
498  H_[n][n] = H_[n][n] + exshift;
499  EValsRe_[n] = H_[n][n];
500  EValsIm_[n] = Zero;
501  --n;
502  iter = 0;
503  }
504  else if (l == n - 1) // Two roots found
505  {
506  w = H_[n][n - 1]*H_[n - 1][n];
507  p = (H_[n - 1][n - 1] - H_[n][n])/2.0;
508  q = p*p + w;
509  z = Foam::sqrt(mag(q));
510  H_[n][n] += exshift;
511  H_[n - 1][n - 1] += exshift;
512  x = H_[n][n];
513 
514  // Scalar pair
515  if (q >= 0)
516  {
517  if (p >= 0)
518  {
519  z = p + z;
520  }
521  else
522  {
523  z = p - z;
524  }
525 
526  EValsRe_[n - 1] = x + z;
527  EValsRe_[n] = EValsRe_[n - 1];
528 
529  if (z != 0.0)
530  {
531  EValsRe_[n] = x - w/z;
532  }
533 
534  EValsIm_[n - 1] = Zero;
535  EValsIm_[n] = Zero;
536  x = H_[n][n - 1];
537  s = mag(x) + mag(z);
538  p = x/s;
539  q = z/s;
540  r = Foam::sqrt(p*p + q*q);
541  p /= r;
542  q /= r;
543 
544  // Row modification
545  for (label j = n - 1; j < nn; ++j)
546  {
547  z = H_[n - 1][j];
548  H_[n - 1][j] = q*z + p*H_[n][j];
549  H_[n][j] = q*H_[n][j] - p*z;
550  }
551 
552  // Column modification
553  for (label i = 0; i <= n; ++i)
554  {
555  z = H_[i][n - 1];
556  H_[i][n - 1] = q*z + p*H_[i][n];
557  H_[i][n] = q*H_[i][n] - p*z;
558  }
559 
560  // Accumulate transformations
561  for (label i = low; i <= high; ++i)
562  {
563  z = EVecs_[i][n - 1];
564  EVecs_[i][n - 1] = q*z + p*EVecs_[i][n];
565  EVecs_[i][n] = q*EVecs_[i][n] - p*z;
566  }
567  }
568  else // Complex pair of roots
569  {
570  EValsRe_[n - 1] = x + p;
571  EValsRe_[n] = x + p;
572  EValsIm_[n - 1] = z;
573  EValsIm_[n] = -z;
574  }
575 
576  n -= 2;
577  iter = 0;
578  }
579  else // No convergence yet
580  {
581  // Form shift
582  x = H_[n][n];
583  y = Zero;
584  w = Zero;
585 
586  if (l < n)
587  {
588  y = H_[n - 1][n - 1];
589  w = H_[n][n - 1]*H_[n - 1][n];
590  }
591 
592  // Wilkinson's original ad-hoc shift
593  if (iter == 10)
594  {
595  exshift += x;
596  for (label i = low; i <= n; ++i)
597  {
598  H_[i][i] -= x;
599  }
600 
601  s = mag(H_[n][n - 1]) + mag(H_[n - 1][n - 2]);
602  x = 0.75*s;
603  y = 0.75*s;
604  w = -0.4375*s*s;
605  }
606 
607  // New ad hoc shift
608  if (iter == 30)
609  {
610  s = (y - x)/2.0;
611  s = s*s + w;
612 
613  if (s > 0)
614  {
615  s = Foam::sqrt(s);
616 
617  if (y < x)
618  {
619  s = -s;
620  }
621 
622  s = x - w/((y - x)/2.0 + s);
623 
624  for (label i = low; i <= n; ++i)
625  {
626  H_[i][i] -= s;
627  }
628 
629  exshift += s;
630  x = 0.964;
631  y = 0.964;
632  w = 0.964;
633  }
634  }
635 
636  iter += 1;
637 
638  // Look for two consecutive SMALL sub-diagonal elements
639  label m = n - 2;
640 
641  while (m >= l)
642  {
643  z = H_[m][m];
644  r = x - z;
645  s = y - z;
646  p = (r*s - w)/H_[m + 1][m] + H_[m][m + 1];
647  q = H_[m + 1][m + 1] - z - r - s;
648  r = H_[m + 2][m + 1];
649  s = mag(p) + mag(q) + mag(r);
650  p /= s;
651  q /= s;
652  r /= s;
653 
654  if (m == l)
655  {
656  break;
657  }
658 
659  if
660  (
661  mag(H_[m][m - 1])*(mag(q) + mag(r))
662  < eps*(mag(p)*(mag(H_[m - 1][m - 1])
663  + mag(z) + mag(H_[m + 1][m + 1])))
664  )
665  {
666  break;
667  }
668 
669  --m;
670  }
671 
672  for (label i = m + 2; i <= n; ++i)
673  {
674  H_[i][i - 2] = Zero;
675 
676  if (i > m + 2)
677  {
678  H_[i][i - 3] = Zero;
679  }
680  }
681 
682  // Double QR step involving rows l:n and columns m:n
683  for (label k = m; k <= n - 1; ++k)
684  {
685  label notlast = (k != n - 1);
686 
687  if (k != m)
688  {
689  p = H_[k][k - 1];
690  q = H_[k + 1][k - 1];
691  r = (notlast ? H_[k + 2][k - 1] : 0.0);
692  x = mag(p) + mag(q) + mag(r);
693 
694  if (x != 0.0)
695  {
696  p /= x;
697  q /= x;
698  r /= x;
699  }
700  }
701 
702  if (x == 0.0)
703  {
704  break;
705  }
706 
707  s = Foam::sqrt(p*p + q*q + r*r);
708 
709  if (p < 0)
710  {
711  s = -s;
712  }
713 
714  if (s != 0)
715  {
716  if (k != m)
717  {
718  H_[k][k - 1] = -s*x;
719  }
720  else if (l != m)
721  {
722  H_[k][k - 1] = -H_[k][k - 1];
723  }
724 
725  p = p + s;
726  x = p/s;
727  y = q/s;
728  z = r/s;
729  q /= p;
730  r /= p;
731 
732  // Row modification
733  for (label j = k; j < nn; ++j)
734  {
735  p = H_[k][j] + q*H_[k + 1][j];
736 
737  if (notlast)
738  {
739  p += r*H_[k + 2][j];
740  H_[k + 2][j] -= p*z;
741  }
742 
743  H_[k][j] -= p*x;
744  H_[k + 1][j] -= p*y;
745  }
746 
747  // Column modification
748  for (label i = 0; i <= min(n, k + 3); ++i)
749  {
750  p = x*H_[i][k] + y*H_[i][k + 1];
751 
752  if (notlast)
753  {
754  p += z*H_[i][k + 2];
755  H_[i][k + 2] -= p*r;
756  }
757 
758  H_[i][k] -= p;
759  H_[i][k + 1] -= p*q;
760  }
761 
762  // Accumulate transformations
763  for (label i = low; i <= high; ++i)
764  {
765  p = x*EVecs_[i][k] + y*EVecs_[i][k + 1];
766 
767  if (notlast)
768  {
769  p += z*EVecs_[i][k + 2];
770  EVecs_[i][k + 2] -= p*r;
771  }
772 
773  EVecs_[i][k] -= p;
774  EVecs_[i][k + 1] -= p*q;
775  }
776  }
777  }
778  }
779  }
780 
781  // Backsubstitute to find vectors of upper triangular form
782  if (norm == 0.0)
783  {
784  return;
785  }
786 
787  for (n = nn-1; n >= 0; --n)
788  {
789  p = EValsRe_[n];
790  q = EValsIm_[n];
791 
792  // scalar vector
793  if (q == 0)
794  {
795  label l = n;
796  H_[n][n] = cmptType(1);
797 
798  for (label i = n-1; i >= 0; --i)
799  {
800  w = H_[i][i] - p;
801  r = Zero;
802 
803  for (label j = l; j <= n; ++j)
804  {
805  r += H_[i][j]*H_[j][n];
806  }
807 
808  if (EValsIm_[i] < 0.0)
809  {
810  z = w;
811  s = r;
812  }
813  else
814  {
815  l = i;
816 
817  if (EValsIm_[i] == 0.0)
818  {
819  if (w != 0.0)
820  {
821  H_[i][n] = -r/w;
822  }
823  else
824  {
825  H_[i][n] = -r/(eps*norm);
826  }
827  }
828  else // Solve real equations
829  {
830  x = H_[i][i + 1];
831  y = H_[i + 1][i];
832  q = (EValsRe_[i] - p)*(EValsRe_[i] - p)
833  + EValsIm_[i]*EValsIm_[i];
834 
835  t = (x*s - z*r)/q;
836  H_[i][n] = t;
837 
838  if (mag(x) > mag(z))
839  {
840  H_[i + 1][n] = (-r - w*t)/x;
841  }
842  else
843  {
844  H_[i + 1][n] = (-s - y*t)/z;
845  }
846  }
847 
848  // Overflow control
849  t = mag(H_[i][n]);
850 
851  if ((eps*t)*t > 1)
852  {
853  for (label j = i; j <= n; ++j)
854  {
855  H_[j][n] /= t;
856  }
857  }
858  }
859  }
860  }
861  else if (q < 0) // Complex vector
862  {
863  label l = n - 1;
864 
865  // Last vector component imaginary so matrix is triangular
866  if (mag(H_[n][n - 1]) > mag(H_[n - 1][n]))
867  {
868  H_[n - 1][n - 1] = q/H_[n][n - 1];
869  H_[n - 1][n] = -(H_[n][n] - p)/H_[n][n - 1];
870  }
871  else
872  {
873  complex cDiv = complex(0.0, -H_[n - 1][n])
874  /complex(H_[n - 1][n-1]-p, q);
875 
876  H_[n - 1][n - 1] = cDiv.real();
877  H_[n - 1][n] = cDiv.imag();
878  }
879 
880  H_[n][n - 1] = Zero;
881  H_[n][n] = cmptType(1);
882 
883  for (label i = n - 2; i >= 0; --i)
884  {
885  cmptType ra, sa, vr, vi;
886  ra = Zero;
887  sa = Zero;
888 
889  for (label j = l; j <= n; ++j)
890  {
891  ra += H_[i][j]*H_[j][n-1];
892  sa += H_[i][j]*H_[j][n];
893  }
894 
895  w = H_[i][i] - p;
896 
897  if (EValsIm_[i] < 0.0)
898  {
899  z = w;
900  r = ra;
901  s = sa;
902  }
903  else
904  {
905  l = i;
906 
907  if (EValsIm_[i] == 0)
908  {
909  complex cDiv = complex(-ra, -sa)/complex(w, q);
910  H_[i][n - 1] = cDiv.real();
911  H_[i][n] = cDiv.imag();
912  }
913  else
914  {
915  // Solve complex equations
916  x = H_[i][i + 1];
917  y = H_[i + 1][i];
918  vr = (EValsRe_[i] - p)*(EValsRe_[i] - p)
919  + EValsIm_[i]*EValsIm_[i] - q*q;
920 
921  vi = 2.0*(EValsRe_[i] - p)*q;
922 
923  if ((vr == 0.0) && (vi == 0.0))
924  {
925  vr = eps*norm*(mag(w) + mag(q) + mag(x) + mag(y)
926  + mag(z));
927  }
928 
929  complex cDiv =
930  complex(x*r - z*ra + q*sa, x*s - z*sa - q*ra)
931  /complex(vr, vi);
932 
933  H_[i][n - 1] = cDiv.real();
934  H_[i][n] = cDiv.imag();
935 
936  if (mag(x) > (mag(z) + mag(q)))
937  {
938  H_[i + 1][n - 1] = (-ra - w*H_[i][n - 1]
939  + q*H_[i][n])/x;
940 
941  H_[i + 1][n] = (-sa - w*H_[i][n]
942  - q*H_[i][n - 1])/x;
943  }
944  else
945  {
946  complex cDiv = complex(-r - y*H_[i][n - 1], -s
947  - y*H_[i][n])/complex(z, q);
948 
949  H_[i + 1][n - 1] = cDiv.real();
950  H_[i + 1][n] = cDiv.imag();
951  }
952  }
953 
954  // Overflow control
955  t = max(mag(H_[i][n - 1]), mag(H_[i][n]));
956 
957  if ((eps*t)*t > 1)
958  {
959  for (label j = i; j <= n; ++j)
960  {
961  H_[j][n - 1] /= t;
962  H_[j][n] /= t;
963  }
964  }
965  }
966  }
967  }
968  }
969 
970  // Vectors of isolated roots
971  for (label i = 0; i < nn; ++i)
972  {
973  if (i < low || i > high)
974  {
975  for (label j = i; j < nn; ++j)
976  {
977  EVecs_[i][j] = H_[i][j];
978  }
979  }
980  }
981 
982  // Back transformation to get eigenvectors of original matrix
983  for (label j = nn - 1; j >= low; --j)
984  {
985  for (label i = low; i <= high; ++i)
986  {
987  z = Zero;
988 
989  for (label k = low; k <= min(j, high); ++k)
990  {
991  z = z + EVecs_[i][k]*H_[k][j];
992  }
993 
994  EVecs_[i][j] = z;
995  }
996  }
997 }
998 
999 
1000 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1001 
1002 template<class cmptType>
1003 Foam::EigenMatrix<cmptType>::EigenMatrix(const SquareMatrix<cmptType>& A)
1004 :
1005  n_(A.n()),
1006  EValsRe_(n_, Zero),
1007  EValsIm_(n_, Zero),
1008  EVecs_(n_, Zero),
1009  H_()
1010 {
1011  if (n_ <= 0)
1012  {
1013  FatalErrorInFunction
1014  << "Input matrix has zero size."
1015  << abort(FatalError);
1016  }
1017 
1018  if (A.symmetric())
1019  {
1020  EVecs_ = A;
1021  tridiagonaliseSymmMatrix();
1022  symmTridiagonalQL();
1023  }
1024  else
1025  {
1026  H_ = A;
1027  Hessenberg();
1028  realSchur();
1029  }
1030 }
1031 
1032 
1033 template<class cmptType>
1034 Foam::EigenMatrix<cmptType>::EigenMatrix
1035 (
1036  const SquareMatrix<cmptType>& A,
1037  bool symmetric
1038 )
1039 :
1040  n_(A.n()),
1041  EValsRe_(n_, Zero),
1042  EValsIm_(n_, Zero),
1043  EVecs_(n_, Zero),
1044  H_()
1045 {
1046  if (n_ <= 0)
1047  {
1048  FatalErrorInFunction
1049  << "Input matrix has zero size."
1050  << abort(FatalError);
1051  }
1052 
1053  if (symmetric)
1054  {
1055  EVecs_ = A;
1056  tridiagonaliseSymmMatrix();
1057  symmTridiagonalQL();
1058  }
1059  else
1060  {
1061  H_ = A;
1062  Hessenberg();
1063  realSchur();
1064  }
1065 }
1066 
1067 
1068 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1069 
1070 template<class cmptType>
1071 const Foam::SquareMatrix<Foam::complex>
1072 Foam::EigenMatrix<cmptType>::complexEVecs() const
1073 {
1074  SquareMatrix<complex> EVecs(EVecs_.n());
1075 
1076  std::transform
1077  (
1078  EVecs_.cbegin(),
1079  EVecs_.cend(),
1080  EVecs.begin(),
1081  [&](const cmptType& x){ return complex(x); }
1082  );
1083 
1084  for (label i = 0; i < EValsIm_.size(); ++i)
1085  {
1086  if (mag(EValsIm_[i]) > VSMALL)
1087  {
1088  for (label j = 0; j < EVecs.m(); ++j)
1089  {
1090  EVecs(j, i) = complex(EVecs_(j, i), EVecs_(j, i+1));
1091  EVecs(j, i+1) = EVecs(j, i).conjugate();
1092  }
1093  ++i;
1094  }
1095  }
1096 
1097  return EVecs;
1098 }
1099 
1100 
1101 // ************************************************************************* //
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::EigenMatrix::EigenMatrix
EigenMatrix()=delete
No default construct.
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:40
Foam::constant::physicoChemical::c2
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:137
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
y
scalar y
Definition: LISASMDCalcMethod1.H:14
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:26
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam::EigenMatrix::complexEVecs
const SquareMatrix< complex > complexEVecs() const
Return right eigenvectors in unpacked form.
Definition: EigenMatrix.C:1065
f
labelList f(nPoints)
Foam::constant::universal::h
const dimensionedScalar h
Planck constant.
Definition: setRegionSolidFields.H:32
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:68
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
n
label n
Definition: TABSMDCalcMethod2.H:31
p
volScalarField & p
Definition: createFieldRefs.H:8
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:521
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::EigenMatrix
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
Definition: EigenMatrix.H:173
Foam::hypot
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
Definition: dimensionedScalar.C:320
Foam::SquareMatrix
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix...
Definition: SquareMatrix.H:59
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127