1 /** @file gsSumOp.h
2 
3     @brief Provides the sum of \a gsLinearOperator s as a \a gsLinearOperator.
4 
5     This file is part of the G+Smo library.
6 
7     This Source Code Form is subject to the terms of the Mozilla Public
8     License, v. 2.0. If a copy of the MPL was not distributed with this
9     file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11     Author(s): S. Takacs
12 */
13 #pragma once
14 
15 #include <gsSolver/gsLinearOperator.h>
16 
17 namespace gismo
18 {
19 
20 /// @brief Class for representing the sum of objects of type \a gsLinearOperator as \a gsLinearOperator
21 ///
22 /// @ingroup Solver
23 template<typename T>
24 class gsSumOp GISMO_FINAL : public gsLinearOperator<T>
25 {
26     typedef typename gsLinearOperator<T>::Ptr BasePtr;
27 public:
28 
29     /// Shared pointer for gsSumOp
30     typedef memory::shared_ptr<gsSumOp> Ptr;
31 
32     /// Unique pointer for gsSumOp
33     typedef memory::unique_ptr<gsSumOp> uPtr;
34 
35     /// Empty constructor. To be filled with addOperator()
gsSumOp()36     gsSumOp() : m_ops(0) {}
37 
38     /// Constructor taking a vector of Linear Operators
gsSumOp(std::vector<BasePtr> ops)39     gsSumOp(std::vector<BasePtr> ops)
40         : m_ops(give(ops))
41     {
42 #ifndef NDEBUG
43         const size_t sz = m_ops.size();
44         for (size_t i=0; i<sz; ++i)
45         {
46             GISMO_ASSERT ( m_ops[0]->rows() == m_ops[i]->rows() && m_ops[0]->cols() == m_ops[i]->cols(),
47                            "Dimensions of the operators do not fit." );
48         }
49 #endif
50     }
51 
52     /// Constructor taking two Linear Operators
gsSumOp(BasePtr op0,BasePtr op1)53     gsSumOp(BasePtr op0, BasePtr op1)
54         : m_ops(2)
55     {
56         GISMO_ASSERT ( op0->rows() == op1->rows() && op0->cols() == op1->cols(), "Dimensions of the operators do not fit." );
57         m_ops[0] = give(op0); m_ops[1] = give(op1);
58     }
59 
60     /// Constructor taking three Linear Operators
gsSumOp(BasePtr op0,BasePtr op1,BasePtr op2)61     gsSumOp(BasePtr op0, BasePtr op1, BasePtr op2 )
62         : m_ops(3)
63     {
64         GISMO_ASSERT ( op0->rows() == op1->rows() && op0->cols() == op1->cols()
65                         && op0->rows() == op2->rows() && op0->cols() == op2->cols(), "Dimensions of the operators do not fit." );
66         m_ops[0] = give(op0); m_ops[1] = give(op1); m_ops[2] = give(op2);
67     }
68 
69     /// Make command returning a smart pointer
make()70     static uPtr make()
71     { return uPtr( new gsSumOp() ); }
72 
73     /// Make command returning a smart pointer
make(std::vector<BasePtr> ops)74     static uPtr make(std::vector<BasePtr> ops)
75     { return uPtr( new gsSumOp(give(ops)) ); }
76 
77     /// Make command returning a smart pointer
make(BasePtr op0,BasePtr op1)78     static uPtr make(BasePtr op0, BasePtr op1)
79     { return uPtr( new gsSumOp(give(op0),give(op1)) ); }
80 
81     /// Make command returning a smart pointer
make(BasePtr op0,BasePtr op1,BasePtr op2)82     static uPtr make(BasePtr op0, BasePtr op1, BasePtr op2)
83     { return uPtr( new gsSumOp(give(op0),give(op1),give(op2)) ); }
84 
85     /// Add another operator
addOperator(BasePtr op)86     void addOperator(BasePtr op)
87     {
88         GISMO_ASSERT ( m_ops.empty() || ( op->rows() == m_ops[0]->rows() && op->cols() == m_ops[0]->cols() ),
89                        "Dimensions of the operators do not fit." );
90         m_ops.push_back(give(op));
91     }
92 
apply(const gsMatrix<T> & input,gsMatrix<T> & x)93     void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
94     {
95         GISMO_ASSERT ( !m_ops.empty(), "gsSumOp::apply does not work for 0 operators." );
96 
97         // Here, we could make a permanently allocated vector
98         gsMatrix<T> tmp;
99         const size_t sz = m_ops.size();
100 
101         m_ops[0]->apply(input,x);
102         for (size_t i=1; i<sz; ++i)
103         {
104             m_ops[i]->apply(input,tmp);
105             x += tmp;
106         }
107     }
108 
rows()109     index_t rows() const
110     {
111         GISMO_ASSERT( !m_ops.empty(), "gsSumOp::rows does not work for 0 operators." );
112         return m_ops[0]->rows();
113     }
114 
cols()115     index_t cols() const
116     {
117         GISMO_ASSERT( !m_ops.empty(), "gsSumOp::cols does not work for 0 operators." );
118         return m_ops[0]->cols();
119     }
120 
121     /// Return a vector of shared pointers to all operators
getOps()122     const std::vector<BasePtr>& getOps() const { return m_ops; }
123 
124 private:
125     std::vector<BasePtr> m_ops;
126 
127 };
128 
129 }
130