perturbation.hpp

// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#ifndef GKO_PUBLIC_CORE_BASE_PERTURBATION_HPP_
#define GKO_PUBLIC_CORE_BASE_PERTURBATION_HPP_
 
 
#include <memory>
 
 
#include <ginkgo/core/base/lin_op.hpp>
#include <ginkgo/core/matrix/dense.hpp>
 
 
namespace gko {
 
 
template <typename ValueType = default_precision>
class Perturbation : public EnableLinOp<Perturbation<ValueType>>,
                     public EnableCreateMethod<Perturbation<ValueType>> {
    friend class EnablePolymorphicObject<Perturbation, LinOp>;
    friend class EnableCreateMethod<Perturbation>;
 
public:
    using value_type = ValueType;
 
    const std::shared_ptr<const LinOp> get_basis() const noexcept
    {
        return basis_;
    }
 
    const std::shared_ptr<const LinOp> get_projector() const noexcept
    {
        return projector_;
    }
 
    const std::shared_ptr<const LinOp> get_scalar() const noexcept
    {
        return scalar_;
    }
 
    Perturbation& operator=(const Perturbation& other);
 
    Perturbation& operator=(Perturbation&& other);
 
    Perturbation(const Perturbation& other);
 
    Perturbation(Perturbation&& other);
 
    static std::unique_ptr<Perturbation> create(
        std::shared_ptr<const Executor> exec);
 
    static std::unique_ptr<Perturbation> create(
        std::shared_ptr<const LinOp> scalar,
        std::shared_ptr<const LinOp> basis);
 
    static std::unique_ptr<Perturbation> create(
        std::shared_ptr<const LinOp> scalar, std::shared_ptr<const LinOp> basis,
        std::shared_ptr<const LinOp> projector);
 
protected:
    explicit Perturbation(std::shared_ptr<const Executor> exec);
 
    explicit Perturbation(std::shared_ptr<const LinOp> scalar,
                          std::shared_ptr<const LinOp> basis);
 
    explicit Perturbation(std::shared_ptr<const LinOp> scalar,
                          std::shared_ptr<const LinOp> basis,
                          std::shared_ptr<const LinOp> projector);
 
    void apply_impl(const LinOp* b, LinOp* x) const override;
 
    void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta,
                    LinOp* x) const override;
 
    void validate_perturbation();
 
private:
    std::shared_ptr<const LinOp> basis_;
    std::shared_ptr<const LinOp> projector_;
    std::shared_ptr<const LinOp> scalar_;
 
    // TODO: solve race conditions when multithreading
    mutable struct cache_struct {
        cache_struct() = default;
        ~cache_struct() = default;
        cache_struct(const cache_struct& other) {}
        cache_struct& operator=(const cache_struct& other) { return *this; }
 
        // allocate linops of cache. The dimension of `intermediate` is
        // (the number of rows of projector, the number of columns of b). Others
        // are 1x1 scalar.
        void allocate(std::shared_ptr<const Executor> exec, dim<2> size)
        {
            using vec = gko::matrix::Dense<ValueType>;
            if (one == nullptr) {
                one = initialize<vec>({gko::one<ValueType>()}, exec);
            }
            if (alpha_scalar == nullptr) {
                alpha_scalar = vec::create(exec, gko::dim<2>(1));
            }
            if (intermediate == nullptr || intermediate->get_size() != size) {
                intermediate = vec::create(exec, size);
            }
        }
 
        std::unique_ptr<LinOp> intermediate;
        std::unique_ptr<LinOp> one;
        std::unique_ptr<LinOp> alpha_scalar;
    } cache_;
};
 
 
}  // namespace gko
 
 
#endif  // GKO_PUBLIC_CORE_BASE_PERTURBATION_HPP_