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 :
34 : #ifndef __H__UG__CPU_ALGEBRA__VECTOR_IMPL__
35 : #define __H__UG__CPU_ALGEBRA__VECTOR_IMPL__
36 :
37 : #include <fstream>
38 : #include <algorithm>
39 : #include "algebra_misc.h"
40 : #include "common/math/ugmath.h"
41 : #include "vector.h" // for urand
42 :
43 : #define prefetchReadWrite(a)
44 :
45 : namespace ug{
46 : template<typename value_type>
47 : inline value_type &Vector<value_type>::operator [] (size_t i)
48 : {
49 : UG_ASSERT(i < m_size, *this << ": tried to access element " << i);
50 0 : return values[i];
51 : }
52 :
53 : template<typename value_type>
54 : inline const value_type &Vector<value_type>::operator [] (size_t i) const
55 : {
56 : UG_ASSERT(i < m_size, *this << ": tried to access element " << i);
57 0 : return values[i];
58 : }
59 :
60 :
61 : // energynorm2 = x*(A*x)
62 : /*inline double Vector<value_type>::energynorm2(const SparseMatrix &A) const
63 : {
64 : double sum=0;
65 : for(size_t i=0; i<m_size; i++) sum += (A[i] * (*this)) * values[i];
66 : //FOR_UNROLL_FWD(i, 0, m_size, UNROLL, sum += A[i] * (*this) * values[i]);
67 : return sum;
68 : }*/
69 :
70 : // dotprod
71 : template<typename value_type>
72 : inline double Vector<value_type>::dotprod(const Vector &w) //const
73 : {
74 : UG_ASSERT(m_size == w.m_size, *this << " has not same size as " << w);
75 :
76 : double sum=0;
77 0 : for(size_t i=0; i<m_size; i++) sum += VecProd(values[i], w[i]);
78 : return sum;
79 : }
80 :
81 : // assign double to whole Vector
82 : template<typename value_type>
83 : inline double Vector<value_type>::operator = (double d)
84 : {
85 0 : for(size_t i=0; i<m_size; i++)
86 0 : values[i] = d;
87 : return d;
88 : }
89 :
90 : template<typename value_type>
91 0 : inline void Vector<value_type>::set_random(double from, double to)
92 : {
93 0 : for(size_t i=0; i<size(); i++)
94 0 : for(size_t j=0; j<GetSize(values[i]); j++)
95 0 : BlockRef(values[i], j) = urand(from, to);
96 0 : }
97 :
98 :
99 : template<typename value_type>
100 0 : inline void Vector<value_type>::operator = (const vector_type &v)
101 : {
102 : resize(v.size());
103 0 : for(size_t i=0; i<m_size; i++)
104 0 : values[i] = v[i];
105 0 : }
106 :
107 : template<typename value_type>
108 0 : inline void Vector<value_type>::operator += (const vector_type &v)
109 : {
110 : UG_ASSERT(v.size() == size(), "vector sizes must match! (" << v.size() << " != " << size() << ")");
111 0 : for(size_t i=0; i<m_size; i++)
112 0 : values[i] += v[i];
113 0 : }
114 :
115 : template<typename value_type>
116 0 : inline void Vector<value_type>::operator -= (const vector_type &v)
117 : {
118 : UG_ASSERT(v.size() == size(), "vector sizes must match! (" << v.size() << " != " << size() << ")");
119 0 : for(size_t i=0; i<m_size; i++)
120 0 : values[i] -= v[i];
121 0 : }
122 :
123 : ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
124 :
125 :
126 : template<typename value_type>
127 0 : Vector<value_type>::Vector () : m_size(0), m_capacity(0), values(NULL)
128 : {
129 : FORCE_CREATION { p(); } // force creation of this rountines for gdb.
130 : }
131 :
132 : template<typename value_type>
133 0 : Vector<value_type>::Vector(size_t size) : m_size(0), m_capacity(0), values(NULL)
134 : {
135 : FORCE_CREATION { p(); } // force creation of this rountines for gdb.
136 0 : create(size);
137 0 : }
138 :
139 : template<typename value_type>
140 0 : Vector<value_type>::~Vector()
141 : {
142 : destroy();
143 0 : }
144 :
145 : template<typename value_type>
146 : void Vector<value_type>::destroy()
147 : {
148 0 : if(values)
149 : {
150 0 : delete [] values;
151 0 : values = NULL;
152 : }
153 : m_size = 0;
154 : }
155 :
156 :
157 : template<typename value_type>
158 0 : void Vector<value_type>::create(size_t size)
159 : {
160 0 : if(m_size == size) return;
161 : destroy();
162 0 : m_size = size;
163 0 : values = new value_type[size];
164 0 : m_capacity = size;
165 : }
166 :
167 :
168 : template<typename value_type>
169 0 : Vector<value_type>* Vector<value_type>::virtual_clone() const
170 : {
171 0 : return new Vector<value_type>(*this);
172 : }
173 :
174 : template<typename value_type>
175 0 : SmartPtr<Vector<value_type> > Vector<value_type>::clone() const
176 : {
177 0 : return SmartPtr<Vector<value_type> >(this->virtual_clone());
178 : }
179 :
180 : template<typename value_type>
181 0 : Vector<value_type>* Vector<value_type>::virtual_clone_without_values() const
182 : {
183 0 : return new Vector<value_type>(this->m_size);
184 : }
185 :
186 : template<typename value_type>
187 0 : SmartPtr<Vector<value_type> > Vector<value_type>::clone_without_values() const
188 : {
189 0 : return SmartPtr<Vector<value_type> >(this->virtual_clone_without_values());
190 : }
191 :
192 : template<typename value_type>
193 0 : void Vector<value_type>::reserve_exactly(size_t newCapacity, bool bCopyValues)
194 : {
195 : UG_ASSERT(newCapacity >= m_size, "use resize, then reserve_exactly");
196 0 : value_type *new_values = new value_type[newCapacity];
197 : // we cannot use memcpy here bcs of variable blocks.
198 0 : if(values != NULL && bCopyValues)
199 : {
200 0 : for(size_t i=0; i<m_size; i++)
201 0 : std::swap(new_values[i], values[i]);
202 0 : for(size_t i=m_size; i<newCapacity; i++)
203 0 : new_values[i] = 0.0;
204 : }
205 0 : if(values) delete [] values;
206 0 : values = new_values;
207 0 : m_capacity = newCapacity;
208 0 : }
209 :
210 :
211 : template<typename value_type>
212 : void Vector<value_type>::reserve_sloppy(size_t newCapacity, bool bCopyValues)
213 : {
214 : if(newCapacity <= m_capacity) return;
215 : reserve_exactly(newCapacity + m_capacity/2, bCopyValues);
216 : }
217 :
218 : template<typename value_type>
219 : void Vector<value_type>::resize_sloppy(size_t newSize, bool bCopyValues)
220 : {
221 0 : if(newSize > m_capacity)
222 : {
223 0 : size_t newCapacity = m_size/2 + newSize;
224 0 : reserve_exactly(newCapacity, true);
225 : }
226 0 : m_size = newSize;
227 : }
228 :
229 : template<typename value_type>
230 : void Vector<value_type>::resize_exactly(size_t newSize, bool bCopyValues)
231 : {
232 0 : if(newSize > m_capacity)
233 0 : reserve_exactly(newSize, true);
234 0 : m_size = newSize;
235 : }
236 :
237 :
238 :
239 : template<typename value_type>
240 : void Vector<value_type>::create(const Vector &v)
241 : {
242 : UG_ASSERT(m_size == 0, *this << " already created");
243 : m_size = v.m_size;
244 : values = new value_type[m_size];
245 : m_capacity = m_size;
246 :
247 : // we cannot use memcpy here bcs of variable blocks.
248 : for(size_t i=0; i<m_size; i++)
249 : values[i] = v.values[i];
250 : }
251 :
252 :
253 : // print
254 : template<typename value_type>
255 0 : void Vector<value_type>::print(const char * const text) const
256 : {
257 :
258 0 : if(text) std::cout << " == " << text;
259 0 : std::cout << " size: " << m_size << " =================" << std::endl;
260 0 : for(size_t i=0; i<m_size; i++)
261 : //cout << values[i] << " ";
262 0 : std::cout << i << ": " << values[i] << std::endl;
263 : std::cout << std::endl;
264 0 : }
265 :
266 :
267 : template<typename value_type>
268 : template<typename V>
269 : void Vector<value_type>::add(const V& u)
270 : {
271 : for(size_t i=0; i < u.size(); i++)
272 : values[u.index(i)] += u[i];
273 : }
274 :
275 : template<typename value_type>
276 : template<typename V>
277 : void Vector<value_type>::set(const V& u)
278 : {
279 : for(size_t i=0; i < u.size(); i++)
280 : values[u.index(i)] = u[i];
281 : }
282 :
283 : template<typename value_type>
284 : template<typename V>
285 : void Vector<value_type>::get(V& u) const
286 : {
287 : for(size_t i=0; i < u.size(); i++)
288 : u[i] = values[u.index(i)];
289 : }
290 :
291 :
292 :
293 :
294 : template<typename value_type>
295 : void Vector<value_type>::add(const value_type *u, const size_t *indices, size_t nr)
296 : {
297 : for(size_t i=0; i < nr; i++)
298 : values[indices[i]] += u[i];
299 : }
300 :
301 : template<typename value_type>
302 : void Vector<value_type>::set(const value_type *u, const size_t *indices, size_t nr)
303 : {
304 : for(size_t i=0; i < nr; i++)
305 : values[indices[i]] = u[i];
306 : }
307 :
308 : template<typename value_type>
309 : void Vector<value_type>::get(value_type *u, const size_t *indices, size_t nr) const
310 : {
311 : for(size_t i=0; i < nr; i++)
312 : u[i] = values[indices[i]] ;
313 : }
314 :
315 :
316 : template<typename value_type>
317 : double operator *(const TRANSPOSED<Vector<value_type> > &x, const Vector<value_type> &y)
318 : {
319 : return x.T().dotprod(y);
320 : }
321 :
322 : template<typename value_type>
323 0 : inline double Vector<value_type>::norm() const
324 : {
325 : double d=0;
326 0 : for(size_t i=0; i<size(); ++i)
327 0 : d+=BlockNorm2(values[i]);
328 0 : return sqrt(d);
329 : }
330 :
331 : template<typename value_type>
332 0 : inline double Vector<value_type>::maxnorm() const
333 : {
334 0 : double d=0;
335 0 : for(size_t i=0; i<size(); ++i)
336 0 : d = std::max(d, BlockMaxNorm(values[i]));
337 0 : return d;
338 : }
339 :
340 : template<typename TValueType>
341 : void CloneVector(Vector<TValueType> &dest, const Vector<TValueType>& src)
342 : {
343 : dest.resize(src.size());
344 : }
345 :
346 : template<typename TValueType, class TOStream>
347 0 : void Serialize(TOStream &buf, const Vector<TValueType> &v)
348 : {
349 0 : Serialize(buf, v.size());
350 0 : for(size_t i=0; i < v.size(); i++)
351 : Serialize(buf, v[i]);
352 0 : }
353 :
354 : template<typename TValueType, class TIStream>
355 0 : void Deserialize(TIStream &buf, Vector<TValueType> &v)
356 : {
357 : size_t s = Deserialize<size_t>(buf);
358 : v.resize(s);
359 0 : for(size_t i=0; i < s; i++)
360 : Deserialize(buf, v[i]);
361 0 : }
362 :
363 : }//namespace ug
364 :
365 : #endif
|
