gsLinearOperator.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsIO/gsOptionList.h>

namespace gismo
{

template<class T>
class gsLinearOperator
{
public:

    typedef memory::shared_ptr<gsLinearOperator> Ptr;

    typedef memory::unique_ptr<gsLinearOperator> uPtr;

    static gsIdentityOp<T> Identity(const index_t dim) {return gsIdentityOp<T>(dim); }

    virtual ~gsLinearOperator() {}

    virtual void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const = 0;

    virtual index_t rows() const = 0;

    virtual index_t cols() const = 0;

    // NOTE: this is rather inefficient and is only provided for debugging and testing purposes
    void toMatrix(gsMatrix<T>& result)
    {
        gsMatrix<T> eye = gsMatrix<T>::Identity(cols(), cols());
        this->apply(eye, result);
    }

    // This implementation provides an empty object
    static gsOptionList defaultOptions()               { return gsOptionList(); }

    // This implementation does not read any input
    virtual void setOptions(const gsOptionList &)  {}

}; // gsLinearOperator

template<class T>
class gsScaledOp GISMO_FINAL : public gsLinearOperator<T>
{
public:
    typedef memory::shared_ptr<gsScaledOp> Ptr;

    typedef memory::unique_ptr<gsScaledOp> uPtr;

    typedef typename gsLinearOperator<T>::Ptr BasePtr;

    gsScaledOp(BasePtr op, T scalar = 1) : m_op(give(op)), m_scalar(scalar)    {}

    static uPtr make(BasePtr op, T scalar = 1)
    { return uPtr( new gsScaledOp(give(op), scalar) ); }

    virtual void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        m_op->apply(input, x);
        x *= m_scalar;
    }

    index_t rows() const {return m_op->rows();}

    index_t cols() const {return m_op->cols();}

private:
    const BasePtr m_op;
    const T m_scalar;
}; // gsScaladOp


template<class T>
class gsIdentityOp GISMO_FINAL : public gsLinearOperator<T>
{
public:

    typedef memory::shared_ptr<gsIdentityOp> Ptr;

    typedef memory::unique_ptr<gsIdentityOp> uPtr;

    gsIdentityOp(index_t dim) : m_dim(dim) {}

    static uPtr make(index_t dim) { return uPtr( new gsIdentityOp(dim) ); }

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        x = input;
    }

    index_t rows() const {return m_dim;}

    index_t cols() const {return m_dim;}

private:
    const index_t m_dim;
};

template<class T, class L>
class gsLinearLambdaOp GISMO_FINAL : public gsLinearOperator<T>
{
public:
    gsLinearLambdaOp(L lambda, index_t rows, index_t cols)
        : m_lambda(lambda), m_rows(rows), m_cols(cols) {}
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        m_lambda(input, x);
    }
    index_t rows() const {return m_rows;}
    index_t cols() const {return m_cols;}
private:
    L m_lambda;
    index_t m_rows, m_cols;
};

template<class T=real_t, class L>
typename gsLinearOperator<T>::uPtr makeLinearOp(L lambda, index_t rows, index_t cols)
{
    return typename gsLinearOperator<T>::uPtr(new gsLinearLambdaOp<T,L>(lambda,rows,cols));
}


} // namespace gismo