Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / QR / SelfAdjointEigenSolver.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_SELFADJOINTEIGENSOLVER_H
26 #define EIGEN_SELFADJOINTEIGENSOLVER_H
27
28 /** \qr_module \ingroup QR_Module
29   * \nonstableyet
30   *
31   * \class SelfAdjointEigenSolver
32   *
33   * \brief Eigen values/vectors solver for selfadjoint matrix
34   *
35   * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
36   *
37   * \note MatrixType must be an actual Matrix type, it can't be an expression type.
38   *
39   * \sa MatrixBase::eigenvalues(), class EigenSolver
40   */
41 template<typename _MatrixType> class SelfAdjointEigenSolver
42 {
43   public:
44
45     enum {Size = _MatrixType::RowsAtCompileTime };
46     typedef _MatrixType MatrixType;
47     typedef typename MatrixType::Scalar Scalar;
48     typedef typename NumTraits<Scalar>::Real RealScalar;
49     typedef std::complex<RealScalar> Complex;
50     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
51     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
52     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
53
54     SelfAdjointEigenSolver()
55         : m_eivec(int(Size), int(Size)),
56           m_eivalues(int(Size))
57     {
58       ei_assert(Size!=Dynamic);
59     }
60
61     SelfAdjointEigenSolver(int size)
62         : m_eivec(size, size),
63           m_eivalues(size)
64     {}
65
66     /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix,
67       * as well as the eigenvectors if \a computeEigenvectors is true.
68       *
69       * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool)
70       */
71     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
72       : m_eivec(matrix.rows(), matrix.cols()),
73         m_eivalues(matrix.cols())
74     {
75       compute(matrix, computeEigenvectors);
76     }
77
78     /** Constructors computing the eigenvalues of the generalized eigen problem
79       * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$
80       * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors
81       * are computed if \a computeEigenvectors is true.
82       *
83       * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool)
84       */
85     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
86       : m_eivec(matA.rows(), matA.cols()),
87         m_eivalues(matA.cols())
88     {
89       compute(matA, matB, computeEigenvectors);
90     }
91
92     void compute(const MatrixType& matrix, bool computeEigenvectors = true);
93
94     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true);
95
96     /** \returns the computed eigen vectors as a matrix of column vectors */
97     MatrixType eigenvectors(void) const
98     {
99       #ifndef NDEBUG
100       ei_assert(m_eigenvectorsOk);
101       #endif
102       return m_eivec;
103     }
104
105     /** \returns the computed eigen values */
106     RealVectorType eigenvalues(void) const { return m_eivalues; }
107
108     /** \returns the positive square root of the matrix
109       *
110       * \note the matrix itself must be positive in order for this to make sense.
111       */
112     MatrixType operatorSqrt() const
113     {
114       return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint();
115     }
116
117     /** \returns the positive inverse square root of the matrix
118       *
119       * \note the matrix itself must be positive definite in order for this to make sense.
120       */
121     MatrixType operatorInverseSqrt() const
122     {
123       return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint();
124     }
125
126
127   protected:
128     MatrixType m_eivec;
129     RealVectorType m_eivalues;
130     #ifndef NDEBUG
131     bool m_eigenvectorsOk;
132     #endif
133 };
134
135 #ifndef EIGEN_HIDE_HEAVY_CODE
136
137 // from Golub's "Matrix Computations", algorithm 5.1.3
138 template<typename Scalar>
139 static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s)
140 {
141   if (b==0)
142   {
143     c = 1; s = 0;
144   }
145   else if (ei_abs(b)>ei_abs(a))
146   {
147     Scalar t = -a/b;
148     s = Scalar(1)/ei_sqrt(1+t*t);
149     c = s * t;
150   }
151   else
152   {
153     Scalar t = -b/a;
154     c = Scalar(1)/ei_sqrt(1+t*t);
155     s = c * t;
156   }
157 }
158
159 /** \internal
160   *
161   * \qr_module
162   *
163   * Performs a QR step on a tridiagonal symmetric matrix represented as a
164   * pair of two vectors \a diag and \a subdiag.
165   *
166   * \param matA the input selfadjoint matrix
167   * \param hCoeffs returned Householder coefficients
168   *
169   * For compilation efficiency reasons, this procedure does not use eigen expression
170   * for its arguments.
171   *
172   * Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
173   * "implicit symmetric QR step with Wilkinson shift"
174   */
175 template<typename RealScalar, typename Scalar>
176 static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n);
177
178 /** Computes the eigenvalues of the selfadjoint matrix \a matrix,
179   * as well as the eigenvectors if \a computeEigenvectors is true.
180   *
181   * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool)
182   */
183 template<typename MatrixType>
184 void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
185 {
186   #ifndef NDEBUG
187   m_eigenvectorsOk = computeEigenvectors;
188   #endif
189   assert(matrix.cols() == matrix.rows());
190   int n = matrix.cols();
191   m_eivalues.resize(n,1);
192
193   if(n==1)
194   {
195     m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
196     m_eivec.setOnes();
197     return;
198   }
199
200   m_eivec = matrix;
201
202   // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?
203   // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
204   // (same for diag and subdiag)
205   RealVectorType& diag = m_eivalues;
206   typename TridiagonalizationType::SubDiagonalType subdiag(n-1);
207   TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors);
208
209   int end = n-1;
210   int start = 0;
211   while (end>0)
212   {
213     for (int i = start; i<end; ++i)
214       if (ei_isMuchSmallerThan(ei_abs(subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
215         subdiag[i] = 0;
216
217     // find the largest unreduced block
218     while (end>0 && subdiag[end-1]==0)
219       end--;
220     if (end<=0)
221       break;
222     start = end - 1;
223     while (start>0 && subdiag[start-1]!=0)
224       start--;
225
226     ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
227   }
228
229   // Sort eigenvalues and corresponding vectors.
230   // TODO make the sort optional ?
231   // TODO use a better sort algorithm !!
232   for (int i = 0; i < n-1; ++i)
233   {
234     int k;
235     m_eivalues.segment(i,n-i).minCoeff(&k);
236     if (k > 0)
237     {
238       std::swap(m_eivalues[i], m_eivalues[k+i]);
239       m_eivec.col(i).swap(m_eivec.col(k+i));
240     }
241   }
242 }
243
244 /** Computes the eigenvalues of the generalized eigen problem
245   * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$
246   * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors
247   * are computed if \a computeEigenvectors is true.
248   *
249   * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool)
250   */
251 template<typename MatrixType>
252 void SelfAdjointEigenSolver<MatrixType>::
253 compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors)
254 {
255   ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
256
257   // Compute the cholesky decomposition of matB = L L'
258   LLT<MatrixType> cholB(matB);
259
260   // compute C = inv(L) A inv(L')
261   MatrixType matC = matA;
262   cholB.matrixL().solveTriangularInPlace(matC);
263   // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' :
264   matC = matC.adjoint().eval();
265   cholB.matrixL().template marked<LowerTriangular>().solveTriangularInPlace(matC);
266   matC = matC.adjoint().eval();
267   // this version works too:
268 //   matC = matC.transpose();
269 //   cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC);
270 //   matC = matC.transpose();
271   // FIXME: this should work: (currently it only does for small matrices)
272 //   Transpose<MatrixType> trMatC(matC);
273 //   cholB.matrixL().conjugate().eval().template marked<LowerTriangular>().solveTriangularInPlace(trMatC);
274
275   compute(matC, computeEigenvectors);
276
277   if (computeEigenvectors)
278   {
279     // transform back the eigen vectors: evecs = inv(U) * evecs
280     cholB.matrixL().adjoint().template marked<UpperTriangular>().solveTriangularInPlace(m_eivec);
281     for (int i=0; i<m_eivec.cols(); ++i)
282       m_eivec.col(i) = m_eivec.col(i).normalized();
283   }
284 }
285
286 #endif // EIGEN_HIDE_HEAVY_CODE
287
288 /** \qr_module
289   *
290   * \returns a vector listing the eigenvalues of this matrix.
291   */
292 template<typename Derived>
293 inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1>
294 MatrixBase<Derived>::eigenvalues() const
295 {
296   ei_assert(Flags&SelfAdjointBit);
297   return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues();
298 }
299
300 template<typename Derived, bool IsSelfAdjoint>
301 struct ei_operatorNorm_selector
302 {
303   static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
304   operatorNorm(const MatrixBase<Derived>& m)
305   {
306     // FIXME if it is really guaranteed that the eigenvalues are already sorted,
307     // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
308     return m.eigenvalues().cwise().abs().maxCoeff();
309   }
310 };
311
312 template<typename Derived> struct ei_operatorNorm_selector<Derived, false>
313 {
314   static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
315   operatorNorm(const MatrixBase<Derived>& m)
316   {
317     typename Derived::PlainMatrixType m_eval(m);
318     // FIXME if it is really guaranteed that the eigenvalues are already sorted,
319     // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
320     return ei_sqrt(
321              (m_eval*m_eval.adjoint())
322              .template marked<SelfAdjoint>()
323              .eigenvalues()
324              .maxCoeff()
325            );
326   }
327 };
328
329 /** \qr_module
330   *
331   * \returns the matrix norm of this matrix.
332   */
333 template<typename Derived>
334 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
335 MatrixBase<Derived>::operatorNorm() const
336 {
337   return ei_operatorNorm_selector<Derived, Flags&SelfAdjointBit>
338        ::operatorNorm(derived());
339 }
340
341 #ifndef EIGEN_EXTERN_INSTANTIATIONS
342 template<typename RealScalar, typename Scalar>
343 static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n)
344 {
345   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
346   RealScalar e2 = ei_abs2(subdiag[end-1]);
347   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2));
348   RealScalar x = diag[start] - mu;
349   RealScalar z = subdiag[start];
350
351   for (int k = start; k < end; ++k)
352   {
353     RealScalar c, s;
354     ei_givens_rotation(x, z, c, s);
355
356     // do T = G' T G
357     RealScalar sdk = s * diag[k] + c * subdiag[k];
358     RealScalar dkp1 = s * subdiag[k] + c * diag[k+1];
359
360     diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]);
361     diag[k+1] = s * sdk + c * dkp1;
362     subdiag[k] = c * sdk - s * dkp1;
363
364     if (k > start)
365       subdiag[k - 1] = c * subdiag[k-1] - s * z;
366
367     x = subdiag[k];
368
369     if (k < end - 1)
370     {
371       z = -s * subdiag[k+1];
372       subdiag[k + 1] = c * subdiag[k+1];
373     }
374
375     // apply the givens rotation to the unit matrix Q = Q * G
376     // G only modifies the two columns k and k+1
377     if (matrixQ)
378     {
379       #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
380       #else
381       int kn = k*n;
382       int kn1 = (k+1)*n;
383       #endif
384       // let's do the product manually to avoid the need of temporaries...
385       for (int i=0; i<n; ++i)
386       {
387         #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
388         Scalar matrixQ_i_k = matrixQ[i*n+k];
389         matrixQ[i*n+k]   = c * matrixQ_i_k - s * matrixQ[i*n+k+1];
390         matrixQ[i*n+k+1] = s * matrixQ_i_k + c * matrixQ[i*n+k+1];
391         #else
392         Scalar matrixQ_i_k = matrixQ[i+kn];
393         matrixQ[i+kn]  = c * matrixQ_i_k - s * matrixQ[i+kn1];
394         matrixQ[i+kn1] = s * matrixQ_i_k + c * matrixQ[i+kn1];
395         #endif
396       }
397     }
398   }
399 }
400 #endif
401
402 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H