$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 static void check_template_parameters() 00355 { 00356 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00357 } 00358 00359 MatrixType m_eivec; 00360 RealVectorType m_eivalues; 00361 typename TridiagonalizationType::SubDiagonalType m_subdiag; 00362 ComputationInfo m_info; 00363 bool m_isInitialized; 00364 bool m_eigenvectorsOk; 00365 }; 00366 00383 namespace internal { 00384 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> 00385 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); 00386 } 00387 00388 template<typename MatrixType> 00389 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> 00390 ::compute(const MatrixType& matrix, int options) 00391 { 00392 check_template_parameters(); 00393 00394 using std::abs; 00395 eigen_assert(matrix.cols() == matrix.rows()); 00396 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00397 && (options&EigVecMask)!=EigVecMask 00398 && "invalid option parameter"); 00399 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 00400 Index n = matrix.cols(); 00401 m_eivalues.resize(n,1); 00402 00403 if(n==1) 00404 { 00405 m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); 00406 if(computeEigenvectors) 00407 m_eivec.setOnes(n,n); 00408 m_info = Success; 00409 m_isInitialized = true; 00410 m_eigenvectorsOk = computeEigenvectors; 00411 return *this; 00412 } 00413 00414 // declare some aliases 00415 RealVectorType& diag = m_eivalues; 00416 MatrixType& mat = m_eivec; 00417 00418 // map the matrix coefficients to [-1:1] to avoid over- and underflow. 00419 mat = matrix.template triangularView<Lower>(); 00420 RealScalar scale = mat.cwiseAbs().maxCoeff(); 00421 if(scale==RealScalar(0)) scale = RealScalar(1); 00422 mat.template triangularView<Lower>() /= scale; 00423 m_subdiag.resize(n-1); 00424 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); 00425 00426 Index end = n-1; 00427 Index start = 0; 00428 Index iter = 0; // total number of iterations 00429 00430 while (end>0) 00431 { 00432 for (Index i = start; i<end; ++i) 00433 if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) 00434 m_subdiag[i] = 0; 00435 00436 // find the largest unreduced block 00437 while (end>0 && m_subdiag[end-1]==0) 00438 { 00439 end--; 00440 } 00441 if (end<=0) 00442 break; 00443 00444 // if we spent too many iterations, we give up 00445 iter++; 00446 if(iter > m_maxIterations * n) break; 00447 00448 start = end - 1; 00449 while (start>0 && m_subdiag[start-1]!=0) 00450 start--; 00451 00452 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); 00453 } 00454 00455 if (iter <= m_maxIterations * n) 00456 m_info = Success; 00457 else 00458 m_info = NoConvergence; 00459 00460 // Sort eigenvalues and corresponding vectors. 00461 // TODO make the sort optional ? 00462 // TODO use a better sort algorithm !! 00463 if (m_info == Success) 00464 { 00465 for (Index i = 0; i < n-1; ++i) 00466 { 00467 Index k; 00468 m_eivalues.segment(i,n-i).minCoeff(&k); 00469 if (k > 0) 00470 { 00471 std::swap(m_eivalues[i], m_eivalues[k+i]); 00472 if(computeEigenvectors) 00473 m_eivec.col(i).swap(m_eivec.col(k+i)); 00474 } 00475 } 00476 } 00477 00478 // scale back the eigen values 00479 m_eivalues *= scale; 00480 00481 m_isInitialized = true; 00482 m_eigenvectorsOk = computeEigenvectors; 00483 return *this; 00484 } 00485 00486 00487 namespace internal { 00488 00489 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues 00490 { 00491 static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) 00492 { eig.compute(A,options); } 00493 }; 00494 00495 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false> 00496 { 00497 typedef typename SolverType::MatrixType MatrixType; 00498 typedef typename SolverType::RealVectorType VectorType; 00499 typedef typename SolverType::Scalar Scalar; 00500 typedef typename MatrixType::Index Index; 00501 00506 static inline void computeRoots(const MatrixType& m, VectorType& roots) 00507 { 00508 using std::sqrt; 00509 using std::atan2; 00510 using std::cos; 00511 using std::sin; 00512 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); 00513 const Scalar s_sqrt3 = sqrt(Scalar(3.0)); 00514 00515 // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The 00516 // eigenvalues are the roots to this equation, all guaranteed to be 00517 // real-valued, because the matrix is symmetric. 00518 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); 00519 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); 00520 Scalar c2 = m(0,0) + m(1,1) + m(2,2); 00521 00522 // Construct the parameters used in classifying the roots of the equation 00523 // and in solving the equation for the roots in closed form. 00524 Scalar c2_over_3 = c2*s_inv3; 00525 Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3; 00526 if(a_over_3<Scalar(0)) 00527 a_over_3 = Scalar(0); 00528 00529 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); 00530 00531 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b; 00532 if(q<Scalar(0)) 00533 q = Scalar(0); 00534 00535 // Compute the eigenvalues by solving for the roots of the polynomial. 00536 Scalar rho = sqrt(a_over_3); 00537 Scalar theta = atan2(sqrt(q),half_b)*s_inv3; // since sqrt(q) > 0, atan2 is in [0, pi] and theta is in [0, pi/3] 00538 Scalar cos_theta = cos(theta); 00539 Scalar sin_theta = sin(theta); 00540 // roots are already sorted, since cos is monotonically decreasing on [0, pi] 00541 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3) 00542 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3) 00543 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta; 00544 } 00545 00546 static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative) 00547 { 00548 using std::abs; 00549 Index i0; 00550 // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal): 00551 mat.diagonal().cwiseAbs().maxCoeff(&i0); 00552 // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector, 00553 // so let's save it: 00554 representative = mat.col(i0); 00555 Scalar n0, n1; 00556 VectorType c0, c1; 00557 n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm(); 00558 n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm(); 00559 if(n0>n1) res = c0/std::sqrt(n0); 00560 else res = c1/std::sqrt(n1); 00561 00562 return true; 00563 } 00564 00565 static inline void run(SolverType& solver, const MatrixType& mat, int options) 00566 { 00567 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); 00568 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00569 && (options&EigVecMask)!=EigVecMask 00570 && "invalid option parameter"); 00571 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 00572 00573 MatrixType& eivecs = solver.m_eivec; 00574 VectorType& eivals = solver.m_eivalues; 00575 00576 // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. 00577 Scalar shift = mat.trace() / Scalar(3); 00578 // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later 00579 MatrixType scaledMat = mat.template selfadjointView<Lower>(); 00580 scaledMat.diagonal().array() -= shift; 00581 Scalar scale = scaledMat.cwiseAbs().maxCoeff(); 00582 if(scale > 0) scaledMat /= scale; // TODO for scale==0 we could save the remaining operations 00583 00584 // compute the eigenvalues 00585 computeRoots(scaledMat,eivals); 00586 00587 // compute the eigenvectors 00588 if(computeEigenvectors) 00589 { 00590 if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon()) 00591 { 00592 // All three eigenvalues are numerically the same 00593 eivecs.setIdentity(); 00594 } 00595 else 00596 { 00597 MatrixType tmp; 00598 tmp = scaledMat; 00599 00600 // Compute the eigenvector of the most distinct eigenvalue 00601 Scalar d0 = eivals(2) - eivals(1); 00602 Scalar d1 = eivals(1) - eivals(0); 00603 Index k(0), l(2); 00604 if(d0 > d1) 00605 { 00606 std::swap(k,l); 00607 d0 = d1; 00608 } 00609 00610 // Compute the eigenvector of index k 00611 { 00612 tmp.diagonal().array () -= eivals(k); 00613 // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector. 00614 extract_kernel(tmp, eivecs.col(k), eivecs.col(l)); 00615 } 00616 00617 // Compute eigenvector of index l 00618 if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1) 00619 { 00620 // If d0 is too small, then the two other eigenvalues are numerically the same, 00621 // and thus we only have to ortho-normalize the near orthogonal vector we saved above. 00622 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l); 00623 eivecs.col(l).normalize(); 00624 } 00625 else 00626 { 00627 tmp = scaledMat; 00628 tmp.diagonal().array () -= eivals(l); 00629 00630 VectorType dummy; 00631 extract_kernel(tmp, eivecs.col(l), dummy); 00632 } 00633 00634 // Compute last eigenvector from the other two 00635 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized(); 00636 } 00637 } 00638 00639 // Rescale back to the original size. 00640 eivals *= scale; 00641 eivals.array() += shift; 00642 00643 solver.m_info = Success; 00644 solver.m_isInitialized = true; 00645 solver.m_eigenvectorsOk = computeEigenvectors; 00646 } 00647 }; 00648 00649 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel 00650 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false> 00651 { 00652 typedef typename SolverType::MatrixType MatrixType; 00653 typedef typename SolverType::RealVectorType VectorType; 00654 typedef typename SolverType::Scalar Scalar; 00655 00656 static inline void computeRoots(const MatrixType& m, VectorType& roots) 00657 { 00658 using std::sqrt; 00659 const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0))); 00660 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1)); 00661 roots(0) = t1 - t0; 00662 roots(1) = t1 + t0; 00663 } 00664 00665 static inline void run(SolverType& solver, const MatrixType& mat, int options) 00666 { 00667 using std::sqrt; 00668 using std::abs; 00669 00670 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); 00671 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00672 && (options&EigVecMask)!=EigVecMask 00673 && "invalid option parameter"); 00674 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 00675 00676 MatrixType& eivecs = solver.m_eivec; 00677 VectorType& eivals = solver.m_eivalues; 00678 00679 // map the matrix coefficients to [-1:1] to avoid over- and underflow. 00680 Scalar scale = mat.cwiseAbs().maxCoeff(); 00681 scale = (std::max)(scale,Scalar(1)); 00682 MatrixType scaledMat = mat / scale; 00683 00684 // Compute the eigenvalues 00685 computeRoots(scaledMat,eivals); 00686 00687 // compute the eigen vectors 00688 if(computeEigenvectors) 00689 { 00690 if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon()) 00691 { 00692 eivecs.setIdentity(); 00693 } 00694 else 00695 { 00696 scaledMat.diagonal().array () -= eivals(1); 00697 Scalar a2 = numext::abs2(scaledMat(0,0)); 00698 Scalar c2 = numext::abs2(scaledMat(1,1)); 00699 Scalar b2 = numext::abs2(scaledMat(1,0)); 00700 if(a2>c2) 00701 { 00702 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); 00703 eivecs.col(1) /= sqrt(a2+b2); 00704 } 00705 else 00706 { 00707 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); 00708 eivecs.col(1) /= sqrt(c2+b2); 00709 } 00710 00711 eivecs.col(0) << eivecs.col(1).unitOrthogonal(); 00712 } 00713 } 00714 00715 // Rescale back to the original size. 00716 eivals *= scale; 00717 00718 solver.m_info = Success; 00719 solver.m_isInitialized = true; 00720 solver.m_eigenvectorsOk = computeEigenvectors; 00721 } 00722 }; 00723 00724 } 00725 00726 template<typename MatrixType> 00727 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> 00728 ::computeDirect(const MatrixType& matrix, int options) 00729 { 00730 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options); 00731 return *this; 00732 } 00733 00734 namespace internal { 00735 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> 00736 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) 00737 { 00738 using std::abs; 00739 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); 00740 RealScalar e = subdiag[end-1]; 00741 // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still 00742 // underflow thus leading to inf/NaN values when using the following commented code: 00743 // RealScalar e2 = numext::abs2(subdiag[end-1]); 00744 // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); 00745 // This explain the following, somewhat more complicated, version: 00746 RealScalar mu = diag[end]; 00747 if(td==0) 00748 mu -= abs(e); 00749 else 00750 { 00751 RealScalar e2 = numext::abs2(subdiag[end-1]); 00752 RealScalar h = numext::hypot(td,e); 00753 if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); 00754 else mu -= e2 / (td + (td>0 ? h : -h)); 00755 } 00756 00757 RealScalar x = diag[start] - mu; 00758 RealScalar z = subdiag[start]; 00759 for (Index k = start; k < end; ++k) 00760 { 00761 JacobiRotation<RealScalar> rot; 00762 rot.makeGivens(x, z); 00763 00764 // do T = G' T G 00765 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; 00766 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; 00767 00768 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); 00769 diag[k+1] = rot.s() * sdk + rot.c() * dkp1; 00770 subdiag[k] = rot.c() * sdk - rot.s() * dkp1; 00771 00772 00773 if (k > start) 00774 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; 00775 00776 x = subdiag[k]; 00777 00778 if (k < end - 1) 00779 { 00780 z = -rot.s() * subdiag[k+1]; 00781 subdiag[k + 1] = rot.c() * subdiag[k+1]; 00782 } 00783 00784 // apply the givens rotation to the unit matrix Q = Q * G 00785 if (matrixQ) 00786 { 00787 // FIXME if StorageOrder == RowMajor this operation is not very efficient 00788 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n); 00789 q.applyOnTheRight(k,k+1,rot); 00790 } 00791 } 00792 } 00793 00794 } // end namespace internal 00795 00796 } // end namespace Eigen 00797 00798 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H