b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2vector_dense.H
1//------------------------------------------------------------------------
2// b2vector_dense.H --
3//
4// A framework for generic linear algebraic (matrix) computations.
5//
6// written by Mathias Doreille
7//
8// Copyright (c) 2004-2012,2017
9// SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// All Rights Reserved. Proprietary source code. The contents of
13// this file may not be disclosed to third parties, copied or
14// duplicated in any form, in whole or in part, without the prior
15// written permission of SMR.
16//------------------------------------------------------------------------
17
18#ifndef __B2VECTOR_DENSE_H__
19#define __B2VECTOR_DENSE_H__
20
21#include "b2linear_algebra_def.H"
22#include "utils/b2csda.H"
23
24namespace b2000 { namespace b2linalg {
25
26struct Vdense {
27 typedef Vincrement_ref base;
28 typedef Vincrement_scale_constref const_base;
29 typedef Vdense copy;
30};
31
32template <typename T>
33class Vector<T, Vdense> : public std::vector<T> {
34public:
35 Vector() {}
36
37 Vector(size_t s) : std::vector<T>(s) {}
38
39 Vector(const Vector& v) : std::vector<T>(v) {}
40
41 template <typename T1>
42 Vector(const Vector<T1, Vdense>& v) : std::vector<T>(v.begin(), v.end()) {}
43
44 template <typename T1>
45 Vector(const Vector<T1, Vdense_constref>& v) : std::vector<T>(v.v, v.v + v.s) {}
46
47 template <typename T1>
48 Vector(Vector<T1, Vdense_ref>& v) : std::vector<T>(v.v, v.v + v.s) {}
49
50 /*
51 Vector(const Vector<T, Vincrement_scale_constref>& v_):
52 std::vector<T>(v_.s) {
53 const T* vi = v_.v;
54 typename std::vector<T>::iterator vo = std::vector<T>::begin();
55 typename std::vector<T>::iterator vo_end = std::vector<T>::end();
56 for (; vo != vo_end; ++vo, vi += v_.incr)
57 *vo = v_.scale * *vi;
58 }
59 */
60
61 Vector(const Vector<T, Vindex_scale_constref>& v_) : std::vector<T>(v_.s) {
62 for (size_t i = 0; i != v_.s; ++i) { (*this)[i] = v_.scale * v_.v[v_.index[i]]; }
63 }
64
65 Vector(const Vector<T, Vindex_constref>& v_) : std::vector<T>(v_.s) {
66 for (size_t i = 0; i != v_.s; ++i) { (*this)[i] = v_.v[v_.index[i]]; }
67 }
68
69 Vector(const Vector<T, Vindex_ref>& v_) : std::vector<T>(v_.s) {
70 for (size_t i = 0; i != v_.s; ++i) { (*this)[i] = v_.v[v_.index[i]]; }
71 }
72
73 /*
74 template<typename STORAGE>
75 Vector(const Vector<T, STORAGE>& v_):
76 std::vector<T>(v_.size()) {
77 Vector<T, Vincrement_ref>(*this) = v_;
78 }
79 */
80
81 bool is_null() const { return this == &null; }
82
83 void set_zero() { std::fill(std::vector<T>::begin(), std::vector<T>::end(), 0); }
84
90 void reset(const size_t size) {
91 this->resize(size);
92 std::ranges::fill(*this, T{});
93 }
94
95 template <typename STORAGE>
96 Vector& operator=(const Vector<T, STORAGE>& v_) {
97 std::vector<T>::resize(v_.size());
98 Vector<T, Vincrement_ref>(*this) = v_;
99 return *this;
100 }
101
102 template <typename STORAGE>
103 Vector& operator+=(const Vector<T, STORAGE>& v_) {
104 Vector<T, Vincrement_ref>(*this) = *this + v_;
105 return *this;
106 }
107
108 template <typename STORAGE>
109 Vector& operator-=(const Vector<T, STORAGE>& v_) {
110 Vector<T, Vincrement_ref>(*this) = *this - v_;
111 return *this;
112 }
113
114 Vector& operator*=(T scale_) {
115 *this = *this * scale_;
116 return *this;
117 }
118
119 template <typename STORAGE>
120 void scale(const Vector<T, STORAGE>& v_) {
121 if (v_.size() != this->size()) {
122 Exception() << "The two vectors do not have the same size." << THROW;
123 }
124 for (size_t i = 0; i != this->size(); ++i) { (*this)[i] *= v_[i]; }
125 }
126
127 template <typename STORAGE>
128 void scale_inv(const Vector<T, STORAGE>& v_) {
129 if (v_.size() != this->size()) {
130 Exception() << "The two vectors do not have the same size." << THROW;
131 }
132 for (size_t i = 0; i != this->size(); ++i) { (*this)[i] /= v_[i]; }
133 }
134
135 Vector<T, Vdense_ref> operator()(const Interval& i) {
136 return Vector<T, Vdense_ref>(i.size(), &(*this)[i[0]]);
137 }
138
139 Vector<T, Vdense_constref> operator()(const Interval& i) const {
140 return Vector<T, Vdense_constref>(i.size(), &(*this)[i[0]]);
141 }
142
143 Vector<T, Vincrement_ref> operator()(const Slice& i) {
144 return Vector<T, Vincrement_ref>(i.size(), &(*this)[i[0]], i.increment());
145 }
146
147 Vector<T, Vincrement_constref> operator()(const Slice& i) const {
148 return Vector<T, Vincrement_constref>(i.size(), &(*this)[i[0]], i.increment());
149 }
150
151 Vector<T, Vindex_ref> operator()(const Index& i) {
152 return Vector<T, Vindex_ref>(i.size(), &(*this)[0], &i[0]);
153 }
154
155 Vector<T, Vindex_constref> operator()(const Index& i) const {
156 return Vector<T, Vindex_constref>(i.size(), &(*this)[0], &i[0]);
157 }
158
159 void normalize() {
160 blas::scal(
161 std::vector<T>::size(), T(1) / blas::nrm2(std::vector<T>::size(), &(*this)[0], 1),
162 &(*this)[0], 1);
163 }
164
165 T norm2() const { return blas::nrm2(std::vector<T>::size(), &(*this)[0], 1); }
166
167 T norm_inf() const {
168 T res = 0;
169 for (size_t i = 0; i != std::vector<T>::size(); ++i) { res = max_abs(res, (*this)[i]); }
170 return res;
171 }
172
173 void remove(Index& idx) {
174 size_t ii = idx[0] + 1;
175 typename std::vector<T>::iterator oo = std::vector<T>::begin() + idx[0];
176 for (size_t i = 1; i < idx.size(); ++i) {
177 oo = std::copy(std::vector<T>::begin() + ii, std::vector<T>::begin() + idx[i], oo);
178 ii = idx[i] + 1;
179 }
180 oo = std::copy(std::vector<T>::begin() + ii, std::vector<T>::end(), oo);
181 std::vector<T>::erase(oo, std::vector<T>::end());
182 }
183
184 static Vector null;
185
186 friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
187 l << "vector ";
188 l.write(v.size(), &v[0], 1);
189 return l;
190 }
191
192private:
193 MVFRIEND;
194};
195
196template <typename T>
197Vector<T, Vdense> Vector<T, Vdense>::null;
198
199struct Vdense_constref {
200 typedef Vincrement_scale_constref base;
201 typedef Vincrement_scale_constref const_base;
202 typedef Vdense copy;
203};
204
205template <typename T>
206class Vector<T, Vdense_constref> {
207public:
208 Vector() : s(0), v(0) {}
209
210 Vector(size_t s_, const T* v_) : s(s_), v(v_) {}
211
212 Vector(const Vector& v_) : s(v_.s), v(v_.v) {}
213
214 Vector(const Vector<T, Vdense>& v_) : s(v_.size()), v(&v_[0]) {}
215
216 Vector(const Vector<T, Vdense_ref>& v_) : s(v_.s), v(v_.v) {}
217
218 explicit Vector(const Matrix<T, Mrectangle_constref>& m_) : s(m_.s1 * m_.s2), v(m_.m) {}
219
220 bool is_null() const { return v == 0; }
221
222 size_t size() const { return s; }
223
224 const T& operator[](size_t i) const { return v[i]; }
225
226 T norm2() const { return blas::nrm2(s, v, 1); }
227
228 T norm_inf() const {
229 T res = 0;
230 for (size_t i = 0; i != s; ++i) { res = max_abs(res, v[i]); }
231 return res;
232 }
233
234 Vector<T, Vdense_constref> operator()(const Interval& i) const {
235 return Vector<T, Vdense_constref>(i.size(), &v[i[0]]);
236 }
237
238 Vector<T, Vincrement_constref> operator()(const Slice& i) const {
239 return Vector<T, Vincrement_constref>(i.size(), &v[i[0]], i.increment());
240 }
241
242 Vector<T, Vindex_constref> operator()(const Index& i) const {
243 return Vector<T, Vindex_constref>(i.size(), &v[0], &i[0]);
244 }
245
246 static Vector null;
247
248 friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
249 l << "vector ";
250 l.write(v.s, v.v, 1);
251 return l;
252 }
253
254private:
255 size_t s;
256 const T* v;
257 MVFRIEND;
258};
259
260template <typename T>
261Vector<T, Vdense_constref> Vector<T, Vdense_constref>::null;
262
263struct Vdense_ref {
264 typedef Vincrement_ref base;
265 typedef Vincrement_scale_constref const_base;
266 typedef Vdense copy;
267};
268
269template <typename T>
270class Vector<T, Vdense_ref> {
271public:
272 Vector() : s(0), v(0) {}
273
274 Vector(size_t s_, T* v_) : s(s_), v(v_) {}
275
276 Vector(const Vector& v_) : s(v_.s), v(v_.v) {}
277
278 Vector(Vector<T, Vdense>& v_) : s(v_.size()), v(&v_[0]) {}
279
280 explicit Vector(const Matrix<T, Mrectangle_ref>& m_) : s(m_.s1 * m_.s2), v(m_.m) {}
281
282 Vector& operator=(const Vector& v_) {
283 if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
284 std::copy(v_.v, v_.v + s, v);
285 return *this;
286 }
287
288 bool is_null() const { return v == 0; }
289
290 size_t size() const { return s; }
291
292 void set_zero() { std::fill_n(v, s, 0); }
293
294 template <typename T1, typename STORAGE1>
295 Vector& operator=(const Vector<T1, STORAGE1>& v_) {
296 if (v_.size() != this->size()) {
297 Exception() << "The two vectors do not have the same size." << THROW;
298 }
299 Vector<T, Vincrement_ref>(*this) = v_;
300 return *this;
301 }
302
303 template <typename STORAGE>
304 Vector& operator+=(const Vector<T, STORAGE>& v_) {
305 *this = *this + v_;
306 return *this;
307 }
308
309 template <typename STORAGE>
310 Vector& operator-=(const Vector<T, STORAGE>& v_) {
311 *this = *this + (-v_);
312 return *this;
313 }
314
315 Vector& operator*=(T scale_) {
316 *this = *this * scale_;
317 return *this;
318 }
319
320 template <typename STORAGE>
321 void scale(const Vector<T, STORAGE>& v_) {
322 if (v_.size() != this->size()) {
323 Exception() << "The two vectors do not have the same size." << THROW;
324 }
325 for (size_t i = 0; i != this->size(); ++i) { (*this)[i] *= v_[i]; }
326 }
327
328 template <typename STORAGE>
329 void scale_inv(const Vector<T, STORAGE>& v_) {
330 if (v_.size() != this->size()) {
331 Exception() << "The two vectors do not have the same size." << THROW;
332 }
333 for (size_t i = 0; i != this->size(); ++i) { (*this)[i] /= v_[i]; }
334 }
335 const T& operator[](size_t i) const { return v[i]; }
336
337 T& operator[](size_t i) { return v[i]; }
338
339 Vector<T, Vdense_ref> operator()(const Interval& i) {
340 return Vector<T, Vdense_ref>(i.size(), &v[i[0]]);
341 }
342
343 Vector<T, Vincrement_ref> operator()(const Slice& i) {
344 return Vector<T, Vincrement_ref>(i.size(), &v[i[0]], i.increment());
345 }
346
347 Vector<T, Vindex_ref> operator()(const Index& i) {
348 return Vector<T, Vindex_ref>(i.size(), &v[0], &i[0]);
349 }
350
351 T norm2() const { return blas::nrm2(s, v, 1); }
352
353 T norm_inf() const {
354 double res = 0;
355 for (size_t i = 0; i != s; ++i) { res = max_abs(res, v[i]); }
356 return res;
357 }
358
359 static Vector null;
360
361 friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
362 l << "vector ";
363 l.write(v.s, v.v, 1);
364 return l;
365 }
366
367private:
368 size_t s;
369 T* v;
370 MVFRIEND;
371};
372
373template <>
374inline b2000::csda<double> Vector<b2000::csda<double>, Vdense>::norm2() const {
375 b2000::csda<double> n = 0;
376 for (auto const& i : *this) { n += std::norm(i); }
377 return n;
378}
379
380template <>
381inline b2000::csda<double> Vector<b2000::csda<double>, Vdense_ref>::norm2() const {
382 b2000::csda<double> n = 0;
383 for (size_t i = 0; i != size(); ++i) { n += std::norm(v[i]); }
384 return n;
385}
386
387template <>
388inline b2000::csda<double> Vector<b2000::csda<double>, Vdense_constref>::norm2() const {
389 b2000::csda<double> n = 0;
390 for (size_t i = 0; i != size(); ++i) { n += std::norm(v[i]); }
391 return n;
392}
393
394template <typename T>
395Vector<T, Vdense_ref> Vector<T, Vdense_ref>::null;
396
397}} // namespace b2000::b2linalg
398#endif
THROW
#define THROW
Definition b2exception.H:198
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32