Tapkee
generalized_eigendecomposition.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2012-2013 Sergey Lisitsyn
00004  */
00005 
00006 #ifndef TAPKEE_GENERALIZED_EIGENDECOMPOSITION_H_
00007 #define TAPKEE_GENERALIZED_EIGENDECOMPOSITION_H_
00008 
00009 /* Tapkee includes */
00010 #ifdef TAPKEE_WITH_ARPACK
00011     #include <tapkee/utils/arpack_wrapper.hpp>
00012 #endif
00013 #include <tapkee/routines/matrix_operations.hpp>
00014 /* End of Tapkee includes */
00015 
00016 namespace tapkee
00017 {
00018 namespace tapkee_internal
00019 {
00020 
00021 #ifdef TAPKEE_WITH_ARPACK
00022 
00023 template <class LMatrixType, class RMatrixType, class MatrixOperationType> 
00024 EigendecompositionResult generalized_eigendecomposition_impl_arpack(const LMatrixType& lhs, 
00025         const RMatrixType& rhs, IndexType target_dimension, unsigned int skip)
00026 {
00027     timed_context context("ARPACK DSXUPD generalized eigendecomposition");
00028 
00029     ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperationType> 
00030         arpack(lhs,rhs,target_dimension+skip,"SM");
00031     
00032     if (arpack.info() == Eigen::Success)
00033     {
00034         std::stringstream ss;
00035         ss << "Took " << arpack.getNbrIterations() << " iterations.";
00036         LoggingSingleton::instance().message_info(ss.str());
00037         DenseMatrix selected_eigenvectors = (arpack.eigenvectors()).rightCols(target_dimension);
00038         return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
00039     }
00040     else
00041     {
00042         throw eigendecomposition_error("eigendecomposition failed");
00043     }
00044     return EigendecompositionResult();
00045 }
00046 #endif
00047 
00049 template <class LMatrixType, class RMatrixType, class MatrixOperationType> 
00050 EigendecompositionResult generalized_eigendecomposition_impl_dense(const LMatrixType& lhs,
00051         const RMatrixType& rhs, IndexType target_dimension, unsigned int skip)
00052 {
00053     timed_context context("Eigen dense generalized eigendecomposition");
00054 
00055     DenseMatrix dense_lhs = lhs;
00056     DenseMatrix dense_rhs = rhs;
00057     Eigen::GeneralizedSelfAdjointEigenSolver<DenseMatrix> solver(dense_lhs, dense_rhs);
00058     if (solver.info() == Eigen::Success)
00059     {
00060         if (MatrixOperationType::largest)
00061         {
00062             assert(skip==0);
00063             DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
00064             return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
00065         } 
00066         else
00067         {
00068             DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
00069             return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
00070         }
00071     }
00072     else
00073     {
00074         throw eigendecomposition_error("eigendecomposition failed");
00075     }
00076 
00077     return EigendecompositionResult();
00078 }
00079 
00080 template <class LMatrixType, class RMatrixType, class MatrixOperationType>
00081 EigendecompositionResult generalized_eigendecomposition(EigenMethod method, const LMatrixType& lhs,
00082                                                         const RMatrixType& rhs,
00083                                                         IndexType target_dimension, unsigned int skip)
00084 {
00085     LoggingSingleton::instance().message_info("Using the " + get_eigen_method_name(method) + " eigendecomposition method.");
00086     switch (method)
00087     {
00088 #ifdef TAPKEE_WITH_ARPACK
00089         case Arpack: 
00090             return generalized_eigendecomposition_impl_arpack<LMatrixType, RMatrixType, MatrixOperationType>(lhs, rhs, target_dimension, skip);
00091 #endif
00092         case Dense:
00093             return generalized_eigendecomposition_impl_dense<LMatrixType, RMatrixType, MatrixOperationType>(lhs, rhs, target_dimension, skip);
00094         case Randomized:
00095             throw unsupported_method_error("Randomized method is not supported for generalized eigenproblems");
00096             return EigendecompositionResult();
00097         default: break;
00098     }
00099     return EigendecompositionResult();
00100 }
00101 
00102 } // End of namespace tapkee_internal
00103 } // End of namespace tapkee
00104 
00105 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines