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