Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
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 :
34 : // include bridge
35 : #include "bridge/bridge.h"
36 : #include "bridge/util.h"
37 : #include "bridge/util_algebra_dependent.h"
38 :
39 : #ifndef BRIDGE_MAT_VEC_OPERATIONS_H_
40 : #define BRIDGE_MAT_VEC_OPERATIONS_H_
41 :
42 : namespace ug{
43 :
44 : namespace bridge{
45 : namespace AlgebraCommon{
46 :
47 : template<typename TAlgebra>
48 0 : class VecScaleAddClass
49 : {
50 : typedef typename TAlgebra::vector_type vector_type;
51 : public:
52 0 : VecScaleAddClass(double scale1, SmartPtr<vector_type> v1)
53 : {
54 0 : scaling.push_back(scale1);
55 0 : vecs.push_back(v1);
56 0 : }
57 0 : VecScaleAddClass(double scale1, SmartPtr<vector_type> v1, double scale2, SmartPtr<vector_type> v2)
58 : {
59 0 : scaling.push_back(scale1);
60 0 : vecs.push_back(v1);
61 0 : scaling.push_back(scale2);
62 0 : vecs.push_back(v2);
63 0 : }
64 :
65 0 : VecScaleAddClass(double scale, SmartPtr<VecScaleAddClass<TAlgebra > > vsac, double scale1, SmartPtr<vector_type> v1)
66 : {
67 0 : scaling.push_back(scale1);
68 0 : vecs.push_back(v1);
69 0 : for(size_t i=0; i<vsac->size(); i++)
70 : {
71 0 : scaling.push_back(scale*vsac->scaling[i]);
72 0 : vecs.push_back(vsac->vecs[i]);
73 : }
74 0 : }
75 :
76 0 : VecScaleAddClass(double scale1, SmartPtr<vector_type> v1, double scale, SmartPtr<VecScaleAddClass<TAlgebra > > vsac)
77 : {
78 0 : scaling.push_back(scale1);
79 0 : vecs.push_back(v1);
80 0 : for(size_t i=0; i<vsac->size(); i++)
81 : {
82 0 : scaling.push_back(scale*vsac->scaling[i]);
83 0 : vecs.push_back(vsac->vecs[i]);
84 : }
85 0 : }
86 :
87 0 : VecScaleAddClass(double scale, SmartPtr<VecScaleAddClass<TAlgebra > > vsac)
88 : {
89 0 : for(size_t i=0; i<vsac->size(); i++)
90 : {
91 0 : scaling.push_back(scale*vsac->scaling[i]);
92 0 : vecs.push_back(vsac->vecs[i]);
93 : }
94 0 : }
95 :
96 0 : VecScaleAddClass() {}
97 :
98 : VecScaleAddClass(const VecScaleAddClass<TAlgebra> &parent) :
99 : scaling(parent.scaling), vecs(parent.vecs)
100 : { }
101 :
102 :
103 : size_t size() const { return scaling.size(); }
104 :
105 0 : SmartPtr<vector_type> eval()
106 : {
107 0 : SmartPtr<vector_type> p = vecs[0]->clone();
108 0 : assign(p);
109 0 : return p;
110 : }
111 :
112 0 : void assign(SmartPtr<vector_type> p)
113 : {
114 : #ifdef UG_PARALLEL
115 : for(size_t i=0; i<size(); i++)
116 : vecs[i]->change_storage_type(PST_CONSISTENT);
117 : #endif
118 0 : if(size() == 1)
119 0 : VecScaleAssign(*p, scaling[0], *vecs[0]);
120 0 : else if(size() == 2)
121 0 : VecScaleAdd(*p, scaling[0], *vecs[0], scaling[1], *vecs[1]);
122 0 : else if(size() == 3)
123 0 : VecScaleAdd(*p, scaling[0], *vecs[0], scaling[1], *vecs[1], scaling[2], *vecs[2]);
124 : else
125 : {
126 0 : UG_THROW("not supported ATM.");
127 : }
128 0 : }
129 : private:
130 : std::vector<double> scaling;
131 : std::vector<SmartPtr<vector_type> > vecs;
132 : };
133 :
134 : template<typename TAlgebra>
135 0 : void Assign(SmartPtr<typename TAlgebra::vector_type> p, SmartPtr<VecScaleAddClass<TAlgebra> > vsac)
136 : {
137 0 : vsac->assign(p);
138 0 : }
139 :
140 : template<typename TAlgebra>
141 0 : SmartPtr<typename TAlgebra::vector_type> Eval(SmartPtr<VecScaleAddClass<TAlgebra> > vsac)
142 : {
143 0 : return vsac->eval();
144 : }
145 :
146 :
147 : template<typename TAlgebra>
148 0 : double VecProd2(SmartPtr<typename TAlgebra::vector_type> v1, SmartPtr<typename TAlgebra::vector_type> v2)
149 : {
150 0 : return v1->dotprod(*v2);
151 : }
152 :
153 : template<typename T>
154 0 : SmartPtr<T> VecProdOp(SmartPtr< ILinearOperator<T, T> > op, SmartPtr<T> v)
155 : {
156 0 : SmartPtr<T> v2 = v->clone();
157 : #ifdef UG_PARALLEL
158 : v->change_storage_type(PST_CONSISTENT);
159 : #endif
160 0 : op->apply(*v2, *v);
161 0 : return v2;
162 : }
163 :
164 :
165 : template<typename TAlgebra>
166 0 : double VecScaleAddProd1(SmartPtr<typename TAlgebra::vector_type> v1, SmartPtr<VecScaleAddClass<TAlgebra> > vsac)
167 : {
168 0 : SmartPtr<typename TAlgebra::vector_type> v2 = vsac->eval();
169 0 : return VecProd2<TAlgebra>(v1, v2);
170 : }
171 :
172 : template<typename TAlgebra>
173 0 : double VecScaleAddProd2(SmartPtr<VecScaleAddClass<TAlgebra> > vsac, SmartPtr<typename TAlgebra::vector_type> v1)
174 : {
175 0 : SmartPtr<typename TAlgebra::vector_type> v2 = vsac->eval();
176 0 : return VecProd2<TAlgebra>(v1, v2);
177 : }
178 :
179 : template<typename TAlgebra>
180 0 : double VecNorm(SmartPtr<typename TAlgebra::vector_type> v)
181 : {
182 0 : return v->norm();
183 : }
184 :
185 :
186 : template<typename TAlgebra>
187 0 : double VecMaxNorm(SmartPtr<typename TAlgebra::vector_type> v)
188 : {
189 0 : return v->maxnorm();
190 : }
191 :
192 :
193 : template<typename TAlgebra>
194 0 : double VecScaleAddNorm(SmartPtr<VecScaleAddClass<TAlgebra> > vsac)
195 : {
196 0 : SmartPtr<typename TAlgebra::vector_type> v = vsac->eval();
197 0 : return v->norm();
198 : }
199 :
200 :
201 :
202 : //! calculates dest = alpha1*v1 + alpha2*v2
203 : template<typename vector_t>
204 0 : inline void VecScaleAdd2(vector_t &dest, double alpha1, const vector_t &v1, double alpha2, const vector_t &v2)
205 : {
206 0 : VecScaleAdd(dest, alpha1, v1, alpha2, v2);
207 0 : }
208 :
209 : //! calculates dest = alpha1*v1 + alpha2*v2 + alpha3*v3
210 : template<typename vector_t>
211 0 : inline void VecScaleAdd3(vector_t &dest, double alpha1, const vector_t &v1, double alpha2, const vector_t &v2, double alpha3, const vector_t &v3)
212 : {
213 0 : VecScaleAdd(dest, alpha1, v1, alpha2, v2, alpha3, v3);
214 0 : }
215 :
216 : }
217 : }
218 : }
219 : #endif /* BRIDGE_MAT_VEC_OPERATIONS_H_ */
|
