gsProductOp.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsLinearOperator.h>

namespace gismo {
template<typename T>
class gsProductOp GISMO_FINAL : public gsLinearOperator<T>
{
    typedef typename gsLinearOperator<T>::Ptr BasePtr;
public:

    typedef memory::shared_ptr<gsProductOp> Ptr;

    typedef memory::unique_ptr<gsProductOp> uPtr;

    gsProductOp() : m_ops() {}

    gsProductOp(std::vector<BasePtr> ops) : m_ops(give(ops))
    {
#ifndef NDEBUG
        const index_t sz = m_ops.size();
        for (index_t i=0; i<sz-1; ++i)
        {
            GISMO_ASSERT ( m_ops[i]->rows() == m_ops[i+1]->cols(),
                           "Dimensions of the operators do not fit." );
        }
#endif
    }

    gsProductOp(BasePtr op0, BasePtr op1) : m_ops(2)
    {
        GISMO_ASSERT ( op0->rows() == op1->cols(),
                       "Dimensions of the operators do not fit." );
        m_ops[0]=give(op0); m_ops[1]=give(op1);
    }

    gsProductOp(BasePtr op0, BasePtr op1, BasePtr op2) : m_ops(3)
    {
        GISMO_ASSERT ( op0->rows() == op1->cols() && op1->rows() == op2->cols(),
                       "Dimensions of the operators do not fit." );
        m_ops[0]=give(op0); m_ops[1]=give(op1); m_ops[2]=give(op2);
    }

    static uPtr make()
    { return uPtr( new gsProductOp() ); }

    static uPtr make(std::vector<BasePtr> ops)
    { return uPtr( new gsProductOp(give(ops)) ); }

    static uPtr make(BasePtr op0, BasePtr op1)
    { return uPtr( new gsProductOp(give(op0),give(op1)) ); }

    static uPtr make(BasePtr op0, BasePtr op1, BasePtr op2)
    { return uPtr( new gsProductOp(give(op0),give(op1),give(op2)) ); }

    void addOperator( BasePtr op )
    {
        GISMO_ASSERT ( m_ops.empty() || m_ops.back()->rows() == op->cols(),
                       "Dimensions of the operators do not fit." );
        m_ops.push_back(give(op));
    }

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        // The product of 0 operators is the identity
        if ( m_ops.empty() ) { x = input; return; }

        // Here, we could make a permanently allocated vector
        gsMatrix<T> temp;
        const size_t sz = m_ops.size();

        m_ops[0]->apply(input,x);
        for(size_t i=1; i<sz;++i)
        {
            temp.swap(x);
            m_ops[i]->apply(temp,x);
        }
    }

    index_t rows() const {
        GISMO_ASSERT( !m_ops.empty(), "gsProductOp::rows does not work for 0 operators.");
        return m_ops.back()->rows();

    }

    index_t cols() const {
        GISMO_ASSERT( !m_ops.empty(), "gsProductOp::cols does not work for 0 operators.");
        return m_ops.front()->cols();
    }

    const std::vector<BasePtr>& getOps() const { return m_ops; }

private:
    std::vector<BasePtr> m_ops;
};

} // namespace gismo