SelfAdjointEigenSolver.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
00012 #define EIGEN_SELFADJOINTEIGENSOLVER_H
00013 
00014 #include "./Tridiagonalization.h"
00015 
00016 namespace Eigen { 
00017 
00018 template<typename _MatrixType>
00019 class GeneralizedSelfAdjointEigenSolver;
00020 
00021 namespace internal {
00022 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
00023 }
00024 
00068 template<typename _MatrixType> class SelfAdjointEigenSolver
00069 {
00070   public:
00071 
00072     typedef _MatrixType MatrixType;
00073     enum {
00074       Size = MatrixType::RowsAtCompileTime,
00075       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00076       Options = MatrixType::Options,
00077       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00078     };
00079     
00081     typedef typename MatrixType::Scalar Scalar;
00082     typedef typename MatrixType::Index Index;
00083 
00090     typedef typename NumTraits<Scalar>::Real RealScalar;
00091     
00092     friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
00093 
00099     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
00100     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
00101 
00112     SelfAdjointEigenSolver()
00113         : m_eivec(),
00114           m_eivalues(),
00115           m_subdiag(),
00116           m_isInitialized(false)
00117     { }
00118 
00131     SelfAdjointEigenSolver(Index size)
00132         : m_eivec(size, size),
00133           m_eivalues(size),
00134           m_subdiag(size > 1 ? size - 1 : 1),
00135           m_isInitialized(false)
00136     {}
00137 
00153     SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
00154       : m_eivec(matrix.rows(), matrix.cols()),
00155         m_eivalues(matrix.cols()),
00156         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00157         m_isInitialized(false)
00158     {
00159       compute(matrix, options);
00160     }
00161 
00192     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
00193     
00208     SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
00209 
00228     const MatrixType& eigenvectors() const
00229     {
00230       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00231       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00232       return m_eivec;
00233     }
00234 
00250     const RealVectorType& eigenvalues() const
00251     {
00252       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00253       return m_eivalues;
00254     }
00255 
00274     MatrixType operatorSqrt() const
00275     {
00276       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00277       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00278       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00279     }
00280 
00299     MatrixType operatorInverseSqrt() const
00300     {
00301       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00302       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00303       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00304     }
00305 
00310     ComputationInfo info() const
00311     {
00312       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00313       return m_info;
00314     }
00315 
00321     static const int m_maxIterations = 30;
00322 
00323     #ifdef EIGEN2_SUPPORT
00324     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
00325       : m_eivec(matrix.rows(), matrix.cols()),
00326         m_eivalues(matrix.cols()),
00327         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00328         m_isInitialized(false)
00329     {
00330       compute(matrix, computeEigenvectors);
00331     }
00332     
00333     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00334         : m_eivec(matA.cols(), matA.cols()),
00335           m_eivalues(matA.cols()),
00336           m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
00337           m_isInitialized(false)
00338     {
00339       static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00340     }
00341     
00342     void compute(const MatrixType& matrix, bool computeEigenvectors)
00343     {
00344       compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00345     }
00346 
00347     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00348     {
00349       compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00350     }
00351     #endif // EIGEN2_SUPPORT
00352 
00353   protected:
00354     MatrixType m_eivec;
00355     RealVectorType m_eivalues;
00356     typename TridiagonalizationType::SubDiagonalType m_subdiag;
00357     ComputationInfo m_info;
00358     bool m_isInitialized;
00359     bool m_eigenvectorsOk;
00360 };
00361 
00378 namespace internal {
00379 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00380 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
00381 }
00382 
00383 template<typename MatrixType>
00384 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00385 ::compute(const MatrixType& matrix, int options)
00386 {
00387   using std::abs;
00388   eigen_assert(matrix.cols() == matrix.rows());
00389   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00390           && (options&EigVecMask)!=EigVecMask
00391           && "invalid option parameter");
00392   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00393   Index n = matrix.cols();
00394   m_eivalues.resize(n,1);
00395 
00396   if(n==1)
00397   {
00398     m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0));
00399     if(computeEigenvectors)
00400       m_eivec.setOnes(n,n);
00401     m_info = Success;
00402     m_isInitialized = true;
00403     m_eigenvectorsOk = computeEigenvectors;
00404     return *this;
00405   }
00406 
00407   // declare some aliases
00408   RealVectorType& diag = m_eivalues;
00409   MatrixType& mat = m_eivec;
00410 
00411   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00412   mat = matrix.template triangularView<Lower>();
00413   RealScalar scale = mat.cwiseAbs().maxCoeff();
00414   if(scale==RealScalar(0)) scale = RealScalar(1);
00415   mat.template triangularView<Lower>() /= scale;
00416   m_subdiag.resize(n-1);
00417   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
00418   
00419   Index end = n-1;
00420   Index start = 0;
00421   Index iter = 0; // total number of iterations
00422 
00423   while (end>0)
00424   {
00425     for (Index i = start; i<end; ++i)
00426       if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1]))))
00427         m_subdiag[i] = 0;
00428 
00429     // find the largest unreduced block
00430     while (end>0 && m_subdiag[end-1]==0)
00431     {
00432       end--;
00433     }
00434     if (end<=0)
00435       break;
00436 
00437     // if we spent too many iterations, we give up
00438     iter++;
00439     if(iter > m_maxIterations * n) break;
00440 
00441     start = end - 1;
00442     while (start>0 && m_subdiag[start-1]!=0)
00443       start--;
00444 
00445     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
00446   }
00447 
00448   if (iter <= m_maxIterations * n)
00449     m_info = Success;
00450   else
00451     m_info = NoConvergence;
00452 
00453   // Sort eigenvalues and corresponding vectors.
00454   // TODO make the sort optional ?
00455   // TODO use a better sort algorithm !!
00456   if (m_info == Success)
00457   {
00458     for (Index i = 0; i < n-1; ++i)
00459     {
00460       Index k;
00461       m_eivalues.segment(i,n-i).minCoeff(&k);
00462       if (k > 0)
00463       {
00464         std::swap(m_eivalues[i], m_eivalues[k+i]);
00465         if(computeEigenvectors)
00466           m_eivec.col(i).swap(m_eivec.col(k+i));
00467       }
00468     }
00469   }
00470   
00471   // scale back the eigen values
00472   m_eivalues *= scale;
00473 
00474   m_isInitialized = true;
00475   m_eigenvectorsOk = computeEigenvectors;
00476   return *this;
00477 }
00478 
00479 
00480 namespace internal {
00481   
00482 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
00483 {
00484   static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
00485   { eig.compute(A,options); }
00486 };
00487 
00488 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
00489 {
00490   typedef typename SolverType::MatrixType MatrixType;
00491   typedef typename SolverType::RealVectorType VectorType;
00492   typedef typename SolverType::Scalar Scalar;
00493   
00494   static inline void computeRoots(const MatrixType& m, VectorType& roots)
00495   {
00496     using std::sqrt;
00497     using std::atan2;
00498     using std::cos;
00499     using std::sin;
00500     const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
00501     const Scalar s_sqrt3 = sqrt(Scalar(3.0));
00502 
00503     // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
00504     // eigenvalues are the roots to this equation, all guaranteed to be
00505     // real-valued, because the matrix is symmetric.
00506     Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
00507     Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
00508     Scalar c2 = m(0,0) + m(1,1) + m(2,2);
00509 
00510     // Construct the parameters used in classifying the roots of the equation
00511     // and in solving the equation for the roots in closed form.
00512     Scalar c2_over_3 = c2*s_inv3;
00513     Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
00514     if (a_over_3 > Scalar(0))
00515       a_over_3 = Scalar(0);
00516 
00517     Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
00518 
00519     Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3;
00520     if (q > Scalar(0))
00521       q = Scalar(0);
00522 
00523     // Compute the eigenvalues by solving for the roots of the polynomial.
00524     Scalar rho = sqrt(-a_over_3);
00525     Scalar theta = atan2(sqrt(-q),half_b)*s_inv3;
00526     Scalar cos_theta = cos(theta);
00527     Scalar sin_theta = sin(theta);
00528     roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
00529     roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
00530     roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
00531 
00532     // Sort in increasing order.
00533     if (roots(0) >= roots(1))
00534       std::swap(roots(0),roots(1));
00535     if (roots(1) >= roots(2))
00536     {
00537       std::swap(roots(1),roots(2));
00538       if (roots(0) >= roots(1))
00539         std::swap(roots(0),roots(1));
00540     }
00541   }
00542   
00543   static inline void run(SolverType& solver, const MatrixType& mat, int options)
00544   {
00545     using std::sqrt;
00546     eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
00547     eigen_assert((options&~(EigVecMask|GenEigMask))==0
00548             && (options&EigVecMask)!=EigVecMask
00549             && "invalid option parameter");
00550     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00551     
00552     MatrixType& eivecs = solver.m_eivec;
00553     VectorType& eivals = solver.m_eivalues;
00554   
00555     // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00556     Scalar scale = mat.cwiseAbs().maxCoeff();
00557     MatrixType scaledMat = mat / scale;
00558 
00559     // compute the eigenvalues
00560     computeRoots(scaledMat,eivals);
00561 
00562     // compute the eigen vectors
00563     if(computeEigenvectors)
00564     {
00565       Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();
00566       safeNorm2 *= safeNorm2;
00567       if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
00568       {
00569         eivecs.setIdentity();
00570       }
00571       else
00572       {
00573         scaledMat = scaledMat.template selfadjointView<Lower>();
00574         MatrixType tmp;
00575         tmp = scaledMat;
00576 
00577         Scalar d0 = eivals(2) - eivals(1);
00578         Scalar d1 = eivals(1) - eivals(0);
00579         int k =  d0 > d1 ? 2 : 0;
00580         d0 = d0 > d1 ? d1 : d0;
00581 
00582         tmp.diagonal().array () -= eivals(k);
00583         VectorType cross;
00584         Scalar n;
00585         n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
00586 
00587         if(n>safeNorm2)
00588           eivecs.col(k) = cross / sqrt(n);
00589         else
00590         {
00591           n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
00592 
00593           if(n>safeNorm2)
00594             eivecs.col(k) = cross / sqrt(n);
00595           else
00596           {
00597             n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
00598 
00599             if(n>safeNorm2)
00600               eivecs.col(k) = cross / sqrt(n);
00601             else
00602             {
00603               // the input matrix and/or the eigenvaues probably contains some inf/NaN,
00604               // => exit
00605               // scale back to the original size.
00606               eivals *= scale;
00607 
00608               solver.m_info = NumericalIssue;
00609               solver.m_isInitialized = true;
00610               solver.m_eigenvectorsOk = computeEigenvectors;
00611               return;
00612             }
00613           }
00614         }
00615 
00616         tmp = scaledMat;
00617         tmp.diagonal().array() -= eivals(1);
00618 
00619         if(d0<=Eigen::NumTraits<Scalar>::epsilon())
00620           eivecs.col(1) = eivecs.col(k).unitOrthogonal();
00621         else
00622         {
00623           n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
00624           if(n>safeNorm2)
00625             eivecs.col(1) = cross / sqrt(n);
00626           else
00627           {
00628             n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
00629             if(n>safeNorm2)
00630               eivecs.col(1) = cross / sqrt(n);
00631             else
00632             {
00633               n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();
00634               if(n>safeNorm2)
00635                 eivecs.col(1) = cross / sqrt(n);
00636               else
00637               {
00638                 // we should never reach this point,
00639                 // if so the last two eigenvalues are likely to ve very closed to each other
00640                 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
00641               }
00642             }
00643           }
00644 
00645           // make sure that eivecs[1] is orthogonal to eivecs[2]
00646           Scalar d = eivecs.col(1).dot(eivecs.col(k));
00647           eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
00648         }
00649 
00650         eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();
00651       }
00652     }
00653     // Rescale back to the original size.
00654     eivals *= scale;
00655     
00656     solver.m_info = Success;
00657     solver.m_isInitialized = true;
00658     solver.m_eigenvectorsOk = computeEigenvectors;
00659   }
00660 };
00661 
00662 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel
00663 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false>
00664 {
00665   typedef typename SolverType::MatrixType MatrixType;
00666   typedef typename SolverType::RealVectorType VectorType;
00667   typedef typename SolverType::Scalar Scalar;
00668   
00669   static inline void computeRoots(const MatrixType& m, VectorType& roots)
00670   {
00671     using std::sqrt;
00672     const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
00673     const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
00674     roots(0) = t1 - t0;
00675     roots(1) = t1 + t0;
00676   }
00677   
00678   static inline void run(SolverType& solver, const MatrixType& mat, int options)
00679   {
00680     using std::sqrt;
00681     eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
00682     eigen_assert((options&~(EigVecMask|GenEigMask))==0
00683             && (options&EigVecMask)!=EigVecMask
00684             && "invalid option parameter");
00685     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00686     
00687     MatrixType& eivecs = solver.m_eivec;
00688     VectorType& eivals = solver.m_eivalues;
00689   
00690     // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00691     Scalar scale = mat.cwiseAbs().maxCoeff();
00692     scale = (std::max)(scale,Scalar(1));
00693     MatrixType scaledMat = mat / scale;
00694     
00695     // Compute the eigenvalues
00696     computeRoots(scaledMat,eivals);
00697     
00698     // compute the eigen vectors
00699     if(computeEigenvectors)
00700     {
00701       scaledMat.diagonal().array () -= eivals(1);
00702       Scalar a2 = numext::abs2(scaledMat(0,0));
00703       Scalar c2 = numext::abs2(scaledMat(1,1));
00704       Scalar b2 = numext::abs2(scaledMat(1,0));
00705       if(a2>c2)
00706       {
00707         eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
00708         eivecs.col(1) /= sqrt(a2+b2);
00709       }
00710       else
00711       {
00712         eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
00713         eivecs.col(1) /= sqrt(c2+b2);
00714       }
00715 
00716       eivecs.col(0) << eivecs.col(1).unitOrthogonal();
00717     }
00718     
00719     // Rescale back to the original size.
00720     eivals *= scale;
00721     
00722     solver.m_info = Success;
00723     solver.m_isInitialized = true;
00724     solver.m_eigenvectorsOk = computeEigenvectors;
00725   }
00726 };
00727 
00728 }
00729 
00730 template<typename MatrixType>
00731 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00732 ::computeDirect(const MatrixType& matrix, int options)
00733 {
00734   internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
00735   return *this;
00736 }
00737 
00738 namespace internal {
00739 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00740 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
00741 {
00742   using std::abs;
00743   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
00744   RealScalar e = subdiag[end-1];
00745   // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
00746   // underflow thus leading to inf/NaN values when using the following commented code:
00747 //   RealScalar e2 = numext::abs2(subdiag[end-1]);
00748 //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
00749   // This explain the following, somewhat more complicated, version:
00750   RealScalar mu = diag[end];
00751   if(td==0)
00752     mu -= abs(e);
00753   else
00754   {
00755     RealScalar e2 = numext::abs2(subdiag[end-1]);
00756     RealScalar h = numext::hypot(td,e);
00757     if(e2==0)  mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
00758     else       mu -= e2 / (td + (td>0 ? h : -h));
00759   }
00760   
00761   RealScalar x = diag[start] - mu;
00762   RealScalar z = subdiag[start];
00763   for (Index k = start; k < end; ++k)
00764   {
00765     JacobiRotation<RealScalar> rot;
00766     rot.makeGivens(x, z);
00767 
00768     // do T = G' T G
00769     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
00770     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
00771 
00772     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
00773     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
00774     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
00775     
00776 
00777     if (k > start)
00778       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
00779 
00780     x = subdiag[k];
00781 
00782     if (k < end - 1)
00783     {
00784       z = -rot.s() * subdiag[k+1];
00785       subdiag[k + 1] = rot.c() * subdiag[k+1];
00786     }
00787     
00788     // apply the givens rotation to the unit matrix Q = Q * G
00789     if (matrixQ)
00790     {
00791       // FIXME if StorageOrder == RowMajor this operation is not very efficient
00792       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
00793       q.applyOnTheRight(k,k+1,rot);
00794     }
00795   }
00796 }
00797 
00798 } // end namespace internal
00799 
00800 } // end namespace Eigen
00801 
00802 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:52