scaled.hpp

Go to the documentation of this file.

#ifndef IRLBA_MATRIX_SCALED_HPP
#define IRLBA_MATRIX_SCALED_HPP
 
#include <memory>
 
#include "../utils.hpp"
#include "interface.hpp"
 
#include "Eigen/Dense"
 
namespace irlba {
 
template<class EigenVector_, class Matrix_, class Scale_>
class ScaledWorkspace final : public Workspace<EigenVector_> {
public:
    ScaledWorkspace(const Matrix_& matrix, const Scale_& scale, const bool column, const bool divide) : 
        my_work(matrix.new_known_workspace()),
        my_scale(scale),
        my_column(column),
        my_divide(divide)
    {}
 
private:
    I<decltype(std::declval<Matrix_>().new_known_workspace())> my_work;
    const Scale_& my_scale;
    bool my_column;
    bool my_divide;
 
    EigenVector_ my_buffer;
 
public:
    void multiply(const EigenVector_& right, EigenVector_& out) {
        if (my_column) {
            if (my_divide) {
                my_buffer = right.cwiseQuotient(my_scale);
            } else {
                my_buffer = right.cwiseProduct(my_scale);
            }
            my_work->multiply(my_buffer, out);
 
        } else {
            my_work->multiply(right, out);
            if (my_divide) {
                out.array() /= my_scale.array();
            } else {
                out.array() *= my_scale.array();
            }
        }
    }
};
 
template<class EigenVector_, class Matrix_, class Scale_>
class ScaledAdjointWorkspace final : public AdjointWorkspace<EigenVector_> {
public:
    ScaledAdjointWorkspace(const Matrix_& matrix, const Scale_& scale, const bool column, const bool divide) :
        my_work(matrix.new_known_adjoint_workspace()),
        my_scale(scale),
        my_column(column),
        my_divide(divide)
    {}
 
private:
    I<decltype(std::declval<Matrix_>().new_known_adjoint_workspace())> my_work;
    const Scale_& my_scale;
    bool my_column;
    bool my_divide;
 
    EigenVector_ my_buffer;
 
public:
    void multiply(const EigenVector_& right, EigenVector_& out) {
        if (my_column) {
            my_work->multiply(right, out);
            if (my_divide) {
                out.array() /= my_scale.array();
            } else {
                out.array() *= my_scale.array();
            }
 
        } else {
            if (my_divide) {
                my_buffer = right.cwiseQuotient(my_scale);
            } else {
                my_buffer = right.cwiseProduct(my_scale);
            }
            my_work->multiply(my_buffer, out);
        }
    }
};
 
template<class EigenMatrix_, class Matrix_, class Scale_>
class ScaledRealizeWorkspace final : public RealizeWorkspace<EigenMatrix_> {
public:
    ScaledRealizeWorkspace(const Matrix_& matrix, const Scale_& scale, const bool column, const bool divide) :
        my_work(matrix.new_known_realize_workspace()),
        my_scale(scale),
        my_column(column),
        my_divide(divide)
    {}
 
private:
    I<decltype(std::declval<Matrix_>().new_known_realize_workspace())> my_work;
    const Scale_& my_scale;
    bool my_column;
    bool my_divide;
 
public:
    const EigenMatrix_& realize(EigenMatrix_& buffer) {
        my_work->realize_copy(buffer);
 
        if (my_column) {
            if (my_divide) {
                buffer.array().rowwise() /= my_scale.adjoint().array();
            } else {
                buffer.array().rowwise() *= my_scale.adjoint().array();
            }
 
        } else {
            if (my_divide) {
                buffer.array().colwise() /= my_scale.array();
            } else {
                buffer.array().colwise() *= my_scale.array();
            }
        }
 
        return buffer;
    }
};
template<
    class EigenVector_,
    class EigenMatrix_,
    class MatrixPointer_ = const Matrix<EigenVector_, EigenMatrix_>*,
    class ScalePointer_ = const EigenVector_*
>
class ScaledMatrix final : public Matrix<EigenVector_, EigenMatrix_> {
public:
    ScaledMatrix(MatrixPointer_ matrix, ScalePointer_ scale, bool column, bool divide) : 
        my_matrix(std::move(matrix)),
        my_scale(std::move(scale)),
        my_column(column),
        my_divide(divide)
    {}
 
private:
    MatrixPointer_ my_matrix;
    ScalePointer_ my_scale;
    bool my_column;
    bool my_divide;
 
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:
    auto new_known_workspace() const {
        return std::make_unique<ScaledWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_scale)> > >(*my_matrix, *my_scale, my_column, my_divide);
    }
 
    auto new_known_adjoint_workspace() const {
        return std::make_unique<ScaledAdjointWorkspace<EigenVector_, I<decltype(*my_matrix)>, I<decltype(*my_scale)> > >(*my_matrix, *my_scale, my_column, my_divide);
    }
 
    auto new_known_realize_workspace() const {
        return std::make_unique<ScaledRealizeWorkspace<EigenMatrix_, I<decltype(*my_matrix)>, I<decltype(*my_scale)> > >(*my_matrix, *my_scale, my_column, my_divide);
    }
};
 
}
 
#endif