condition_number_utility.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#if !defined(KRATOS_COND_NUMBER_UTILITY )
#define  KRATOS_COND_NUMBER_UTILITY
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "spaces/ublas_space.h"
#include "utilities/math_utils.h"
 
// Linear solvers
#include "linear_solvers/linear_solver.h"
#include "linear_solvers/skyline_lu_factorization_solver.h"
#include "linear_solvers/power_iteration_eigenvalue_solver.h"
#include "linear_solvers/power_iteration_highest_eigenvalue_solver.h"
 
namespace Kratos
{
 
 
 
 
 
 
class ConditionNumberUtility
{
public:
 
 
    KRATOS_CLASS_POINTER_DEFINITION( ConditionNumberUtility );
 
    typedef std::size_t SizeType;
 
    typedef std::size_t IndexType;
 
    typedef UblasSpace<double, CompressedMatrix, boost::numeric::ublas::vector<double>> SparseSpaceType;
 
    typedef UblasSpace<double, Matrix, Vector> LocalSpaceType;
 
    typedef SparseSpaceType::MatrixType SparseMatrixType;
 
    typedef SparseSpaceType::VectorType VectorType;
 
    typedef LocalSpaceType::MatrixType DenseMatrixType;
 
    typedef LocalSpaceType::VectorType DenseVectorType;
 
    typedef LinearSolver<SparseSpaceType, LocalSpaceType> LinearSolverType;
 
    typedef Reorderer<SparseSpaceType,  LocalSpaceType > ReordererType;
 
    typedef SkylineLUFactorizationSolver<SparseSpaceType,  LocalSpaceType, ReordererType > SkylineLUFactorizationSolverType;
 
    typedef PowerIterationHighestEigenvalueSolver<SparseSpaceType,  LocalSpaceType, LinearSolverType > PowerIterationHighestEigenvalueSolverType;
 
    typedef PowerIterationEigenvalueSolver<SparseSpaceType,  LocalSpaceType, LinearSolverType > PowerIterationEigenvalueSolverType;
 
 
 
    ConditionNumberUtility()
    {
        // definition of empty parameters
        Parameters empty_parameters = Parameters(R"({})" );
        LinearSolverType::Pointer psolver = LinearSolverType::Pointer( new SkylineLUFactorizationSolverType() );
        mpEigenSolverMax = LinearSolverType::Pointer( new PowerIterationHighestEigenvalueSolverType(empty_parameters, psolver) );
        mpEigenSolverMin = LinearSolverType::Pointer( new PowerIterationEigenvalueSolverType(empty_parameters, psolver) );
    }
 
    ConditionNumberUtility(
        LinearSolverType::Pointer pEigenSolverMax,
        LinearSolverType::Pointer pEigenSolverMin
        ) :mpEigenSolverMax(pEigenSolverMax),
           mpEigenSolverMin(pEigenSolverMin)
    {}
 
    virtual ~ConditionNumberUtility(){}
 
 
    double GetConditionNumber(SparseMatrixType& rInputMatrix)
    {
        KRATOS_ERROR_IF(mpEigenSolverMax == nullptr || mpEigenSolverMin == nullptr) << "ERROR:: PLEASE DEFINE THE EigenSolvers" << std::endl;
        return GetConditionNumber(rInputMatrix, mpEigenSolverMax, mpEigenSolverMin);
    }
 
    double GetConditionNumber(
        SparseMatrixType& rInputMatrix,
        LinearSolverType::Pointer pEigenSolverMax,
        LinearSolverType::Pointer pEigenSolverMin
        )
    {
        // The eigen system
        DenseVectorType eigen_values;
        DenseMatrixType eigen_vectors;
 
        const SizeType size_matrix = rInputMatrix.size1();
 
        SparseMatrixType identity_matrix(size_matrix, size_matrix);
        for (IndexType i = 0; i < size_matrix; ++i)
            identity_matrix.push_back(i, i, 1.0);
 
        pEigenSolverMax->Solve(rInputMatrix, identity_matrix, eigen_values, eigen_vectors);
        const double max_lambda = eigen_values[0];
 
        pEigenSolverMin->Solve(rInputMatrix, identity_matrix, eigen_values, eigen_vectors);
        const double min_lambda = eigen_values[0];
 
        KRATOS_ERROR_IF(min_lambda < std::numeric_limits<double>::epsilon()) << "ERROR:: NOT POSSIBLE TO COMPUTE CONDITION NUMBER. ZERO EIGENVALUE" << std::endl;
 
        const double condition_number = std::abs(max_lambda)/std::abs(min_lambda);
 
        return condition_number;
    }
 
 
 
 
 
 
 
private:
 
 
 
    LinearSolverType::Pointer mpEigenSolverMax; 
    LinearSolverType::Pointer mpEigenSolverMin; 
 
 
 
 
 
 
 
}; /* Class ConditionNumberUtility */
 
 
 
 
}  /* namespace Kratos.*/
 
#endif /* KRATOS_COND_NUMBER_UTILITY  defined */