gsSumOp.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsLinearOperator.h>

namespace gismo
{

template<typename T>
class gsSumOp GISMO_FINAL : public gsLinearOperator<T>
{
    typedef typename gsLinearOperator<T>::Ptr BasePtr;
public:

    typedef memory::shared_ptr<gsSumOp> Ptr;

    typedef memory::unique_ptr<gsSumOp> uPtr;

    gsSumOp() : m_ops(0) {}

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

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

    gsSumOp(BasePtr op0, BasePtr op1, BasePtr op2 )
        : m_ops(3)
    {
        GISMO_ASSERT ( op0->rows() == op1->rows() && op0->cols() == op1->cols()
                        && op0->rows() == op2->rows() && op0->cols() == 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 gsSumOp() ); }

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

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

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

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

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        GISMO_ASSERT ( !m_ops.empty(), "gsSumOp::apply does not work for 0 operators." );

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

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

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

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

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

private:
    std::vector<BasePtr> m_ops;

};

}