Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / QR / EigenSolver.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_EIGENSOLVER_H
26 #define EIGEN_EIGENSOLVER_H
27
28 /** \ingroup QR_Module
29   * \nonstableyet
30   *
31   * \class EigenSolver
32   *
33   * \brief Eigen values/vectors solver for non selfadjoint matrices
34   *
35   * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
36   *
37   * Currently it only support real matrices.
38   *
39   * \note this code was adapted from JAMA (public domain)
40   *
41   * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver
42   */
43 template<typename _MatrixType> class EigenSolver
44 {
45   public:
46
47     typedef _MatrixType MatrixType;
48     typedef typename MatrixType::Scalar Scalar;
49     typedef typename NumTraits<Scalar>::Real RealScalar;
50     typedef std::complex<RealScalar> Complex;
51     typedef Matrix<Complex, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
52     typedef Matrix<Complex, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> EigenvectorType;
53     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
54     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
55
56     /** 
57     * \brief Default Constructor.
58     *
59     * The default constructor is useful in cases in which the user intends to
60     * perform decompositions via EigenSolver::compute(const MatrixType&).
61     */
62     EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
63
64     EigenSolver(const MatrixType& matrix)
65       : m_eivec(matrix.rows(), matrix.cols()),
66         m_eivalues(matrix.cols()),
67         m_isInitialized(false)
68     {
69       compute(matrix);
70     }
71
72
73     EigenvectorType eigenvectors(void) const;
74
75     /** \returns a real matrix V of pseudo eigenvectors.
76       *
77       * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks,
78       * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D
79       * and V satisfy A*V = V*D.
80       *
81       * More precisely, if the diagonal matrix of the eigen values is:\n
82       * \f$
83       * \left[ \begin{array}{cccccc}
84       * u+iv &      &      &      &   &   \\
85       *      & u-iv &      &      &   &   \\
86       *      &      & a+ib &      &   &   \\
87       *      &      &      & a-ib &   &   \\
88       *      &      &      &      & x &   \\
89       *      &      &      &      &   & y \\
90       * \end{array} \right]
91       * \f$ \n
92       * then, we have:\n
93       * \f$
94       * D =\left[ \begin{array}{cccccc}
95       *  u & v &    &   &   &   \\
96       * -v & u &    &   &   &   \\
97       *    &   &  a & b &   &   \\
98       *    &   & -b & a &   &   \\
99       *    &   &    &   & x &   \\
100       *    &   &    &   &   & y \\
101       * \end{array} \right]
102       * \f$
103       *
104       * \sa pseudoEigenvalueMatrix()
105       */
106     const MatrixType& pseudoEigenvectors() const 
107     { 
108       ei_assert(m_isInitialized && "EigenSolver is not initialized.");
109       return m_eivec; 
110     }
111
112     MatrixType pseudoEigenvalueMatrix() const;
113
114     /** \returns the eigenvalues as a column vector */
115     EigenvalueType eigenvalues() const 
116     { 
117       ei_assert(m_isInitialized && "EigenSolver is not initialized.");
118       return m_eivalues; 
119     }
120
121     void compute(const MatrixType& matrix);
122
123   private:
124
125     void orthes(MatrixType& matH, RealVectorType& ort);
126     void hqr2(MatrixType& matH);
127
128   protected:
129     MatrixType m_eivec;
130     EigenvalueType m_eivalues;
131     bool m_isInitialized;
132 };
133
134 /** \returns the real block diagonal matrix D of the eigenvalues.
135   *
136   * See pseudoEigenvectors() for the details.
137   */
138 template<typename MatrixType>
139 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
140 {
141   ei_assert(m_isInitialized && "EigenSolver is not initialized.");
142   int n = m_eivec.cols();
143   MatrixType matD = MatrixType::Zero(n,n);
144   for (int i=0; i<n; ++i)
145   {
146     if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
147       matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
148     else
149     {
150       matD.template block<2,2>(i,i) <<  ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
151                                        -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
152       ++i;
153     }
154   }
155   return matD;
156 }
157
158 /** \returns the normalized complex eigenvectors as a matrix of column vectors.
159   *
160   * \sa eigenvalues(), pseudoEigenvectors()
161   */
162 template<typename MatrixType>
163 typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigenvectors(void) const
164 {
165   ei_assert(m_isInitialized && "EigenSolver is not initialized.");
166   int n = m_eivec.cols();
167   EigenvectorType matV(n,n);
168   for (int j=0; j<n; ++j)
169   {
170     if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j)))))
171     {
172       // we have a real eigen value
173       matV.col(j) = m_eivec.col(j).template cast<Complex>();
174     }
175     else
176     {
177       // we have a pair of complex eigen values
178       for (int i=0; i<n; ++i)
179       {
180         matV.coeffRef(i,j)   = Complex(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
181         matV.coeffRef(i,j+1) = Complex(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
182       }
183       matV.col(j).normalize();
184       matV.col(j+1).normalize();
185       ++j;
186     }
187   }
188   return matV;
189 }
190
191 template<typename MatrixType>
192 void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
193 {
194   assert(matrix.cols() == matrix.rows());
195   int n = matrix.cols();
196   m_eivalues.resize(n,1);
197
198   MatrixType matH = matrix;
199   RealVectorType ort(n);
200
201   // Reduce to Hessenberg form.
202   orthes(matH, ort);
203
204   // Reduce Hessenberg to real Schur form.
205   hqr2(matH);
206
207   m_isInitialized = true;
208 }
209
210 // Nonsymmetric reduction to Hessenberg form.
211 template<typename MatrixType>
212 void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
213 {
214   //  This is derived from the Algol procedures orthes and ortran,
215   //  by Martin and Wilkinson, Handbook for Auto. Comp.,
216   //  Vol.ii-Linear Algebra, and the corresponding
217   //  Fortran subroutines in EISPACK.
218
219   int n = m_eivec.cols();
220   int low = 0;
221   int high = n-1;
222
223   for (int m = low+1; m <= high-1; ++m)
224   {
225     // Scale column.
226     RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
227     if (scale != 0.0)
228     {
229       // Compute Householder transformation.
230       RealScalar h = 0.0;
231       // FIXME could be rewritten, but this one looks better wrt cache
232       for (int i = high; i >= m; i--)
233       {
234         ort.coeffRef(i) = matH.coeff(i,m-1)/scale;
235         h += ort.coeff(i) * ort.coeff(i);
236       }
237       RealScalar g = ei_sqrt(h);
238       if (ort.coeff(m) > 0)
239         g = -g;
240       h = h - ort.coeff(m) * g;
241       ort.coeffRef(m) = ort.coeff(m) - g;
242
243       // Apply Householder similarity transformation
244       // H = (I-u*u'/h)*H*(I-u*u')/h)
245       int bSize = high-m+1;
246       matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h)
247         * (ort.segment(m, bSize).transpose() *  matH.block(m, m, bSize, n-m)).lazy()).lazy();
248
249       matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)).lazy()
250         * (ort.segment(m, bSize)/h).transpose()).lazy();
251
252       ort.coeffRef(m) = scale*ort.coeff(m);
253       matH.coeffRef(m,m-1) = scale*g;
254     }
255   }
256
257   // Accumulate transformations (Algol's ortran).
258   m_eivec.setIdentity();
259
260   for (int m = high-1; m >= low+1; m--)
261   {
262     if (matH.coeff(m,m-1) != 0.0)
263     {
264       ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
265
266       int bSize = high-m+1;
267       m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) /  (matH.coeff(m,m-1) * ort.coeff(m) ) )
268         * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)).lazy());
269     }
270   }
271 }
272
273 // Complex scalar division.
274 template<typename Scalar>
275 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
276 {
277   Scalar r,d;
278   if (ei_abs(yr) > ei_abs(yi))
279   {
280       r = yi/yr;
281       d = yr + r*yi;
282       return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
283   }
284   else
285   {
286       r = yr/yi;
287       d = yi + r*yr;
288       return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
289   }
290 }
291
292
293 // Nonsymmetric reduction from Hessenberg to real Schur form.
294 template<typename MatrixType>
295 void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
296 {
297   //  This is derived from the Algol procedure hqr2,
298   //  by Martin and Wilkinson, Handbook for Auto. Comp.,
299   //  Vol.ii-Linear Algebra, and the corresponding
300   //  Fortran subroutine in EISPACK.
301
302   // Initialize
303   int nn = m_eivec.cols();
304   int n = nn-1;
305   int low = 0;
306   int high = nn-1;
307   Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
308   Scalar exshift = 0.0;
309   Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y;
310
311   // Store roots isolated by balanc and compute matrix norm
312   // FIXME to be efficient the following would requires a triangular reduxion code
313   // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum();
314   Scalar norm = 0.0;
315   for (int j = 0; j < nn; ++j)
316   {
317     // FIXME what's the purpose of the following since the condition is always false
318     if ((j < low) || (j > high))
319     {
320       m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0);
321     }
322     norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum();
323   }
324
325   // Outer loop over eigenvalue index
326   int iter = 0;
327   while (n >= low)
328   {
329     // Look for single small sub-diagonal element
330     int l = n;
331     while (l > low)
332     {
333       s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l));
334       if (s == 0.0)
335           s = norm;
336       if (ei_abs(matH.coeff(l,l-1)) < eps * s)
337         break;
338       l--;
339     }
340
341     // Check for convergence
342     // One root found
343     if (l == n)
344     {
345       matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
346       m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0);
347       n--;
348       iter = 0;
349     }
350     else if (l == n-1) // Two roots found
351     {
352       w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
353       p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5);
354       q = p * p + w;
355       z = ei_sqrt(ei_abs(q));
356       matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
357       matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift;
358       x = matH.coeff(n,n);
359
360       // Scalar pair
361       if (q >= 0)
362       {
363         if (p >= 0)
364           z = p + z;
365         else
366           z = p - z;
367
368         m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0);
369         m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0);
370
371         x = matH.coeff(n,n-1);
372         s = ei_abs(x) + ei_abs(z);
373         p = x / s;
374         q = z / s;
375         r = ei_sqrt(p * p+q * q);
376         p = p / r;
377         q = q / r;
378
379         // Row modification
380         for (int j = n-1; j < nn; ++j)
381         {
382           z = matH.coeff(n-1,j);
383           matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
384           matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z;
385         }
386
387         // Column modification
388         for (int i = 0; i <= n; ++i)
389         {
390           z = matH.coeff(i,n-1);
391           matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
392           matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z;
393         }
394
395         // Accumulate transformations
396         for (int i = low; i <= high; ++i)
397         {
398           z = m_eivec.coeff(i,n-1);
399           m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
400           m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z;
401         }
402       }
403       else // Complex pair
404       {
405         m_eivalues.coeffRef(n-1) = Complex(x + p, z);
406         m_eivalues.coeffRef(n)   = Complex(x + p, -z);
407       }
408       n = n - 2;
409       iter = 0;
410     }
411     else // No convergence yet
412     {
413       // Form shift
414       x = matH.coeff(n,n);
415       y = 0.0;
416       w = 0.0;
417       if (l < n)
418       {
419           y = matH.coeff(n-1,n-1);
420           w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
421       }
422
423       // Wilkinson's original ad hoc shift
424       if (iter == 10)
425       {
426         exshift += x;
427         for (int i = low; i <= n; ++i)
428           matH.coeffRef(i,i) -= x;
429         s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
430         x = y = Scalar(0.75) * s;
431         w = Scalar(-0.4375) * s * s;
432       }
433
434       // MATLAB's new ad hoc shift
435       if (iter == 30)
436       {
437         s = Scalar((y - x) / 2.0);
438         s = s * s + w;
439         if (s > 0)
440         {
441           s = ei_sqrt(s);
442           if (y < x)
443             s = -s;
444           s = Scalar(x - w / ((y - x) / 2.0 + s));
445           for (int i = low; i <= n; ++i)
446             matH.coeffRef(i,i) -= s;
447           exshift += s;
448           x = y = w = Scalar(0.964);
449         }
450       }
451
452       iter = iter + 1;   // (Could check iteration count here.)
453
454       // Look for two consecutive small sub-diagonal elements
455       int m = n-2;
456       while (m >= l)
457       {
458         z = matH.coeff(m,m);
459         r = x - z;
460         s = y - z;
461         p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1);
462         q = matH.coeff(m+1,m+1) - z - r - s;
463         r = matH.coeff(m+2,m+1);
464         s = ei_abs(p) + ei_abs(q) + ei_abs(r);
465         p = p / s;
466         q = q / s;
467         r = r / s;
468         if (m == l) {
469           break;
470         }
471         if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
472           eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) +
473           ei_abs(matH.coeff(m+1,m+1)))))
474         {
475           break;
476         }
477         m--;
478       }
479
480       for (int i = m+2; i <= n; ++i)
481       {
482         matH.coeffRef(i,i-2) = 0.0;
483         if (i > m+2)
484           matH.coeffRef(i,i-3) = 0.0;
485       }
486
487       // Double QR step involving rows l:n and columns m:n
488       for (int k = m; k <= n-1; ++k)
489       {
490         int notlast = (k != n-1);
491         if (k != m) {
492           p = matH.coeff(k,k-1);
493           q = matH.coeff(k+1,k-1);
494           r = notlast ? matH.coeff(k+2,k-1) : Scalar(0);
495           x = ei_abs(p) + ei_abs(q) + ei_abs(r);
496           if (x != 0.0)
497           {
498             p = p / x;
499             q = q / x;
500             r = r / x;
501           }
502         }
503
504         if (x == 0.0)
505           break;
506
507         s = ei_sqrt(p * p + q * q + r * r);
508
509         if (p < 0)
510           s = -s;
511
512         if (s != 0)
513         {
514           if (k != m)
515             matH.coeffRef(k,k-1) = -s * x;
516           else if (l != m)
517             matH.coeffRef(k,k-1) = -matH.coeff(k,k-1);
518
519           p = p + s;
520           x = p / s;
521           y = q / s;
522           z = r / s;
523           q = q / p;
524           r = r / p;
525
526           // Row modification
527           for (int j = k; j < nn; ++j)
528           {
529             p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
530             if (notlast)
531             {
532               p = p + r * matH.coeff(k+2,j);
533               matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z;
534             }
535             matH.coeffRef(k,j) = matH.coeff(k,j) - p * x;
536             matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y;
537           }
538
539           // Column modification
540           for (int i = 0; i <= std::min(n,k+3); ++i)
541           {
542             p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
543             if (notlast)
544             {
545               p = p + z * matH.coeff(i,k+2);
546               matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r;
547             }
548             matH.coeffRef(i,k) = matH.coeff(i,k) - p;
549             matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q;
550           }
551
552           // Accumulate transformations
553           for (int i = low; i <= high; ++i)
554           {
555             p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
556             if (notlast)
557             {
558               p = p + z * m_eivec.coeff(i,k+2);
559               m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r;
560             }
561             m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p;
562             m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q;
563           }
564         }  // (s != 0)
565       }  // k loop
566     }  // check convergence
567   }  // while (n >= low)
568
569   // Backsubstitute to find vectors of upper triangular form
570   if (norm == 0.0)
571   {
572       return;
573   }
574
575   for (n = nn-1; n >= 0; n--)
576   {
577     p = m_eivalues.coeff(n).real();
578     q = m_eivalues.coeff(n).imag();
579
580     // Scalar vector
581     if (q == 0)
582     {
583       int l = n;
584       matH.coeffRef(n,n) = 1.0;
585       for (int i = n-1; i >= 0; i--)
586       {
587         w = matH.coeff(i,i) - p;
588         r = (matH.row(i).segment(l,n-l+1) * matH.col(n).segment(l, n-l+1))(0,0);
589
590         if (m_eivalues.coeff(i).imag() < 0.0)
591         {
592           z = w;
593           s = r;
594         }
595         else
596         {
597           l = i;
598           if (m_eivalues.coeff(i).imag() == 0.0)
599           {
600             if (w != 0.0)
601               matH.coeffRef(i,n) = -r / w;
602             else
603               matH.coeffRef(i,n) = -r / (eps * norm);
604           }
605           else // Solve real equations
606           {
607             x = matH.coeff(i,i+1);
608             y = matH.coeff(i+1,i);
609             q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
610             t = (x * s - z * r) / q;
611             matH.coeffRef(i,n) = t;
612             if (ei_abs(x) > ei_abs(z))
613               matH.coeffRef(i+1,n) = (-r - w * t) / x;
614             else
615               matH.coeffRef(i+1,n) = (-s - y * t) / z;
616           }
617
618           // Overflow control
619           t = ei_abs(matH.coeff(i,n));
620           if ((eps * t) * t > 1)
621             matH.col(n).end(nn-i) /= t;
622         }
623       }
624     }
625     else if (q < 0) // Complex vector
626     {
627       std::complex<Scalar> cc;
628       int l = n-1;
629
630       // Last vector component imaginary so matrix is triangular
631       if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n)))
632       {
633         matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1);
634         matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1);
635       }
636       else
637       {
638         cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q);
639         matH.coeffRef(n-1,n-1) = ei_real(cc);
640         matH.coeffRef(n-1,n) = ei_imag(cc);
641       }
642       matH.coeffRef(n,n-1) = 0.0;
643       matH.coeffRef(n,n) = 1.0;
644       for (int i = n-2; i >= 0; i--)
645       {
646         Scalar ra,sa,vr,vi;
647         ra = (matH.block(i,l, 1, n-l+1) * matH.block(l,n-1, n-l+1, 1)).lazy()(0,0);
648         sa = (matH.block(i,l, 1, n-l+1) * matH.block(l,n, n-l+1, 1)).lazy()(0,0);
649         w = matH.coeff(i,i) - p;
650
651         if (m_eivalues.coeff(i).imag() < 0.0)
652         {
653           z = w;
654           r = ra;
655           s = sa;
656         }
657         else
658         {
659           l = i;
660           if (m_eivalues.coeff(i).imag() == 0)
661           {
662             cc = cdiv(-ra,-sa,w,q);
663             matH.coeffRef(i,n-1) = ei_real(cc);
664             matH.coeffRef(i,n) = ei_imag(cc);
665           }
666           else
667           {
668             // Solve complex equations
669             x = matH.coeff(i,i+1);
670             y = matH.coeff(i+1,i);
671             vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
672             vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
673             if ((vr == 0.0) && (vi == 0.0))
674               vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z));
675
676             cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
677             matH.coeffRef(i,n-1) = ei_real(cc);
678             matH.coeffRef(i,n) = ei_imag(cc);
679             if (ei_abs(x) > (ei_abs(z) + ei_abs(q)))
680             {
681               matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x;
682               matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x;
683             }
684             else
685             {
686               cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q);
687               matH.coeffRef(i+1,n-1) = ei_real(cc);
688               matH.coeffRef(i+1,n) = ei_imag(cc);
689             }
690           }
691
692           // Overflow control
693           t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n)));
694           if ((eps * t) * t > 1)
695             matH.block(i, n-1, nn-i, 2) /= t;
696
697         }
698       }
699     }
700   }
701
702   // Vectors of isolated roots
703   for (int i = 0; i < nn; ++i)
704   {
705     // FIXME again what's the purpose of this test ?
706     // in this algo low==0 and high==nn-1 !!
707     if (i < low || i > high)
708     {
709       m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i);
710     }
711   }
712
713   // Back transformation to get eigenvectors of original matrix
714   int bRows = high-low+1;
715   for (int j = nn-1; j >= low; j--)
716   {
717     int bSize = std::min(j,high)-low+1;
718     m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize));
719   }
720 }
721
722 #endif // EIGEN_EIGENSOLVER_H