Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : #ifndef __H__UG__LIB_ALGEBRA__OPERATIONS_VEC__
34 : #define __H__UG__LIB_ALGEBRA__OPERATIONS_VEC__
35 :
36 : #include <stddef.h> // size_t
37 : #include <cmath> // log, exp
38 :
39 : namespace ug
40 : {
41 :
42 : // operations for doubles
43 : //-----------------------------------------------------------------------------
44 : // todo: check if we might replace double with template<T>
45 :
46 : // VecScale: These function calculate dest = sum_i alpha_i v_i
47 :
48 : //! calculates dest = alpha1*v1. for doubles
49 : inline void VecScaleAssign(double &dest, double alpha1, const double &v1)
50 : {
51 0 : dest = alpha1*v1;
52 : }
53 :
54 : //! calculates dest = alpha1*v1 + alpha2*v2. for doubles
55 : inline void VecScaleAdd(double &dest, double alpha1, const double &v1, double alpha2, const double &v2)
56 : {
57 0 : dest = alpha1*v1 + alpha2*v2;
58 : }
59 :
60 : //! calculates dest = alpha1*v1 + alpha2*v2 + alpha3*v3. for doubles
61 : inline void VecScaleAdd(double &dest, double alpha1, const double &v1, double alpha2, const double &v2, double alpha3, const double &v3)
62 : {
63 0 : dest = alpha1*v1 + alpha2*v2 + alpha3*v3;
64 : }
65 :
66 :
67 : // VecProd
68 :
69 : //! calculates s += scal<a, b>
70 : inline void VecProdAdd(const double &a, const double &b, double &s)
71 : {
72 0 : s += a*b;
73 : }
74 :
75 : template<typename vector_t>
76 : inline void VecProdAdd(const vector_t &a, const vector_t &b, double &s)
77 : {
78 0 : for(size_t i=0; i<a.size(); i++)
79 : VecProdAdd(a[i], b[i], s);
80 : }
81 :
82 :
83 : //! returns scal<a, b>
84 : inline double VecProd(const double &a, const double &b)
85 : {
86 0 : return a*b;
87 : }
88 :
89 :
90 : //! computes scal<a, b>
91 : inline void VecProd(const double &a, const double &b, double &s)
92 : {
93 : s = a*b;
94 : }
95 :
96 :
97 : // VecNorm
98 :
99 : //! returns norm_2^2(a)
100 : inline double VecNormSquared(const double &a)
101 : {
102 : return a*a;
103 : }
104 :
105 : //! calculates s += norm_2^2(a)
106 : inline void VecNormSquaredAdd(const double &a, double &s)
107 : {
108 0 : s += a*a;
109 : }
110 :
111 : // Elementwise (Hadamard) product of two vectors
112 :
113 : //! calculates s = a * b (the Hadamard product)
114 : inline void VecHadamardProd(double &dest, const double &v1, const double &v2)
115 : {
116 0 : dest = v1 * v2;
117 : }
118 :
119 : // Some unary mathematical elementwise operations on vectors
120 :
121 : //! calculates elementwise exp
122 : inline void VecExp(double &dest, const double &v)
123 : {
124 0 : dest = exp (v);
125 : }
126 :
127 : //! calculates elementwise log (natural logarithm)
128 : inline void VecLog(double &dest, const double &v)
129 : {
130 0 : dest = log (v);
131 : }
132 :
133 : // templated
134 :
135 : // operations for vectors
136 : //-----------------------------------------------------------------------------
137 : // these functions execute vector operations by using the operations on the elements of the vector
138 :
139 : // todo: change vector_t to TE_VEC<vector_t>
140 :
141 :
142 : // VecScale: These function calculate dest = sum_i alpha_i v_i
143 :
144 : //! calculates dest = alpha1*v1
145 : template<typename vector_t>
146 0 : inline void VecScaleAssign(vector_t &dest, double alpha1, const vector_t &v1)
147 : {
148 0 : for(size_t i=0; i<dest.size(); i++)
149 : VecScaleAssign(dest[i], alpha1, v1[i]);
150 0 : }
151 :
152 : //! sets dest = v1 entrywise
153 : template<typename vector_t>
154 0 : inline void VecAssign(vector_t &dest, const vector_t &v1)
155 : {
156 0 : for(size_t i=0; i<dest.size(); i++)
157 0 : dest[i] = v1[i];
158 0 : }
159 :
160 : //! calculates dest = alpha1*v1 + alpha2*v2
161 : template<typename vector_t, template <class T> class TE_VEC>
162 0 : inline void VecScaleAdd(TE_VEC<vector_t> &dest, double alpha1, const TE_VEC<vector_t> &v1, double alpha2, const TE_VEC<vector_t> &v2)
163 : {
164 0 : for(size_t i=0; i<dest.size(); i++)
165 : VecScaleAdd(dest[i], alpha1, v1[i], alpha2, v2[i]);
166 0 : }
167 :
168 : //! calculates dest = alpha1*v1 + alpha2*v2 + alpha3*v3
169 : template<typename vector_t, template <class T> class TE_VEC>
170 0 : inline void VecScaleAdd(TE_VEC<vector_t> &dest, double alpha1, const TE_VEC<vector_t> &v1, double alpha2, const TE_VEC<vector_t> &v2, double alpha3, const TE_VEC<vector_t> &v3)
171 : {
172 0 : for(size_t i=0; i<dest.size(); i++)
173 : VecScaleAdd(dest[i], alpha1, v1[i], alpha2, v2[i], alpha3, v3[i]);
174 0 : }
175 :
176 :
177 : // VecProd
178 :
179 : //! calculates s += scal<a, b>
180 : template<typename vector_t>
181 : inline void VecProd(const vector_t &a, const vector_t &b, double &sum)
182 : {
183 : for(size_t i=0; i<a.size(); i++)
184 : VecProdAdd(a[i], b[i], sum);
185 : }
186 :
187 : //! returns scal<a, b>
188 : template<typename vector_t>
189 : inline double VecProd(const vector_t &a, const vector_t &b)
190 : {
191 : double sum=0;
192 : VecProdAdd(a, b, sum);
193 : return sum;
194 : }
195 :
196 :
197 : //! calculates s += norm_2^2(a)
198 : template<typename vector_t>
199 : inline void VecNormSquaredAdd(const vector_t &a, double &sum)
200 : {
201 0 : for(size_t i=0; i<a.size(); i++)
202 : VecNormSquaredAdd(a[i], sum);
203 : }
204 :
205 : //! returns norm_2^2(a)
206 : template<typename vector_t>
207 : inline double VecNormSquared(const vector_t &a)
208 : {
209 : double sum=0;
210 : VecNormSquaredAdd(a, sum);
211 : return sum;
212 : }
213 :
214 : // Elementwise (Hadamard) product of two vectors
215 : template<typename vector_t>
216 0 : inline void VecHadamardProd(vector_t &dest, const vector_t &v1, const vector_t &v2)
217 : {
218 0 : for(size_t i=0; i<dest.size(); i++)
219 : VecHadamardProd(dest[i], v1[i], v2[i]);
220 0 : }
221 :
222 : // Elementwise exp on a vector
223 : template<typename vector_t>
224 0 : inline void VecExp(vector_t &dest, const vector_t &v)
225 : {
226 0 : for(size_t i=0; i<dest.size(); i++)
227 : VecExp(dest[i], v[i]);
228 0 : }
229 :
230 : // Elementwise log (natural logarithm) on a vector
231 : template<typename vector_t>
232 0 : inline void VecLog(vector_t &dest, const vector_t &v)
233 : {
234 0 : for(size_t i=0; i<dest.size(); i++)
235 : VecLog(dest[i], v[i]);
236 0 : }
237 :
238 : } // namespace ug
239 :
240 : #endif /* __H__UG__LIB_ALGEBRA__OPERATIONS_VEC__ */
|
