constrained_generalized_hermitian_eigen_solver.hpp

Go to the documentation of this file.

/*
* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

#ifndef __mast_slepc_constrained_generalized_hermitian_eigen_solver_h__
#define __mast_slepc_constrained_generalized_hermitian_eigen_solver_h__

// C++ includes
#include <iomanip>

// MAST includes
#include <mast/base/mast_data_types.h>
#include <mast/base/exceptions.hpp>

// Eigen includes
#include <Eigen/Eigenvalues>


namespace MAST {
namespace Solvers {
namespace EigenWrapper {

template <typename ScalarType,
          typename EigenSolverType,
          typename MatType,
          typename VectorType>
class ConstrainedGeneralizedHermitianEigenSolver {
    
public:

    using scalar_t = ScalarType;
    
    ConstrainedGeneralizedHermitianEigenSolver(const std::vector<uint_t> &dofs):
    _dofs         (dofs),
    _initialized  (false),
    _n            (0) { }
    
    virtual ~ConstrainedGeneralizedHermitianEigenSolver() { }
        
    inline ScalarType eig(uint_t i) {
        
        Assert0(_initialized, "solver not initialized");
        Assert2(i < _n,
                i, _n,
                "Eigenvalue index must be less than matrix size");
        
        return _solver.eigenvalues()(i);
    }
    
    
    inline void getEigenVector(uint_t        i,
                               VectorType   &x) {
        
        Assert0(_initialized, "solver not initialized");
        Assert2(i < _n,
                i, _n,
                "Eigenvalue index must be less than matrix size");

        Assert2(x.size() == _n,
                x.size(), _n,
                "Vector must have dimension of original matrix");
        

        x.setZero();
        
        for (uint_t j=0; j<_dofs.size(); j++) {
            x(_dofs[j]) = _solver.eigenvectors().col(i)(j);
        }
    }
    
    inline void solve(MatType          &A_mat,
                      MatType          &B_mat,
                      bool              computeEigenvectors) {
        
        Assert0(!_initialized, "solver already initialized");
        
        Assert2(A_mat.rows() == A_mat.cols(),
                A_mat.rows(), A_mat.cols(),
                "Matrix must be square");
        Assert2(B_mat.rows() == B_mat.cols(),
                B_mat.rows(), B_mat.cols(),
                "Matrix must be square");
        Assert2(A_mat.rows() == B_mat.rows(),
                A_mat.rows(), B_mat.rows(),
                "Matrices must have same dimensions");

        _n = A_mat.rows();
        
        // initialize the matrices
        _init_sub_matrices(A_mat, B_mat, _A_sub, _B_sub);
        
        // create the solver context
        if (computeEigenvectors)
            _solver.compute(_A_sub, _B_sub, Eigen::ComputeEigenvectors|Eigen::Ax_lBx);
        else
            _solver.compute(_A_sub, _B_sub, Eigen::EigenvaluesOnly);
        
        _initialized = true;
    }
    

    inline scalar_t sensitivity_solve(MatType      &B,
                                      MatType      &A_sens,
                                      MatType      &B_sens,
                                      uint_t       i) {

        VectorType
        v1 = VectorType::Zero(_n);
        
        this->getEigenVector(i, v1);
        
        scalar_t
        eig   = this->eig(i);

        eig =
        v1.dot( (A_sens * v1) - (eig * B_sens * v1) )/ // numerator
        v1.dot(B*v1); // denominator
        
        return eig;
    }
    

protected:

    inline void
    _init_sub_matrices(const Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic> &A,
                       const Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic> &B,
                       Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>       &A_sub,
                       Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>       &B_sub) {
        
        Assert0(!_initialized, "solver already initialized");
        
        A_sub.setZero(_dofs.size(), _dofs.size());
        B_sub.setZero(_dofs.size(), _dofs.size());
        
        for (uint_t i=0; i<_dofs.size(); i++) {
            for (uint_t j=0; j<_dofs.size(); j++) {
                A_sub(i, j) = A(_dofs[i], _dofs[j]);
                B_sub(i, j) = B(_dofs[i], _dofs[j]);
            }
        }
    }
    
    
    const std::vector<uint_t>    &_dofs;
    bool                          _initialized;
    int                           _n;
    MatType                       _A_sub;
    MatType                       _B_sub;
    EigenSolverType               _solver;
};

}  // namespace EigenWrapper
}  // namespace Solvers
}  // namespace MAST

#endif // __mast_slepc_constrained_generalized_hermitian_eigen_solver_h__