SearchSpace.h

Go to the documentation of this file.

// Copyright (C) 2020 Netherlands eScience Center <n.renauld@esciencecenter.nl>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#ifndef SPECTRA_SEARCH_SPACE_H
#define SPECTRA_SEARCH_SPACE_H
 
#include <Eigen/Core>
 
#include "RitzPairs.h"
#include "Orthogonalization.h"
 
namespace Spectra {
 
template <typename Scalar>
class SearchSpace
{
private:
    using Index = Eigen::Index;
    using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
 
    Matrix m_basis_vectors;
    Matrix m_op_basis_product;
 
    void append_new_vectors_to_basis(const Matrix& new_vect)
    {
        Index num_update = new_vect.cols();
        m_basis_vectors.conservativeResize(Eigen::NoChange, m_basis_vectors.cols() + num_update);
        m_basis_vectors.rightCols(num_update).noalias() = new_vect;
    }
 
public:
    SearchSpace() = default;
 
    Index size() const { return m_basis_vectors.cols(); }
 
    void initialize_search_space(const Eigen::Ref<const Matrix>& initial_vectors)
    {
        m_basis_vectors = initial_vectors;
        m_op_basis_product = Matrix(initial_vectors.rows(), 0);
    }
 
    template <typename OpType>
    void update_operator_basis_product(OpType& op)
    {
        Index nvec = m_basis_vectors.cols() - m_op_basis_product.cols();
        m_op_basis_product.conservativeResize(Eigen::NoChange, m_basis_vectors.cols());
        m_op_basis_product.rightCols(nvec).noalias() = op * m_basis_vectors.rightCols(nvec);
    }
 
    void restart(const RitzPairs<Scalar>& ritz_pairs, Index size)
    {
        m_basis_vectors = ritz_pairs.ritz_vectors().leftCols(size);
        m_op_basis_product = m_op_basis_product * ritz_pairs.small_ritz_vectors().leftCols(size);
    }
 
    void extend_basis(const Matrix& new_vect)
    {
        Index left_cols_to_skip = size();
        append_new_vectors_to_basis(new_vect);
        twice_is_enough_orthogonalisation(m_basis_vectors, left_cols_to_skip);
    }
 
    const Matrix& basis_vectors() const { return m_basis_vectors; }
 
    const Matrix& operator_basis_product() const { return m_op_basis_product; }
};
 
}  // namespace Spectra
 
#endif  // SPECTRA_SEARCH_SPACE_H