nullspace.hpp

// SPDX-FileCopyrightText: 2026 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#ifndef GKO_PUBLIC_CORE_SOLVER_NULLSPACE_HPP_
#define GKO_PUBLIC_CORE_SOLVER_NULLSPACE_HPP_
 
 
#include <memory>
#include <vector>
 
#include <ginkgo/core/base/array.hpp>
#include <ginkgo/core/base/dense_cache.hpp>
#include <ginkgo/core/base/lin_op.hpp>
#include <ginkgo/core/base/math.hpp>
#include <ginkgo/core/base/polymorphic_object.hpp>
#include <ginkgo/core/base/utils.hpp>
#include <ginkgo/core/matrix/dense.hpp>
 
 
namespace gko {
namespace solver {
 
 
template <typename ValueType = default_precision>
class Nullspace : public LinOp, public EnableCloneable<Nullspace<ValueType>> {
    friend class EnableCloneable<Nullspace>;
    GKO_ASSERT_SUPPORTED_VALUE_TYPE;
 
public:
    using EnableCloneable<Nullspace>::convert_to;
    using EnableCloneable<Nullspace>::move_to;
 
    using value_type = ValueType;
    using absolute_type = remove_complex<ValueType>;
 
    static std::unique_ptr<Nullspace> create(
        std::shared_ptr<const Executor> exec,
        std::vector<std::shared_ptr<const LinOp>> basis,
        bool contains_constant = false);
 
    static std::unique_ptr<Nullspace> create_from_constant(
        std::shared_ptr<const Executor> exec, dim<2> size = {});
 
    bool contains_constant() const noexcept;
 
    size_type get_dimension() const noexcept;
 
    std::shared_ptr<const LinOp> get_basis() const noexcept;
 
    void project(ptr_param<LinOp> v) const;
 
    Nullspace& operator=(const Nullspace& other);
 
    Nullspace& operator=(Nullspace&& other);
 
    Nullspace(const Nullspace& other);
 
    Nullspace(Nullspace&& other);
 
protected:
    explicit Nullspace(std::shared_ptr<const Executor> exec);
 
    Nullspace(std::shared_ptr<const Executor> exec, dim<2> size,
              std::vector<std::shared_ptr<const LinOp>> basis,
              bool contains_constant);
 
    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;
 
private:
    template <typename VectorType>
    void setup_basis(const std::vector<std::shared_ptr<const LinOp>>& basis);
 
    template <typename VectorType>
    const matrix::Dense<ValueType>* get_local_basis(const VectorType* v) const;
 
    // stores [mean(v); V^H v] in coefficients_
    template <typename VectorType>
    void compute_components(const VectorType* v,
                            const matrix::Dense<ValueType>* basis_local,
                            bool has_constant) const;
 
    // result = v - mean - V C with the coefficients from compute_components;
    // result may be v
    template <typename VectorType>
    void subtract_components(const VectorType* v, VectorType* result,
                             const matrix::Dense<ValueType>* basis_local,
                             bool has_constant) const;
 
    template <typename VectorType>
    void remove_components(VectorType* v,
                           const matrix::Dense<ValueType>* basis_local,
                           bool has_constant) const;
 
    bool contains_constant_;
    // orthonormal explicit basis (n x k), Dense or distributed::Vector
    std::shared_ptr<const LinOp> basis_;
    // scratch: (1 + k) x nrhs coefficients (means, then V^H x)
    gko::detail::DenseCache<ValueType> coefficients_;
    gko::detail::DenseCache<ValueType> host_coefficients_;
    mutable array<char> reduction_tmp_;
};
 
 
namespace detail {
 
 
template <typename ValueType>
std::shared_ptr<const Nullspace<ValueType>> prepare_nullspace(
    std::shared_ptr<const LinOp> param, const char* name, const LinOp* op,
    size_type size, std::shared_ptr<const Executor> exec);
 
 
}  // namespace detail
}  // namespace solver
}  // namespace gko
 
 
#endif  // GKO_PUBLIC_CORE_SOLVER_NULLSPACE_HPP_