Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Nägel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__PRODUCT__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__PRODUCT__
35 :
36 : #include "lib_algebra/operator/interface/linear_iterator.h"
37 : #include "common/util/smart_pointer.h"
38 : #include <vector>
39 :
40 : #ifdef UG_PARALLEL
41 : #include "lib_algebra/parallelization/parallelization.h"
42 : #endif
43 :
44 : namespace ug{
45 :
46 : /** Base class for ILinearIterators build from other ILinearIterators */
47 : template <typename X, typename Y>
48 : class CombinedLinearIterator : public ILinearIterator<X,Y>
49 : {
50 : public:
51 0 : CombinedLinearIterator() {};
52 :
53 0 : CombinedLinearIterator(const std::vector<SmartPtr<ILinearIterator<X,Y> > >& vIterator)
54 0 : : m_vIterator(vIterator)
55 0 : {};
56 :
57 : // Name of Iterator
58 : virtual const char* name() const = 0;
59 :
60 : /// returns if parallel solving is supported
61 0 : virtual bool supports_parallel() const
62 : {
63 0 : for(size_t i = 0; i < m_vIterator.size(); ++i)
64 0 : if(!m_vIterator[i]->supports_parallel())
65 : return false;
66 : return true;
67 : }
68 :
69 : // Prepare for Operator J(u) and linearization point u (current solution)
70 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > J, const Y& u) = 0;
71 :
72 : // Prepare for Linear Operator L
73 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > L) = 0;
74 :
75 : // Compute correction
76 : virtual bool apply(Y& c, const X& d) = 0;
77 :
78 : // Compute correction and update defect
79 : virtual bool apply_update_defect(Y& c, X& d) = 0;
80 :
81 0 : void add_iterator(SmartPtr<ILinearIterator<X,Y> > I) {m_vIterator.push_back(I);}
82 :
83 0 : void add_iterator(SmartPtr<ILinearIterator<X,Y> > I,size_t nr) { for (size_t i=0;i<nr;i++) m_vIterator.push_back(I);}
84 :
85 : // Clone
86 : virtual SmartPtr<ILinearIterator<X,Y> > clone() = 0;
87 :
88 : protected:
89 : std::vector<SmartPtr<ILinearIterator<X,Y> > > m_vIterator;
90 : };
91 :
92 :
93 : /** This operator is a product of ILinearIterator (multiplicative composition).*/
94 : template <typename X, typename Y>
95 0 : class LinearIteratorProduct : public CombinedLinearIterator<X,Y>
96 : {
97 : protected:
98 : typedef CombinedLinearIterator<X,Y> base_type;
99 : using base_type::m_vIterator;
100 :
101 : public:
102 0 : LinearIteratorProduct() {};
103 :
104 0 : LinearIteratorProduct(const std::vector<SmartPtr<ILinearIterator<X,Y> > >& vIterator)
105 0 : : CombinedLinearIterator<X,Y>(vIterator)
106 : {};
107 :
108 : // Name of Iterator
109 0 : virtual const char* name() const {return "IteratorProduct";}
110 :
111 : // Prepare for Operator J(u) and linearization point u (current solution)
112 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > J, const Y& u) {
113 : bool bRes = true;
114 0 : for(size_t i = 0; i < m_vIterator.size(); i++){
115 0 : if ((i>0) && (m_vIterator[i]==m_vIterator[i-1])) continue;
116 0 : bRes &= m_vIterator[i]->init(J,u);
117 : }
118 0 : return bRes;
119 : }
120 :
121 : // Prepare for Linear Operator L
122 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > L)
123 : {
124 : bool bRes = true;
125 0 : for(size_t i = 0; i < m_vIterator.size(); i++) {
126 0 : if ((i>0) && (m_vIterator[i]==m_vIterator[i-1])) continue;
127 0 : bRes &= m_vIterator[i]->init(L);
128 : }
129 0 : return bRes;
130 : }
131 :
132 : // Compute correction
133 0 : virtual bool apply(Y& c, const X& d)
134 : {
135 : // create temporary defect and forward request
136 0 : SmartPtr<X> spDTmp = d.clone();
137 0 : return apply_update_defect(c, *spDTmp);
138 : }
139 :
140 : // Compute correction and update defect
141 0 : virtual bool apply_update_defect(Y& c, X& d)
142 : {
143 0 : SmartPtr<Y> spCTmp = c.clone_without_values();
144 :
145 : bool bRes = true;
146 : c.set(0.0);
147 0 : for(size_t i = 0; i < m_vIterator.size(); i++)
148 : {
149 0 : bRes &= m_vIterator[i]->apply_update_defect(*spCTmp,d);
150 0 : c += (*spCTmp);
151 : }
152 0 : return bRes;
153 : }
154 :
155 : // Clone
156 0 : virtual SmartPtr<ILinearIterator<X,Y> > clone() {
157 0 : return make_sp(new LinearIteratorProduct(*this));
158 : }
159 : };
160 :
161 : /** This operator is a sum of ILinearIterator (additive composition). */
162 : template <typename X, typename Y>
163 : class LinearIteratorSum : public CombinedLinearIterator<X,Y>
164 : {
165 : protected:
166 : typedef CombinedLinearIterator<X,Y> base_type;
167 : using base_type::m_vIterator;
168 :
169 : public:
170 0 : LinearIteratorSum() {};
171 :
172 0 : LinearIteratorSum(const std::vector<SmartPtr<ILinearIterator<X,Y> > >& vIterator)
173 0 : : CombinedLinearIterator<X,Y>(vIterator)
174 : {};
175 :
176 : // Name of Iterator
177 0 : virtual const char* name() const {return "IteratorProduct";}
178 :
179 : // Prepare for Operator J(u) and linearization point u (current solution)
180 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > J, const Y& u) {
181 0 : m_spOp = J;
182 : bool bRes = true;
183 0 : for(size_t i = 0; i < m_vIterator.size(); i++){
184 0 : if ((i>0) && (m_vIterator[i]==m_vIterator[i-1])) continue;
185 0 : bRes &= m_vIterator[i]->init(J,u);
186 : }
187 0 : return bRes;
188 : }
189 :
190 : // Prepare for Linear Operator L
191 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > L)
192 : {
193 0 : m_spOp = L;
194 : bool bRes = true;
195 0 : for(size_t i = 0; i < m_vIterator.size(); i++) {
196 0 : if ((i>0) && (m_vIterator[i]==m_vIterator[i-1])) continue;
197 0 : bRes &= m_vIterator[i]->init(L);
198 : }
199 0 : return bRes;
200 : }
201 :
202 : // Compute correction
203 0 : virtual bool apply(Y& c, const X& d)
204 : {
205 : // create temporary correction
206 0 : SmartPtr<Y> spCTmp = c.clone_without_values();
207 :
208 : // this part computes all corrections independently
209 : c.set(0.0);
210 :
211 : bool bRes = true;
212 0 : for(size_t i = 0; i < m_vIterator.size(); i++)
213 : {
214 0 : bRes &= m_vIterator[i]->apply(*spCTmp,d);
215 0 : c += (*spCTmp);
216 : }
217 :
218 0 : return bRes;
219 : }
220 :
221 : // Compute correction and update defect
222 0 : virtual bool apply_update_defect(Y& c, X& d)
223 : {
224 0 : bool bRet = apply(c, d);
225 0 : m_spOp->apply_sub(d, c);
226 0 : return bRet;
227 : }
228 :
229 : // Clone
230 0 : virtual SmartPtr<ILinearIterator<X,Y> > clone() {
231 0 : return make_sp(new LinearIteratorSum(*this));
232 : }
233 :
234 : protected:
235 : SmartPtr<ILinearOperator<Y,X> > m_spOp;
236 : };
237 :
238 :
239 : } // end namespace ug
240 :
241 : #endif
|
