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 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>
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>
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  label iter = 0;
225 
226  do
227  {
228  iter += 1;
229 
230  // Compute implicit shift
231  cmptType g = EValsRe_[l];
232  cmptType p = (EValsRe_[l + 1] - g)/(2.0*EValsIm_[l]);
233  cmptType r = Foam::hypot(p, cmptType(1));
234 
235  if (p < 0)
236  {
237  r = -r;
238  }
239 
240  EValsRe_[l] = EValsIm_[l]/(p + r);
241  EValsRe_[l + 1] = EValsIm_[l]*(p + r);
242  cmptType dl1 = EValsRe_[l + 1];
243  cmptType h = g - EValsRe_[l];
244 
245  for (label i = l + 2; i < n_; ++i)
246  {
247  EValsRe_[i] -= h;
248  }
249 
250  f += h;
251 
252  // Implicit QL transformation.
253  p = EValsRe_[m];
254  cmptType c = cmptType(1);
255  cmptType c2 = c;
256  cmptType c3 = c;
257  cmptType el1 = EValsIm_[l + 1];
258  cmptType s = Zero;
259  cmptType s2 = Zero;
260 
261  for (label i = m - 1; i >= l; --i)
262  {
263  c3 = c2;
264  c2 = c;
265  s2 = s;
266  g = c*EValsIm_[i];
267  h = c*p;
268  r = std::hypot(p, EValsIm_[i]);
269  EValsIm_[i + 1] = s*r;
270  s = EValsIm_[i]/r;
271  c = p/r;
272  p = c*EValsRe_[i] - s*g;
273  EValsRe_[i + 1] = h + s*(c*g + s*EValsRe_[i]);
274 
275  // Accumulate transformation
276  for (label k = 0; k < n_; ++k)
277  {
278  h = EVecs_[k][i + 1];
279  EVecs_[k][i + 1] = s*EVecs_[k][i] + c*h;
280  EVecs_[k][i] = c*EVecs_[k][i] - s*h;
281  }
282  }
283 
284  p = -s*s2*c3*el1*EValsIm_[l]/dl1;
285  EValsIm_[l] = s*p;
286  EValsRe_[l] = c*p;
287  }
288  while (mag(EValsIm_[l]) > eps*tst1); // Convergence check
289  }
290 
291  EValsRe_[l] = EValsRe_[l] + f;
292  EValsIm_[l] = Zero;
293  }
294 
295  // Sorting eigenvalues and corresponding vectors
296  for (label i = 0; i < n_ - 1; ++i)
297  {
298  label k = i;
299  cmptType p = EValsRe_[i];
300 
301  for (label j = i + 1; j < n_; ++j)
302  {
303  if (EValsRe_[j] < p)
304  {
305  k = j;
306  p = EValsRe_[j];
307  }
308  }
309 
310  if (k != i)
311  {
312  EValsRe_[k] = EValsRe_[i];
313  EValsRe_[i] = p;
314 
315  for (label j = 0; j < n_; ++j)
316  {
317  p = EVecs_[j][i];
318  EVecs_[j][i] = EVecs_[j][k];
319  EVecs_[j][k] = p;
320  }
321  }
322  }
323 }
324 
325 
326 template<class cmptType>
328 {
329  DiagonalMatrix<cmptType> orth_(n_);
330 
331  label low = 0;
332  label high = n_ - 1;
333 
334  for (label m = low + 1; m <= high - 1; ++m)
335  {
336  // Scale column
337  cmptType scale = Zero;
338 
339  for (label i = m; i <= high; ++i)
340  {
341  scale = scale + mag(H_[i][m - 1]);
342  }
343 
344  if (scale != 0.0)
345  {
346  // Compute Householder transformation
347  cmptType h = Zero;
348 
349  for (label i = high; i >= m; --i)
350  {
351  orth_[i] = H_[i][m - 1]/scale;
352  h += orth_[i]*orth_[i];
353  }
354 
355  cmptType g = Foam::sqrt(h);
356 
357  if (orth_[m] > 0)
358  {
359  g = -g;
360  }
361 
362  h -= orth_[m]*g;
363  orth_[m] -= g;
364 
365  // Apply Householder similarity transformation
366  // H = (I - u*u'/h)*H*(I - u*u')/h)
367  for (label j = m; j < n_; ++j)
368  {
369  cmptType f = Zero;
370 
371  for (label i = high; i >= m; --i)
372  {
373  f += orth_[i]*H_[i][j];
374  }
375 
376  f /= h;
377 
378  for (label i = m; i <= high; ++i)
379  {
380  H_[i][j] -= f*orth_[i];
381  }
382  }
383 
384  for (label i = 0; i <= high; ++i)
385  {
386  cmptType f = Zero;
387 
388  for (label j = high; j >= m; --j)
389  {
390  f += orth_[j]*H_[i][j];
391  }
392 
393  f /= h;
394 
395  for (label j = m; j <= high; ++j)
396  {
397  H_[i][j] -= f*orth_[j];
398  }
399  }
400 
401  orth_[m] = scale*orth_[m];
402  H_[m][m-1] = scale*g;
403  }
404  }
405 
406  // Accumulate transformations
407  for (label i = 0; i < n_; ++i)
408  {
409  for (label j = 0; j < n_; ++j)
410  {
411  EVecs_[i][j] = (i == j ? 1.0 : 0.0);
412  }
413  }
414 
415  for (label m = high - 1; m >= low + 1; --m)
416  {
417  if (H_[m][m - 1] != 0.0)
418  {
419  for (label i = m + 1; i <= high; ++i)
420  {
421  orth_[i] = H_[i][m - 1];
422  }
423 
424  for (label j = m; j <= high; ++j)
425  {
426  cmptType g = Zero;
427 
428  for (label i = m; i <= high; ++i)
429  {
430  g += orth_[i]*EVecs_[i][j];
431  }
432 
433  // Double division avoids possible underflow
434  g = (g/orth_[m])/H_[m][m - 1];
435 
436  for (label i = m; i <= high; ++i)
437  {
438  EVecs_[i][j] += g*orth_[i];
439  }
440  }
441  }
442  }
443 }
444 
445 
446 template<class cmptType>
448 {
449  // Initialize
450  label nn = n_;
451  label n = nn - 1;
452  label low = 0;
453  label high = nn - 1;
454  cmptType eps = pow(2.0, -52.0);
455  cmptType exshift = Zero;
456  cmptType p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
457 
458  // Store roots isolated by balance and compute matrix norm
459  cmptType norm = Zero;
460 
461  for (label i = 0; i < nn; ++i)
462  {
463  if ((i < low) || (i > high))
464  {
465  EValsRe_[i] = H_[i][i];
466  EValsIm_[i] = Zero;
467  }
468 
469  for (label j = max(i - 1, 0); j < nn; ++j)
470  {
471  norm += mag(H_[i][j]);
472  }
473  }
474 
475  label iter = 0;
476 
477  while (n >= low)
478  {
479  // Look for single SMALL sub-diagonal element
480  label l = n;
481 
482  while (l > low)
483  {
484  s = mag(H_[l - 1][l - 1]) + mag(H_[l][l]);
485 
486  if (s == 0.0)
487  {
488  s = norm;
489  }
490 
491  if (mag(H_[l][l - 1]) < eps*s)
492  {
493  break;
494  }
495 
496  l--;
497  }
498 
499  // Check for convergence
500  if (l == n) // One root found
501  {
502  H_[n][n] = H_[n][n] + exshift;
503  EValsRe_[n] = H_[n][n];
504  EValsIm_[n] = Zero;
505  --n;
506  iter = 0;
507  }
508  else if (l == n - 1) // Two roots found
509  {
510  w = H_[n][n - 1]*H_[n - 1][n];
511  p = (H_[n - 1][n - 1] - H_[n][n])/2.0;
512  q = p*p + w;
513  z = Foam::sqrt(mag(q));
514  H_[n][n] += exshift;
515  H_[n - 1][n - 1] += exshift;
516  x = H_[n][n];
517 
518  // Scalar pair
519  if (q >= 0)
520  {
521  if (p >= 0)
522  {
523  z = p + z;
524  }
525  else
526  {
527  z = p - z;
528  }
529 
530  EValsRe_[n - 1] = x + z;
531  EValsRe_[n] = EValsRe_[n - 1];
532 
533  if (z != 0.0)
534  {
535  EValsRe_[n] = x - w/z;
536  }
537 
538  EValsIm_[n - 1] = Zero;
539  EValsIm_[n] = Zero;
540  x = H_[n][n - 1];
541  s = mag(x) + mag(z);
542  p = x/s;
543  q = z/s;
544  r = Foam::sqrt(p*p + q*q);
545  p /= r;
546  q /= r;
547 
548  // Row modification
549  for (label j = n - 1; j < nn; ++j)
550  {
551  z = H_[n - 1][j];
552  H_[n - 1][j] = q*z + p*H_[n][j];
553  H_[n][j] = q*H_[n][j] - p*z;
554  }
555 
556  // Column modification
557  for (label i = 0; i <= n; ++i)
558  {
559  z = H_[i][n - 1];
560  H_[i][n - 1] = q*z + p*H_[i][n];
561  H_[i][n] = q*H_[i][n] - p*z;
562  }
563 
564  // Accumulate transformations
565  for (label i = low; i <= high; ++i)
566  {
567  z = EVecs_[i][n - 1];
568  EVecs_[i][n - 1] = q*z + p*EVecs_[i][n];
569  EVecs_[i][n] = q*EVecs_[i][n] - p*z;
570  }
571  }
572  else // Complex pair of roots
573  {
574  EValsRe_[n - 1] = x + p;
575  EValsRe_[n] = x + p;
576  EValsIm_[n - 1] = z;
577  EValsIm_[n] = -z;
578  }
579 
580  n -= 2;
581  iter = 0;
582  }
583  else // No convergence yet
584  {
585  // Form shift
586  x = H_[n][n];
587  y = Zero;
588  w = Zero;
589 
590  if (l < n)
591  {
592  y = H_[n - 1][n - 1];
593  w = H_[n][n - 1]*H_[n - 1][n];
594  }
595 
596  // Wilkinson's original ad-hoc shift
597  if (iter == 10)
598  {
599  exshift += x;
600  for (label i = low; i <= n; ++i)
601  {
602  H_[i][i] -= x;
603  }
604 
605  s = mag(H_[n][n - 1]) + mag(H_[n - 1][n - 2]);
606  x = 0.75*s;
607  y = 0.75*s;
608  w = -0.4375*s*s;
609  }
610 
611  // New ad hoc shift
612  if (iter == 30)
613  {
614  s = (y - x)/2.0;
615  s = s*s + w;
616 
617  if (s > 0)
618  {
619  s = Foam::sqrt(s);
620 
621  if (y < x)
622  {
623  s = -s;
624  }
625 
626  s = x - w/((y - x)/2.0 + s);
627 
628  for (label i = low; i <= n; ++i)
629  {
630  H_[i][i] -= s;
631  }
632 
633  exshift += s;
634  x = 0.964;
635  y = 0.964;
636  w = 0.964;
637  }
638  }
639 
640  iter += 1;
641 
642  // Look for two consecutive SMALL sub-diagonal elements
643  label m = n - 2;
644 
645  while (m >= l)
646  {
647  z = H_[m][m];
648  r = x - z;
649  s = y - z;
650  p = (r*s - w)/H_[m + 1][m] + H_[m][m + 1];
651  q = H_[m + 1][m + 1] - z - r - s;
652  r = H_[m + 2][m + 1];
653  s = mag(p) + mag(q) + mag(r);
654  p /= s;
655  q /= s;
656  r /= s;
657 
658  if (m == l)
659  {
660  break;
661  }
662 
663  if
664  (
665  mag(H_[m][m - 1])*(mag(q) + mag(r))
666  < eps*(mag(p)*(mag(H_[m - 1][m - 1])
667  + mag(z) + mag(H_[m + 1][m + 1])))
668  )
669  {
670  break;
671  }
672 
673  --m;
674  }
675 
676  for (label i = m + 2; i <= n; ++i)
677  {
678  H_[i][i - 2] = Zero;
679 
680  if (i > m + 2)
681  {
682  H_[i][i - 3] = Zero;
683  }
684  }
685 
686  // Double QR step involving rows l:n and columns m:n
687  for (label k = m; k <= n - 1; ++k)
688  {
689  label notlast = (k != n - 1);
690 
691  if (k != m)
692  {
693  p = H_[k][k - 1];
694  q = H_[k + 1][k - 1];
695  r = (notlast ? H_[k + 2][k - 1] : 0.0);
696  x = mag(p) + mag(q) + mag(r);
697 
698  if (x != 0.0)
699  {
700  p /= x;
701  q /= x;
702  r /= x;
703  }
704  }
705 
706  if (x == 0.0)
707  {
708  break;
709  }
710 
711  s = Foam::sqrt(p*p + q*q + r*r);
712 
713  if (p < 0)
714  {
715  s = -s;
716  }
717 
718  if (s != 0)
719  {
720  if (k != m)
721  {
722  H_[k][k - 1] = -s*x;
723  }
724  else if (l != m)
725  {
726  H_[k][k - 1] = -H_[k][k - 1];
727  }
728 
729  p = p + s;
730  x = p/s;
731  y = q/s;
732  z = r/s;
733  q /= p;
734  r /= p;
735 
736  // Row modification
737  for (label j = k; j < nn; ++j)
738  {
739  p = H_[k][j] + q*H_[k + 1][j];
740 
741  if (notlast)
742  {
743  p += r*H_[k + 2][j];
744  H_[k + 2][j] -= p*z;
745  }
746 
747  H_[k][j] -= p*x;
748  H_[k + 1][j] -= p*y;
749  }
750 
751  // Column modification
752  for (label i = 0; i <= min(n, k + 3); ++i)
753  {
754  p = x*H_[i][k] + y*H_[i][k + 1];
755 
756  if (notlast)
757  {
758  p += z*H_[i][k + 2];
759  H_[i][k + 2] -= p*r;
760  }
761 
762  H_[i][k] -= p;
763  H_[i][k + 1] -= p*q;
764  }
765 
766  // Accumulate transformations
767  for (label i = low; i <= high; ++i)
768  {
769  p = x*EVecs_[i][k] + y*EVecs_[i][k + 1];
770 
771  if (notlast)
772  {
773  p += z*EVecs_[i][k + 2];
774  EVecs_[i][k + 2] -= p*r;
775  }
776 
777  EVecs_[i][k] -= p;
778  EVecs_[i][k + 1] -= p*q;
779  }
780  }
781  }
782  }
783  }
784 
785  // Backsubstitute to find vectors of upper triangular form
786  if (norm == 0.0)
787  {
788  return;
789  }
790 
791  for (n = nn-1; n >= 0; --n)
792  {
793  p = EValsRe_[n];
794  q = EValsIm_[n];
795 
796  // scalar vector
797  if (q == 0)
798  {
799  label l = n;
800  H_[n][n] = cmptType(1);
801 
802  for (label i = n-1; i >= 0; --i)
803  {
804  w = H_[i][i] - p;
805  r = Zero;
806 
807  for (label j = l; j <= n; ++j)
808  {
809  r += H_[i][j]*H_[j][n];
810  }
811 
812  if (EValsIm_[i] < 0.0)
813  {
814  z = w;
815  s = r;
816  }
817  else
818  {
819  l = i;
820 
821  if (EValsIm_[i] == 0.0)
822  {
823  if (w != 0.0)
824  {
825  H_[i][n] = -r/w;
826  }
827  else
828  {
829  H_[i][n] = -r/(eps*norm);
830  }
831  }
832  else // Solve real equations
833  {
834  x = H_[i][i + 1];
835  y = H_[i + 1][i];
836  q = (EValsRe_[i] - p)*(EValsRe_[i] - p)
837  + EValsIm_[i]*EValsIm_[i];
838 
839  t = (x*s - z*r)/q;
840  H_[i][n] = t;
841 
842  if (mag(x) > mag(z))
843  {
844  H_[i + 1][n] = (-r - w*t)/x;
845  }
846  else
847  {
848  H_[i + 1][n] = (-s - y*t)/z;
849  }
850  }
851 
852  // Overflow control
853  t = mag(H_[i][n]);
854 
855  if ((eps*t)*t > 1)
856  {
857  for (label j = i; j <= n; ++j)
858  {
859  H_[j][n] /= t;
860  }
861  }
862  }
863  }
864  }
865  else if (q < 0) // Complex vector
866  {
867  label l = n - 1;
868 
869  // Last vector component imaginary so matrix is triangular
870  if (mag(H_[n][n - 1]) > mag(H_[n - 1][n]))
871  {
872  H_[n - 1][n - 1] = q/H_[n][n - 1];
873  H_[n - 1][n] = -(H_[n][n] - p)/H_[n][n - 1];
874  }
875  else
876  {
877  complex cDiv = complex(0.0, -H_[n - 1][n])
878  /complex(H_[n - 1][n-1]-p, q);
879 
880  H_[n - 1][n - 1] = cDiv.Re();
881  H_[n - 1][n] = cDiv.Im();
882  }
883 
884  H_[n][n - 1] = Zero;
885  H_[n][n] = cmptType(1);
886 
887  for (label i = n - 2; i >= 0; --i)
888  {
889  cmptType ra, sa, vr, vi;
890  ra = Zero;
891  sa = Zero;
892 
893  for (label j = l; j <= n; ++j)
894  {
895  ra += H_[i][j]*H_[j][n-1];
896  sa += H_[i][j]*H_[j][n];
897  }
898 
899  w = H_[i][i] - p;
900 
901  if (EValsIm_[i] < 0.0)
902  {
903  z = w;
904  r = ra;
905  s = sa;
906  }
907  else
908  {
909  l = i;
910 
911  if (EValsIm_[i] == 0)
912  {
913  complex cDiv = complex(-ra, -sa)/complex(w, q);
914  H_[i][n - 1] = cDiv.Re();
915  H_[i][n] = cDiv.Im();
916  }
917  else
918  {
919  // Solve complex equations
920  x = H_[i][i + 1];
921  y = H_[i + 1][i];
922  vr = (EValsRe_[i] - p)*(EValsRe_[i] - p)
923  + EValsIm_[i]*EValsIm_[i] - q*q;
924 
925  vi = 2.0*(EValsRe_[i] - p)*q;
926 
927  if ((vr == 0.0) && (vi == 0.0))
928  {
929  vr = eps*norm*(mag(w) + mag(q) + mag(x) + mag(y)
930  + mag(z));
931  }
932 
933  complex cDiv =
934  complex(x*r - z*ra + q*sa, x*s - z*sa - q*ra)
935  /complex(vr, vi);
936 
937  H_[i][n - 1] = cDiv.Re();
938  H_[i][n] = cDiv.Im();
939 
940  if (mag(x) > (mag(z) + mag(q)))
941  {
942  H_[i + 1][n - 1] = (-ra - w*H_[i][n - 1]
943  + q*H_[i][n])/x;
944 
945  H_[i + 1][n] = (-sa - w*H_[i][n]
946  - q*H_[i][n - 1])/x;
947  }
948  else
949  {
950  complex cDiv = complex(-r - y*H_[i][n - 1], -s
951  - y*H_[i][n])/complex(z, q);
952 
953  H_[i + 1][n - 1] = cDiv.Re();
954  H_[i + 1][n] = cDiv.Im();
955  }
956  }
957 
958  // Overflow control
959  t = max(mag(H_[i][n - 1]), mag(H_[i][n]));
960 
961  if ((eps*t)*t > 1)
962  {
963  for (label j = i; j <= n; ++j)
964  {
965  H_[j][n - 1] /= t;
966  H_[j][n] /= t;
967  }
968  }
969  }
970  }
971  }
972  }
973 
974  // Vectors of isolated roots
975  for (label i = 0; i < nn; ++i)
976  {
977  if (i < low || i > high)
978  {
979  for (label j = i; j < nn; ++j)
980  {
981  EVecs_[i][j] = H_[i][j];
982  }
983  }
984  }
985 
986  // Back transformation to get eigenvectors of original matrix
987  for (label j = nn - 1; j >= low; --j)
988  {
989  for (label i = low; i <= high; ++i)
990  {
991  z = Zero;
992 
993  for (label k = low; k <= min(j, high); ++k)
994  {
995  z = z + EVecs_[i][k]*H_[k][j];
996  }
997 
998  EVecs_[i][j] = z;
999  }
1000  }
1001 }
1002 
1003 
1004 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1005 
1006 template<class cmptType>
1008 :
1009  n_(A.n()),
1010  EValsRe_(n_, Zero),
1011  EValsIm_(n_, Zero),
1012  EVecs_(n_, Zero),
1013  H_()
1014 {
1015  if (n_ <= 0)
1016  {
1018  << "Input matrix has zero size."
1019  << abort(FatalError);
1020  }
1021 
1022  if (A.symmetric())
1023  {
1024  EVecs_ = A;
1025  tridiagonaliseSymmMatrix();
1026  symmTridiagonalQL();
1027  }
1028  else
1029  {
1030  H_ = A;
1031  Hessenberg();
1032  realSchur();
1033  }
1034 }
1035 
1036 
1037 template<class cmptType>
1040  const SquareMatrix<cmptType>& A,
1041  bool symmetric
1042 )
1043 :
1044  n_(A.n()),
1045  EValsRe_(n_, Zero),
1046  EValsIm_(n_, Zero),
1047  EVecs_(n_, Zero),
1048  H_()
1049 {
1050  if (n_ <= 0)
1051  {
1053  << "Input matrix has zero size."
1054  << abort(FatalError);
1055  }
1056 
1057  if (symmetric)
1058  {
1059  EVecs_ = A;
1060  tridiagonaliseSymmMatrix();
1061  symmTridiagonalQL();
1062  }
1063  else
1064  {
1065  H_ = A;
1066  Hessenberg();
1067  realSchur();
1068  }
1069 }
1070 
1071 
1072 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1073 
1074 template<class cmptType>
1077 {
1078  SquareMatrix<complex> EVecs(EVecs_.n());
1079 
1081  (
1082  EVecs_.cbegin(),
1083  EVecs_.cend(),
1084  EVecs.begin(),
1085  [&](const cmptType& x){ return complex(x); }
1086  );
1087 
1088  for (label i = 0; i < EValsIm_.size(); ++i)
1089  {
1090  if (mag(EValsIm_[i]) > VSMALL)
1091  {
1092  for (label j = 0; j < EVecs.m(); ++j)
1093  {
1094  EVecs(j, i) = complex(EVecs_(j, i), EVecs_(j, i+1));
1095  EVecs(j, i+1) = EVecs(j, i).conjugate();
1096  }
1097  ++i;
1098  }
1099  }
1100 
1101  return EVecs;
1102 }
1103 
1104 
1105 // ************************************************************************* //
Foam::roots::complex
Definition: Roots.H:57
p
volScalarField & p
Definition: createFieldRefs.H:8
EigenMatrix.H
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::EigenMatrix::EigenMatrix
EigenMatrix()=delete
No default construct.
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:519
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::hypot
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
Definition: dimensionedScalar.C:327
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::constant::universal::h
const dimensionedScalar h
Planck constant.
Definition: setRegionSolidFields.H:33
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::EigenMatrix::complexEVecs
const SquareMatrix< complex > complexEVecs() const
Return right eigenvectors in unpacked form.
Definition: EigenMatrix.C:1076
Foam::FatalError
error FatalError
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:24
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::constant::physicoChemical::c2
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Foam::SquareMatrix< cmptType >
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:381
f
labelList f(nPoints)
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::EigenMatrix
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
Definition: EigenMatrix.H:178
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
y
scalar y
Definition: LISASMDCalcMethod1.H:14