power_iteration_eigenvalue_solver.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Pooyan Dadvand
//  Collaborator:    Vicente Mataix Ferrandiz
//
//
 
#if !defined(KRATOS_POWER_ITERATION_EIGENVALUE_SOLVER_H_INCLUDED )
#define  KRATOS_POWER_ITERATION_EIGENVALUE_SOLVER_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "spaces/ublas_space.h"
#include "includes/define.h"
#include "linear_solvers/iterative_solver.h"
#include "utilities/random_initializer_utility.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType, class TLinearSolverType,
         class TPreconditionerType = Preconditioner<TSparseSpaceType, TDenseSpaceType>,
         class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType> >
class PowerIterationEigenvalueSolver
    : public IterativeSolver<TSparseSpaceType, TDenseSpaceType, TPreconditionerType, TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(PowerIterationEigenvalueSolver);
 
    typedef IterativeSolver<TSparseSpaceType, TDenseSpaceType, TPreconditionerType, TReordererType> BaseType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
    typedef typename TDenseSpaceType::VectorType DenseVectorType;
 
    typedef std::size_t SizeType;
 
    typedef std::size_t IndexType;
 
 
    PowerIterationEigenvalueSolver() {}
 
    PowerIterationEigenvalueSolver(
        double MaxTolerance,
        unsigned int MaxIterationNumber,
        unsigned int RequiredEigenvalueNumber,
        typename TLinearSolverType::Pointer pLinearSolver
    ): BaseType(MaxTolerance, MaxIterationNumber),
       mRequiredEigenvalueNumber(RequiredEigenvalueNumber),
       mpLinearSolver(pLinearSolver)
    {
 
    }
 
    PowerIterationEigenvalueSolver(
        Parameters ThisParameters,
        typename TLinearSolverType::Pointer pLinearSolver
        ): mpLinearSolver(pLinearSolver)
    {
        Parameters DefaultParameters = Parameters(R"(
        {
            "solver_type"             : "power_iteration_eigenvalue_solver",
            "max_iteration"           : 10000,
            "tolerance"               : 1e-8,
            "required_eigen_number"   : 1,
            "shifting_convergence"    : 0.25,
            "verbosity"               : 1,
            "linear_solver_settings"  : {}
        })" );
 
        ThisParameters.ValidateAndAssignDefaults(DefaultParameters);
 
        mRequiredEigenvalueNumber = ThisParameters["required_eigen_number"].GetInt();
        mEchoLevel = ThisParameters["verbosity"].GetInt();
        BaseType::SetTolerance( ThisParameters["tolerance"].GetDouble() );
        BaseType::SetMaxIterationsNumber( ThisParameters["max_iteration"].GetInt() );
    }
 
    PowerIterationEigenvalueSolver(const PowerIterationEigenvalueSolver& Other) : BaseType(Other)
    {
 
    }
 
 
    ~PowerIterationEigenvalueSolver() override
    {
 
    }
 
 
 
    PowerIterationEigenvalueSolver& operator=(const PowerIterationEigenvalueSolver& Other)
    {
        BaseType::operator=(Other);
        return *this;
    }
 
 
    void Solve(
        SparseMatrixType& K,
        SparseMatrixType& M,
        DenseVectorType& Eigenvalues,
        DenseMatrixType& Eigenvectors
        ) override
    {
 
        const SizeType size = K.size1();
        const SizeType max_iteration = BaseType::GetMaxIterationsNumber();
        const double tolerance = BaseType::GetTolerance();
 
        VectorType x = boost::numeric::ublas::zero_vector<double>(size);
        VectorType y = boost::numeric::ublas::zero_vector<double>(size);
 
        RandomInitializeUtility<double>::RandomInitialize(K, y);
 
        if(Eigenvalues.size() < 1)
            Eigenvalues.resize(1, 0.0);
 
        // Starting with first step
        double beta = 0.0;
        double rho = 0.0;
        double old_rho = Eigenvalues[0];
        VectorType y_old = boost::numeric::ublas::zero_vector<double>(size);
 
        for(SizeType i = 0 ; i < max_iteration ; i++) {
            // K*x = y
            mpLinearSolver->Solve(K, x, y);
 
            rho = inner_prod(y, x);
 
            // y = M*x
            TSparseSpaceType::Mult(M, x, y);
            beta = inner_prod(x, y);
 
            KRATOS_ERROR_IF(beta <= 0.0) << "M is not Positive-definite. beta = " << beta << std::endl;
 
            rho /= beta;
            beta = std::sqrt(beta);
            TSparseSpaceType::InplaceMult(y, 1.0/beta);
 
            KRATOS_ERROR_IF(rho == 0.0) << "Perpendicular eigenvector to M" << std::endl;
 
            const double convergence_rho = std::abs((rho - old_rho) / rho);
            const double convergence_norm = TSparseSpaceType::TwoNorm(y - y_old)/TSparseSpaceType::TwoNorm(y);
 
            if (mEchoLevel > 1)
                KRATOS_INFO("Power Iterator Eigenvalue Solver: ") << "Iteration: " << i << " \t beta: " << beta << "\trho: " << rho << " \tConvergence norm: " << convergence_norm << " \tConvergence rho: " << convergence_rho << std::endl;
 
            if(convergence_norm < tolerance || convergence_rho < tolerance)
                break;
 
            old_rho = rho;
            TSparseSpaceType::Assign(y_old, 1.0, y);
        }
 
        if (mEchoLevel > 0) {
            KRATOS_INFO("rho: ") << rho << std::endl;
            KRATOS_INFO("y: ") << y << std::endl;
        }
 
        Eigenvalues[0] = rho;
 
        if((Eigenvectors.size1() < 1) || (Eigenvectors.size2() < size))
            Eigenvectors.resize(1,size);
 
        for(SizeType i = 0 ; i < size ; i++)
            Eigenvectors(0,i) = y[i];
    }
 
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Power iteration eigenvalue solver with " << BaseType::GetPreconditioner()->Info();
        return  buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        BaseType::PrintData(rOStream);
    }
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
    unsigned int mRequiredEigenvalueNumber;             
 
    unsigned int mEchoLevel;                            
 
    typename TLinearSolverType::Pointer mpLinearSolver; 
 
 
 
 
 
 
 
 
 
 
 
 
}; // Class PowerIterationEigenvalueSolver
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::istream& operator >> (std::istream& IStream,
                                  PowerIterationEigenvalueSolver<TSparseSpaceType, TDenseSpaceType,
                                  TPreconditionerType, TReordererType>& rThis)
{
    return IStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::ostream& operator << (std::ostream& OStream,
                                  const PowerIterationEigenvalueSolver<TSparseSpaceType, TDenseSpaceType,
                                  TPreconditionerType, TReordererType>& rThis)
{
    rThis.PrintInfo(OStream);
    OStream << std::endl;
    rThis.PrintData(OStream);
 
    return OStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_POWER_ITERATION_EIGENVALUE_SOLVER_H_INCLUDED defined