1 /** @file gsLinearOperator.h
2 
3     @brief Simple abstract class for (discrete) linear operators.
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): J. Sogn, A. Manzaflaris, C. Hofreither, S. Takacs, C. Hofer
12 */
13 #pragma once
14 
15 #include <gsCore/gsLinearAlgebra.h>
16 #include <gsIO/gsOptionList.h>
17 
18 namespace gismo
19 {
20 
21 /// @brief Simple abstract class for discrete operators.
22 ///
23 /// Simple abstract class for discrete operators.
24 /// The derived classes have to contain the functions: apply(), cols(), and rows().
25 ///
26 /// \ingroup Solver
27 template<class T>
28 class gsLinearOperator
29 {
30 public:
31 
32     /// Shared pointer for gsLinearOperator
33     typedef memory::shared_ptr<gsLinearOperator> Ptr;
34 
35     /// Unique pointer for gsLinearOperator
36     typedef memory::unique_ptr<gsLinearOperator> uPtr;
37 
38     /// Identity operator
Identity(const index_t dim)39     static gsIdentityOp<T> Identity(const index_t dim) {return gsIdentityOp<T>(dim); }
40 
~gsLinearOperator()41     virtual ~gsLinearOperator() {}
42 
43     /**
44      * @brief apply the operator on the input vector and store the result in x
45      * @param input Input vector
46      * @param x     result vector
47      */
48     virtual void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const = 0;
49 
50     /// Returns the number of rows of the operator
51     virtual index_t rows() const = 0;
52 
53     /// Returns the number of columns of the operator
54     virtual index_t cols() const = 0;
55 
56     // NOTE: this is rather inefficient and is only provided for debugging and testing purposes
toMatrix(gsMatrix<T> & result)57     void toMatrix(gsMatrix<T>& result)
58     {
59         gsMatrix<T> eye = gsMatrix<T>::Identity(cols(), cols());
60         this->apply(eye, result);
61     }
62 
63     /// Get the default options as gsOptionList object
64     // This implementation provides an empty object
defaultOptions()65     static gsOptionList defaultOptions()               { return gsOptionList(); }
66 
67     /// Set options based on a gsOptionList object
68     // This implementation does not read any input
setOptions(const gsOptionList &)69     virtual void setOptions(const gsOptionList &)  {}
70 
71 }; // gsLinearOperator
72 
73 /// @brief Allows an operator to be multiplied with a scalar
74 ///
75 /// \ingroup Solver
76 template<class T>
77 class gsScaledOp GISMO_FINAL : public gsLinearOperator<T>
78 {
79 public:
80     /// Shared pointer for gsScaledOp
81     typedef memory::shared_ptr<gsScaledOp> Ptr;
82 
83     /// Unique pointer for gsScaledOp
84     typedef memory::unique_ptr<gsScaledOp> uPtr;
85 
86     /// Shared pointer for gsLinearOperator
87     typedef typename gsLinearOperator<T>::Ptr BasePtr;
88 
89     /// Constructor taking a shared pointer to a linear operator and a scalar
m_op(give (op))90     gsScaledOp(BasePtr op, T scalar = 1) : m_op(give(op)), m_scalar(scalar)    {}
91 
92     /// Make function returning a smart pointer
93     static uPtr make(BasePtr op, T scalar = 1)
94     { return uPtr( new gsScaledOp(give(op), scalar) ); }
95 
apply(const gsMatrix<T> & input,gsMatrix<T> & x)96     virtual void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
97     {
98         m_op->apply(input, x);
99         x *= m_scalar;
100     }
101 
102     ///Returns the number of rows in the preconditioner
rows()103     index_t rows() const {return m_op->rows();}
104 
105     ///Returns the number of columns in the preconditioner
cols()106     index_t cols() const {return m_op->cols();}
107 
108 private:
109     const BasePtr m_op;
110     const T m_scalar;
111 }; // gsScaladOp
112 
113 
114 /// @brief Identity operator
115 ///
116 /// \ingroup Solver
117 template<class T>
118 class gsIdentityOp GISMO_FINAL : public gsLinearOperator<T>
119 {
120 public:
121 
122     /// Shared pointer for gsIdentityOp
123     typedef memory::shared_ptr<gsIdentityOp> Ptr;
124 
125     /// Unique pointer for gsIdentityOp
126     typedef memory::unique_ptr<gsIdentityOp> uPtr;
127 
128     /// Constructor taking the dimension of the identity operator
gsIdentityOp(index_t dim)129     gsIdentityOp(index_t dim) : m_dim(dim) {}
130 
131     /// Make function returning a smart pointer
make(index_t dim)132     static uPtr make(index_t dim) { return uPtr( new gsIdentityOp(dim) ); }
133 
apply(const gsMatrix<T> & input,gsMatrix<T> & x)134     void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
135     {
136         x = input;
137     }
138 
rows()139     index_t rows() const {return m_dim;}
140 
cols()141     index_t cols() const {return m_dim;}
142 
143 private:
144     const index_t m_dim;
145 };
146 
147 } // namespace gismo
148