centered.hpp

Go to the documentation of this file.

#ifndef IRLBA_MATRIX_CENTERED_HPP
#define IRLBA_MATRIX_CENTERED_HPP
 
#include <memory>
 
#include "../utils.hpp"
#include "interface.hpp"
 
#include "Eigen/Dense"
 
namespace irlba {
 
template<class EigenVector_, class Matrix_, class Center_>
class CenteredWorkspace final : public Workspace<EigenVector_> {
public:
    CenteredWorkspace(const Matrix_& matrix, const Center_& center) : 
        my_work(matrix.new_known_workspace()),
        my_center(center)
    {}
private:
    I<decltype(std::declval<Matrix_>().new_known_workspace())> my_work;
    const Center_& my_center;
 
public:
    void multiply(const EigenVector_& right, EigenVector_& out) {
        my_work->multiply(right, out);
        const auto beta = right.dot(my_center);
        for (auto& o : out) {
            o -= beta;
        }
    }
};
 
template<class EigenVector_, class Matrix_, class Center_>
class CenteredAdjointWorkspace final : public AdjointWorkspace<EigenVector_> {
public:
    CenteredAdjointWorkspace(const Matrix_& matrix, const Center_& center) :
        my_work(matrix.new_known_adjoint_workspace()),
        my_center(center)
    {}
private:
    I<decltype(std::declval<Matrix_>().new_known_adjoint_workspace())> my_work;
    const Center_& my_center;
 
public:
    void multiply(const EigenVector_& right, EigenVector_& out) {
        my_work->multiply(right, out);
        const auto beta = right.sum();
        out -= beta * my_center;
    }
};
 
template<class EigenMatrix_, class Matrix_, class Center_>
class CenteredRealizeWorkspace final : public RealizeWorkspace<EigenMatrix_> {
public:
    CenteredRealizeWorkspace(const Matrix_& matrix, const Center_& center) :
        my_work(matrix.new_known_realize_workspace()),
        my_center(center)
    {}
private:
    I<decltype(std::declval<Matrix_>().new_known_realize_workspace())> my_work;
    const Center_& my_center;
 
public:
    const EigenMatrix_& realize(EigenMatrix_& buffer) {
        my_work->realize_copy(buffer);
        buffer.array().rowwise() -= my_center.adjoint().array();
        return buffer;
    }
};
 
template<
    class EigenVector_,
    class EigenMatrix_,
    class MatrixPointer_ = const Matrix<EigenVector_, EigenMatrix_>*,
    class CenterPointer_ = const EigenVector_* 
>
class CenteredMatrix final : public Matrix<EigenVector_, EigenMatrix_> {
public:
    CenteredMatrix(MatrixPointer_ matrix, CenterPointer_ center) :
        my_matrix(std::move(matrix)),
        my_center(std::move(center))
    {}
 
private:
    MatrixPointer_ my_matrix;
    CenterPointer_ my_center;
 
public:
    Eigen::Index rows() const {
        return my_matrix->rows();
    }
 
    Eigen::Index cols() const {
        return my_matrix->cols();
    }
 
public:
    std::unique_ptr<Workspace<EigenVector_> > new_workspace() const {
        return new_known_workspace();
    }
 
    std::unique_ptr<AdjointWorkspace<EigenVector_> > new_adjoint_workspace() const {
        return new_known_adjoint_workspace();
    }
 
    std::unique_ptr<RealizeWorkspace<EigenMatrix_> > new_realize_workspace() const {
        return new_known_realize_workspace();
    }
 
public:
    std::unique_ptr<CenteredWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > > new_known_workspace() const {
        return std::make_unique<CenteredWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > >(*my_matrix, *my_center);
    }
 
    std::unique_ptr<CenteredAdjointWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > > new_known_adjoint_workspace() const {
        return std::make_unique<CenteredAdjointWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > >(*my_matrix, *my_center);
    }
 
    std::unique_ptr<CenteredRealizeWorkspace<EigenMatrix_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > > new_known_realize_workspace() const {
        return std::make_unique<CenteredRealizeWorkspace<EigenMatrix_, I<decltype(*my_matrix)>, I<decltype(*my_center)> > >(*my_matrix, *my_center);
    }
};
 
}
 
#endif