Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / smoke / intern / EIGENVALUE_HELPER.cpp
1 /** \file \ingroup smoke
2  */
3
4 #include "EIGENVALUE_HELPER.h"
5
6
7 void Eigentred2(sEigenvalue& eval) {
8
9    //  This is derived from the Algol procedures tred2 by
10    //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
11    //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
12    //  Fortran subroutine in EISPACK.
13
14           int n=eval.n;
15
16       for (int j = 0; j < n; j++) {
17          eval.d[j] = eval.V[n-1][j];
18       }
19
20       // Householder reduction to tridiagonal form.
21    
22       for (int i = n-1; i > 0; i--) {
23    
24          // Scale to avoid under/overflow.
25    
26          float scale = 0.0;
27          float h = 0.0;
28          for (int k = 0; k < i; k++) {
29             scale = scale + fabs(eval.d[k]);
30          }
31          if (scale == 0.0f) {
32             eval.e[i] = eval.d[i-1];
33             for (int j = 0; j < i; j++) {
34                eval.d[j] = eval.V[i-1][j];
35                eval.V[i][j] = 0.0;
36                eval.V[j][i] = 0.0;
37             }
38          } else {
39    
40             // Generate Householder vector.
41    
42             for (int k = 0; k < i; k++) {
43                eval.d[k] /= scale;
44                h += eval.d[k] * eval.d[k];
45             }
46             float f = eval.d[i-1];
47             float g = sqrt(h);
48             if (f > 0) {
49                g = -g;
50             }
51             eval.e[i] = scale * g;
52             h = h - f * g;
53             eval.d[i-1] = f - g;
54             for (int j = 0; j < i; j++) {
55                eval.e[j] = 0.0;
56             }
57    
58             // Apply similarity transformation to remaining columns.
59    
60             for (int j = 0; j < i; j++) {
61                f = eval.d[j];
62                eval.V[j][i] = f;
63                g = eval.e[j] + eval.V[j][j] * f;
64                for (int k = j+1; k <= i-1; k++) {
65                   g += eval.V[k][j] * eval.d[k];
66                   eval.e[k] += eval.V[k][j] * f;
67                }
68                eval.e[j] = g;
69             }
70             f = 0.0;
71             for (int j = 0; j < i; j++) {
72                eval.e[j] /= h;
73                f += eval.e[j] * eval.d[j];
74             }
75             float hh = f / (h + h);
76             for (int j = 0; j < i; j++) {
77                eval.e[j] -= hh * eval.d[j];
78             }
79             for (int j = 0; j < i; j++) {
80                f = eval.d[j];
81                g = eval.e[j];
82                for (int k = j; k <= i-1; k++) {
83                   eval.V[k][j] -= (f * eval.e[k] + g * eval.d[k]);
84                }
85                eval.d[j] = eval.V[i-1][j];
86                eval.V[i][j] = 0.0;
87             }
88          }
89          eval.d[i] = h;
90       }
91    
92       // Accumulate transformations.
93    
94       for (int i = 0; i < n-1; i++) {
95          eval.V[n-1][i] = eval.V[i][i];
96          eval.V[i][i] = 1.0;
97          float h = eval.d[i+1];
98          if (h != 0.0f) {
99             for (int k = 0; k <= i; k++) {
100                eval.d[k] = eval.V[k][i+1] / h;
101             }
102             for (int j = 0; j <= i; j++) {
103                float g = 0.0;
104                for (int k = 0; k <= i; k++) {
105                   g += eval.V[k][i+1] * eval.V[k][j];
106                }
107                for (int k = 0; k <= i; k++) {
108                   eval.V[k][j] -= g * eval.d[k];
109                }
110             }
111          }
112          for (int k = 0; k <= i; k++) {
113             eval.V[k][i+1] = 0.0;
114          }
115       }
116       for (int j = 0; j < n; j++) {
117          eval.d[j] = eval.V[n-1][j];
118          eval.V[n-1][j] = 0.0;
119       }
120       eval.V[n-1][n-1] = 1.0;
121       eval.e[0] = 0.0;
122 }
123
124 void Eigencdiv(sEigenvalue& eval, float xr, float xi, float yr, float yi) {
125       float r,d;
126       if (fabs(yr) > fabs(yi)) {
127          r = yi/yr;
128          d = yr + r*yi;
129          eval.cdivr = (xr + r*xi)/d;
130          eval.cdivi = (xi - r*xr)/d;
131       } else {
132          r = yr/yi;
133          d = yi + r*yr;
134          eval.cdivr = (r*xr + xi)/d;
135          eval.cdivi = (r*xi - xr)/d;
136       }
137    }
138
139 void Eigentql2 (sEigenvalue& eval) {
140
141    //  This is derived from the Algol procedures tql2, by
142    //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
143    //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
144    //  Fortran subroutine in EISPACK.
145    
146           int n=eval.n;
147
148       for (int i = 1; i < n; i++) {
149          eval.e[i-1] = eval.e[i];
150       }
151       eval.e[n-1] = 0.0;
152    
153       float f = 0.0;
154       float tst1 = 0.0;
155       float eps = pow(2.0,-52.0);
156       for (int l = 0; l < n; l++) {
157
158          // Find small subdiagonal element
159    
160          tst1 = max(tst1,fabs(eval.d[l]) + fabs(eval.e[l]));
161          int m = l;
162
163         // Original while-loop from Java code
164          while (m < n) {
165             if (fabs(eval.e[m]) <= eps*tst1) {
166                break;
167             }
168             m++;
169          }
170
171    
172          // If m == l, d[l] is an eigenvalue,
173          // otherwise, iterate.
174    
175          if (m > l) {
176             int iter = 0;
177             do {
178                iter = iter + 1;  // (Could check iteration count here.)
179    
180                // Compute implicit shift
181
182                float g = eval.d[l];
183                float p = (eval.d[l+1] - g) / (2.0f * eval.e[l]);
184                float r = hypot(p,1.0);
185                if (p < 0) {
186                   r = -r;
187                }
188                eval.d[l] = eval.e[l] / (p + r);
189                eval.d[l+1] = eval.e[l] * (p + r);
190                float dl1 = eval.d[l+1];
191                float h = g - eval.d[l];
192                for (int i = l+2; i < n; i++) {
193                   eval.d[i] -= h;
194                }
195                f = f + h;
196    
197                // Implicit QL transformation.
198    
199                p = eval.d[m];
200                float c = 1.0;
201                float c2 = c;
202                float c3 = c;
203                float el1 = eval.e[l+1];
204                float s = 0.0;
205                float s2 = 0.0;
206                for (int i = m-1; i >= l; i--) {
207                   c3 = c2;
208                   c2 = c;
209                   s2 = s;
210                   g = c * eval.e[i];
211                   h = c * p;
212                   r = hypot(p,eval.e[i]);
213                   eval.e[i+1] = s * r;
214                   s = eval.e[i] / r;
215                   c = p / r;
216                   p = c * eval.d[i] - s * g;
217                   eval.d[i+1] = h + s * (c * g + s * eval.d[i]);
218    
219                   // Accumulate transformation.
220    
221                   for (int k = 0; k < n; k++) {
222                      h = eval.V[k][i+1];
223                      eval.V[k][i+1] = s * eval.V[k][i] + c * h;
224                      eval.V[k][i] = c * eval.V[k][i] - s * h;
225                   }
226                }
227                p = -s * s2 * c3 * el1 * eval.e[l] / dl1;
228                eval.e[l] = s * p;
229                eval.d[l] = c * p;
230    
231                // Check for convergence.
232    
233             } while (fabs(eval.e[l]) > eps*tst1);
234          }
235          eval.d[l] = eval.d[l] + f;
236          eval.e[l] = 0.0;
237       }
238      
239       // Sort eigenvalues and corresponding vectors.
240    
241       for (int i = 0; i < n-1; i++) {
242          int k = i;
243          float p = eval.d[i];
244          for (int j = i+1; j < n; j++) {
245             if (eval.d[j] < p) {
246                k = j;
247                p = eval.d[j];
248             }
249          }
250          if (k != i) {
251             eval.d[k] = eval.d[i];
252             eval.d[i] = p;
253             for (int j = 0; j < n; j++) {
254                p = eval.V[j][i];
255                eval.V[j][i] = eval.V[j][k];
256                eval.V[j][k] = p;
257             }
258          }
259       }
260 }
261
262 void Eigenorthes (sEigenvalue& eval) {
263    
264       //  This is derived from the Algol procedures orthes and ortran,
265       //  by Martin and Wilkinson, Handbook for Auto. Comp.,
266       //  Vol.ii-Linear Algebra, and the corresponding
267       //  Fortran subroutines in EISPACK.
268    
269           int n=eval.n;
270
271       int low = 0;
272       int high = n-1;
273    
274       for (int m = low+1; m <= high-1; m++) {
275    
276          // Scale column.
277    
278          float scale = 0.0;
279          for (int i = m; i <= high; i++) {
280             scale = scale + fabs(eval.H[i][m-1]);
281          }
282          if (scale != 0.0f) {
283    
284             // Compute Householder transformation.
285    
286             float h = 0.0;
287             for (int i = high; i >= m; i--) {
288                eval.ort[i] = eval.H[i][m-1]/scale;
289                h += eval.ort[i] * eval.ort[i];
290             }
291             float g = sqrt(h);
292             if (eval.ort[m] > 0) {
293                g = -g;
294             }
295             h = h - eval.ort[m] * g;
296             eval.ort[m] = eval.ort[m] - g;
297    
298             // Apply Householder similarity transformation
299             // H = (I-u*u'/h)*H*(I-u*u')/h)
300    
301             for (int j = m; j < n; j++) {
302                float f = 0.0;
303                for (int i = high; i >= m; i--) {
304                   f += eval.ort[i]*eval.H[i][j];
305                }
306                f = f/h;
307                for (int i = m; i <= high; i++) {
308                   eval.H[i][j] -= f*eval.ort[i];
309                }
310            }
311    
312            for (int i = 0; i <= high; i++) {
313                float f = 0.0;
314                for (int j = high; j >= m; j--) {
315                   f += eval.ort[j]*eval.H[i][j];
316                }
317                f = f/h;
318                for (int j = m; j <= high; j++) {
319                   eval.H[i][j] -= f*eval.ort[j];
320                }
321             }
322             eval.ort[m] = scale*eval.ort[m];
323             eval.H[m][m-1] = scale*g;
324          }
325       }
326    
327       // Accumulate transformations (Algol's ortran).
328
329       for (int i = 0; i < n; i++) {
330          for (int j = 0; j < n; j++) {
331             eval.V[i][j] = (i == j ? 1.0 : 0.0);
332          }
333       }
334
335       for (int m = high-1; m >= low+1; m--) {
336          if (eval.H[m][m-1] != 0.0f) {
337             for (int i = m+1; i <= high; i++) {
338                eval.ort[i] = eval.H[i][m-1];
339             }
340             for (int j = m; j <= high; j++) {
341                float g = 0.0;
342                for (int i = m; i <= high; i++) {
343                   g += eval.ort[i] * eval.V[i][j];
344                }
345                // Double division avoids possible underflow
346                g = (g / eval.ort[m]) / eval.H[m][m-1];
347                for (int i = m; i <= high; i++) {
348                   eval.V[i][j] += g * eval.ort[i];
349                }
350             }
351          }
352       }
353    }
354
355 void Eigenhqr2 (sEigenvalue& eval) {
356    
357       //  This is derived from the Algol procedure hqr2,
358       //  by Martin and Wilkinson, Handbook for Auto. Comp.,
359       //  Vol.ii-Linear Algebra, and the corresponding
360       //  Fortran subroutine in EISPACK.
361    
362       // Initialize
363    
364       int nn = eval.n;
365       int n = nn-1;
366       int low = 0;
367       int high = nn-1;
368       float eps = pow(2.0,-52.0);
369       float exshift = 0.0;
370       float p=0,q=0,r=0,s=0,z=0,t,w,x,y;
371    
372       // Store roots isolated by balanc and compute matrix norm
373    
374       float norm = 0.0;
375       for (int i = 0; i < nn; i++) {
376          if ((i < low) || (i > high)) {
377             eval.d[i] = eval.H[i][i];
378             eval.e[i] = 0.0;
379          }
380          for (int j = max(i-1,0); j < nn; j++) {
381             norm = norm + fabs(eval.H[i][j]);
382          }
383       }
384    
385       // Outer loop over eigenvalue index
386    
387       int iter = 0;
388                 int totIter = 0;
389       while (n >= low) {
390
391                         // NT limit no. of iterations
392                         totIter++;
393                         if(totIter>100) {
394                                 //if(totIter>15) std::cout<<"!!!!iter ABORT !!!!!!! "<<totIter<<"\n"; 
395                                 // NT hack/fix, return large eigenvalues
396                                 for (int i = 0; i < nn; i++) {
397                                         eval.d[i] = 10000.;
398                                         eval.e[i] = 10000.;
399                                 }
400                                 return;
401                         }
402    
403          // Look for single small sub-diagonal element
404    
405          int l = n;
406          while (l > low) {
407             s = fabs(eval.H[l-1][l-1]) + fabs(eval.H[l][l]);
408             if (s == 0.0f) {
409                s = norm;
410             }
411             if (fabs(eval.H[l][l-1]) < eps * s) {
412                break;
413             }
414             l--;
415          }
416        
417          // Check for convergence
418          // One root found
419    
420          if (l == n) {
421             eval.H[n][n] = eval.H[n][n] + exshift;
422             eval.d[n] = eval.H[n][n];
423             eval.e[n] = 0.0;
424             n--;
425             iter = 0;
426    
427          // Two roots found
428    
429          } else if (l == n-1) {
430             w = eval.H[n][n-1] * eval.H[n-1][n];
431             p = (eval.H[n-1][n-1] - eval.H[n][n]) / 2.0f;
432             q = p * p + w;
433             z = sqrt(fabs(q));
434             eval.H[n][n] = eval.H[n][n] + exshift;
435             eval.H[n-1][n-1] = eval.H[n-1][n-1] + exshift;
436             x = eval.H[n][n];
437    
438             // float pair
439    
440             if (q >= 0) {
441                if (p >= 0) {
442                   z = p + z;
443                } else {
444                   z = p - z;
445                }
446                eval.d[n-1] = x + z;
447                eval.d[n] = eval.d[n-1];
448                if (z != 0.0f) {
449                   eval.d[n] = x - w / z;
450                }
451                eval.e[n-1] = 0.0;
452                eval.e[n] = 0.0;
453                x = eval.H[n][n-1];
454                s = fabs(x) + fabs(z);
455                p = x / s;
456                q = z / s;
457                r = sqrt(p * p+q * q);
458                p = p / r;
459                q = q / r;
460    
461                // Row modification
462    
463                for (int j = n-1; j < nn; j++) {
464                   z = eval.H[n-1][j];
465                   eval.H[n-1][j] = q * z + p * eval.H[n][j];
466                   eval.H[n][j] = q * eval.H[n][j] - p * z;
467                }
468    
469                // Column modification
470    
471                for (int i = 0; i <= n; i++) {
472                   z = eval.H[i][n-1];
473                   eval.H[i][n-1] = q * z + p * eval.H[i][n];
474                   eval.H[i][n] = q * eval.H[i][n] - p * z;
475                }
476    
477                // Accumulate transformations
478    
479                for (int i = low; i <= high; i++) {
480                   z = eval.V[i][n-1];
481                   eval.V[i][n-1] = q * z + p * eval.V[i][n];
482                   eval.V[i][n] = q * eval.V[i][n] - p * z;
483                }
484    
485             // Complex pair
486    
487             } else {
488                eval.d[n-1] = x + p;
489                eval.d[n] = x + p;
490                eval.e[n-1] = z;
491                eval.e[n] = -z;
492             }
493             n = n - 2;
494             iter = 0;
495    
496          // No convergence yet
497    
498          } else {
499    
500             // Form shift
501    
502             x = eval.H[n][n];
503             y = 0.0;
504             w = 0.0;
505             if (l < n) {
506                y = eval.H[n-1][n-1];
507                w = eval.H[n][n-1] * eval.H[n-1][n];
508             }
509    
510             // Wilkinson's original ad hoc shift
511    
512             if (iter == 10) {
513                exshift += x;
514                for (int i = low; i <= n; i++) {
515                   eval.H[i][i] -= x;
516                }
517                s = fabs(eval.H[n][n-1]) + fabs(eval.H[n-1][n-2]);
518                x = y = 0.75f * s;
519                w = -0.4375f * s * s;
520             }
521
522             // MATLAB's new ad hoc shift
523
524             if (iter == 30) {
525                 s = (y - x) / 2.0f;
526                 s = s * s + w;
527                 if (s > 0) {
528                     s = sqrt(s);
529                     if (y < x) {
530                        s = -s;
531                     }
532                     s = x - w / ((y - x) / 2.0f + s);
533                     for (int i = low; i <= n; i++) {
534                        eval.H[i][i] -= s;
535                     }
536                     exshift += s;
537                     x = y = w = 0.964;
538                 }
539             }
540    
541             iter = iter + 1;   // (Could check iteration count here.)
542    
543             // Look for two consecutive small sub-diagonal elements
544    
545             int m = n-2;
546             while (m >= l) {
547                z = eval.H[m][m];
548                r = x - z;
549                s = y - z;
550                p = (r * s - w) / eval.H[m+1][m] + eval.H[m][m+1];
551                q = eval.H[m+1][m+1] - z - r - s;
552                r = eval.H[m+2][m+1];
553                s = fabs(p) + fabs(q) + fabs(r);
554                p = p / s;
555                q = q / s;
556                r = r / s;
557                if (m == l) {
558                   break;
559                }
560                if (fabs(eval.H[m][m-1]) * (fabs(q) + fabs(r)) <
561                   eps * (fabs(p) * (fabs(eval.H[m-1][m-1]) + fabs(z) +
562                   fabs(eval.H[m+1][m+1])))) {
563                      break;
564                }
565                m--;
566             }
567    
568             for (int i = m+2; i <= n; i++) {
569                eval.H[i][i-2] = 0.0;
570                if (i > m+2) {
571                   eval.H[i][i-3] = 0.0;
572                }
573             }
574    
575             // Double QR step involving rows l:n and columns m:n
576    
577             for (int k = m; k <= n-1; k++) {
578                int notlast = (k != n-1);
579                if (k != m) {
580                   p = eval.H[k][k-1];
581                   q = eval.H[k+1][k-1];
582                   r = (notlast ? eval.H[k+2][k-1] : 0.0f);
583                   x = fabs(p) + fabs(q) + fabs(r);
584                   if (x != 0.0f) {
585                      p = p / x;
586                      q = q / x;
587                      r = r / x;
588                   }
589                }
590                if (x == 0.0f) {
591                   break;
592                }
593                s = sqrt(p * p + q * q + r * r);
594                if (p < 0) {
595                   s = -s;
596                }
597                if (s != 0) {
598                   if (k != m) {
599                      eval.H[k][k-1] = -s * x;
600                   } else if (l != m) {
601                      eval.H[k][k-1] = -eval.H[k][k-1];
602                   }
603                   p = p + s;
604                   x = p / s;
605                   y = q / s;
606                   z = r / s;
607                   q = q / p;
608                   r = r / p;
609    
610                   // Row modification
611    
612                   for (int j = k; j < nn; j++) {
613                      p = eval.H[k][j] + q * eval.H[k+1][j];
614                      if (notlast) {
615                         p = p + r * eval.H[k+2][j];
616                         eval.H[k+2][j] = eval.H[k+2][j] - p * z;
617                      }
618                      eval.H[k][j] = eval.H[k][j] - p * x;
619                      eval.H[k+1][j] = eval.H[k+1][j] - p * y;
620                   }
621    
622                   // Column modification
623    
624                   for (int i = 0; i <= min(n,k+3); i++) {
625                      p = x * eval.H[i][k] + y * eval.H[i][k+1];
626                      if (notlast) {
627                         p = p + z * eval.H[i][k+2];
628                         eval.H[i][k+2] = eval.H[i][k+2] - p * r;
629                      }
630                      eval.H[i][k] = eval.H[i][k] - p;
631                      eval.H[i][k+1] = eval.H[i][k+1] - p * q;
632                   }
633    
634                   // Accumulate transformations
635    
636                   for (int i = low; i <= high; i++) {
637                      p = x * eval.V[i][k] + y * eval.V[i][k+1];
638                      if (notlast) {
639                         p = p + z * eval.V[i][k+2];
640                         eval.V[i][k+2] = eval.V[i][k+2] - p * r;
641                      }
642                      eval.V[i][k] = eval.V[i][k] - p;
643                      eval.V[i][k+1] = eval.V[i][k+1] - p * q;
644                   }
645                }  // (s != 0)
646             }  // k loop
647          }  // check convergence
648       }  // while (n >= low)
649                 //if(totIter>15) std::cout<<"!!!!iter "<<totIter<<"\n";
650       
651       // Backsubstitute to find vectors of upper triangular form
652
653       if (norm == 0.0f) {
654          return;
655       }
656    
657       for (n = nn-1; n >= 0; n--) {
658          p = eval.d[n];
659          q = eval.e[n];
660    
661          // float vector
662    
663          if (q == 0) {
664             int l = n;
665             eval.H[n][n] = 1.0;
666             for (int i = n-1; i >= 0; i--) {
667                w = eval.H[i][i] - p;
668                r = 0.0;
669                for (int j = l; j <= n; j++) {
670                   r = r + eval.H[i][j] * eval.H[j][n];
671                }
672                if (eval.e[i] < 0.0f) {
673                   z = w;
674                   s = r;
675                } else {
676                   l = i;
677                   if (eval.e[i] == 0.0f) {
678                      if (w != 0.0f) {
679                         eval.H[i][n] = -r / w;
680                      } else {
681                         eval.H[i][n] = -r / (eps * norm);
682                      }
683    
684                   // Solve real equations
685    
686                   } else {
687                      x = eval.H[i][i+1];
688                      y = eval.H[i+1][i];
689                      q = (eval.d[i] - p) * (eval.d[i] - p) + eval.e[i] * eval.e[i];
690                      t = (x * s - z * r) / q;
691                      eval.H[i][n] = t;
692                      if (fabs(x) > fabs(z)) {
693                         eval.H[i+1][n] = (-r - w * t) / x;
694                      } else {
695                         eval.H[i+1][n] = (-s - y * t) / z;
696                      }
697                   }
698    
699                   // Overflow control
700    
701                   t = fabs(eval.H[i][n]);
702                   if ((eps * t) * t > 1) {
703                      for (int j = i; j <= n; j++) {
704                         eval.H[j][n] = eval.H[j][n] / t;
705                      }
706                   }
707                }
708             }
709    
710          // Complex vector
711    
712          } else if (q < 0) {
713             int l = n-1;
714
715             // Last vector component imaginary so matrix is triangular
716    
717             if (fabs(eval.H[n][n-1]) > fabs(eval.H[n-1][n])) {
718                eval.H[n-1][n-1] = q / eval.H[n][n-1];
719                eval.H[n-1][n] = -(eval.H[n][n] - p) / eval.H[n][n-1];
720             } else {
721                Eigencdiv(eval, 0.0,-eval.H[n-1][n],eval.H[n-1][n-1]-p,q);
722                eval.H[n-1][n-1] = eval.cdivr;
723                eval.H[n-1][n] = eval.cdivi;
724             }
725             eval.H[n][n-1] = 0.0;
726             eval.H[n][n] = 1.0;
727             for (int i = n-2; i >= 0; i--) {
728                float ra,sa,vr,vi;
729                ra = 0.0;
730                sa = 0.0;
731                for (int j = l; j <= n; j++) {
732                   ra = ra + eval.H[i][j] * eval.H[j][n-1];
733                   sa = sa + eval.H[i][j] * eval.H[j][n];
734                }
735                w = eval.H[i][i] - p;
736    
737                if (eval.e[i] < 0.0f) {
738                   z = w;
739                   r = ra;
740                   s = sa;
741                } else {
742                   l = i;
743                   if (eval.e[i] == 0) {
744                      Eigencdiv(eval,-ra,-sa,w,q);
745                      eval.H[i][n-1] = eval.cdivr;
746                      eval.H[i][n] = eval.cdivi;
747                   } else {
748    
749                      // Solve complex equations
750    
751                      x = eval.H[i][i+1];
752                      y = eval.H[i+1][i];
753                      vr = (eval.d[i] - p) * (eval.d[i] - p) + eval.e[i] * eval.e[i] - q * q;
754                      vi = (eval.d[i] - p) * 2.0f * q;
755                      if ((vr == 0.0f) && (vi == 0.0f)) {
756                         vr = eps * norm * (fabs(w) + fabs(q) +
757                         fabs(x) + fabs(y) + fabs(z));
758                      }
759                      Eigencdiv(eval, x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
760                      eval.H[i][n-1] = eval.cdivr;
761                      eval.H[i][n] = eval.cdivi;
762                      if (fabs(x) > (fabs(z) + fabs(q))) {
763                         eval.H[i+1][n-1] = (-ra - w * eval.H[i][n-1] + q * eval.H[i][n]) / x;
764                         eval.H[i+1][n] = (-sa - w * eval.H[i][n] - q * eval.H[i][n-1]) / x;
765                      } else {
766                         Eigencdiv(eval, -r-y*eval.H[i][n-1],-s-y*eval.H[i][n],z,q);
767                         eval.H[i+1][n-1] = eval.cdivr;
768                         eval.H[i+1][n] = eval.cdivi;
769                      }
770                   }
771    
772                   // Overflow control
773
774                   t = max(fabs(eval.H[i][n-1]),fabs(eval.H[i][n]));
775                   if ((eps * t) * t > 1) {
776                      for (int j = i; j <= n; j++) {
777                         eval.H[j][n-1] = eval.H[j][n-1] / t;
778                         eval.H[j][n] = eval.H[j][n] / t;
779                      }
780                   }
781                }
782             }
783          }
784       }
785    
786       // Vectors of isolated roots
787    
788       for (int i = 0; i < nn; i++) {
789          if (i < low || i > high) {
790             for (int j = i; j < nn; j++) {
791                eval.V[i][j] = eval.H[i][j];
792             }
793          }
794       }
795    
796       // Back transformation to get eigenvectors of original matrix
797    
798       for (int j = nn-1; j >= low; j--) {
799          for (int i = low; i <= high; i++) {
800             z = 0.0;
801             for (int k = low; k <= min(j,high); k++) {
802                z = z + eval.V[i][k] * eval.H[k][j];
803             }
804             eval.V[i][j] = z;
805          }
806       }
807 }
808
809
810
811 int computeEigenvalues3x3(
812                 float dout[3], 
813                 float a[3][3])
814 {
815   /*TNT::Array2D<float> A = TNT::Array2D<float>(3,3, &a[0][0]);
816   TNT::Array1D<float> eig = TNT::Array1D<float>(3);
817   TNT::Array1D<float> eigImag = TNT::Array1D<float>(3);
818   JAMA::Eigenvalue<float> jeig = JAMA::Eigenvalue<float>(A);*/
819
820         sEigenvalue jeig;
821
822         // Compute the values
823         {
824                 jeig.n = 3;
825                 int n=3;
826       //V = Array2D<float>(n,n);
827       //d = Array1D<float>(n);
828       //e = Array1D<float>(n);
829                 for (int y=0; y<3; y++)
830                  {
831                          jeig.d[y]=0.0f;
832                          jeig.e[y]=0.0f;
833                          for (int t=0; t<3; t++) jeig.V[y][t]=0.0f;
834                  }
835
836       jeig.issymmetric = 1;
837       for (int j = 0; (j < 3) && jeig.issymmetric; j++) {
838          for (int i = 0; (i < 3) && jeig.issymmetric; i++) {
839             jeig.issymmetric = (a[i][j] == a[j][i]);
840          }
841       }
842
843       if (jeig.issymmetric) {
844          for (int i = 0; i < 3; i++) {
845             for (int j = 0; j < 3; j++) {
846                jeig.V[i][j] = a[i][j];
847             }
848          }
849    
850          // Tridiagonalize.
851          Eigentred2(jeig);
852    
853          // Diagonalize.
854          Eigentql2(jeig);
855
856       } else {
857          //H = TNT::Array2D<float>(n,n);
858              for (int y=0; y<3; y++)
859                  {
860                          jeig.ort[y]=0.0f;
861                          for (int t=0; t<3; t++) jeig.H[y][t]=0.0f;
862                  }
863          //ort = TNT::Array1D<float>(n);
864          
865          for (int j = 0; j < n; j++) {
866             for (int i = 0; i < n; i++) {
867                jeig.H[i][j] = a[i][j];
868             }
869          }
870    
871          // Reduce to Hessenberg form.
872          Eigenorthes(jeig);
873    
874          // Reduce Hessenberg to real Schur form.
875          Eigenhqr2(jeig);
876       }
877    }
878
879   //jeig.getfloatEigenvalues(eig);
880
881   // complex ones
882   //jeig.getImagEigenvalues(eigImag);
883   dout[0]  = sqrt(jeig.d[0]*jeig.d[0] + jeig.e[0]*jeig.e[0]);
884   dout[1]  = sqrt(jeig.d[1]*jeig.d[1] + jeig.e[1]*jeig.e[1]);
885   dout[2]  = sqrt(jeig.d[2]*jeig.d[2] + jeig.e[2]*jeig.e[2]);
886   return 0;
887 }