source: trunk/yat/utility/VectorBase.h

Last change on this file was 3871, checked in by Peter, 17 months ago

closes #932

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
  • Property svndigest:ignore set to 1008
File size: 9.6 KB
Line 
1#ifndef _theplu_yat_utility_vector_base_
2#define _theplu_yat_utility_vector_base_
3
4// $Id: VectorBase.h 3871 2020-03-01 22:46:57Z peter $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
10  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
11  Copyright (C) 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
12  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
13  Copyright (C) 2009, 2012, 2017, 2020 Peter Johansson
14
15  This file is part of the yat library, http://dev.thep.lu.se/trac/yat
16
17  The yat library is free software; you can redistribute it and/or
18  modify it under the terms of the GNU General Public License as
19  published by the Free Software Foundation; either version 3 of the
20  License, or (at your option) any later version.
21
22  The yat library is distributed in the hope that it will be useful,
23  but WITHOUT ANY WARRANTY; without even the implied warranty of
24  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25  General Public License for more details.
26
27  You should have received a copy of the GNU General Public License
28  along with yat. If not, see <http://www.gnu.org/licenses/>.
29*/
30
31#include "BasicVector.h"
32#include "BLAS_level1.h"// include these as they are expected when using Vector
33#include "BLAS_level2.h"// include these as they are expected when using Vector
34#include "StrideIterator.h"
35#include "VectorExpression.h"
36
37#include <iosfwd>
38#include <vector>
39#include <cstddef> // size_t
40
41#include <gsl/gsl_vector.h>
42
43namespace theplu {
44namespace yat {
45namespace utility {
46
47  class Vector;
48
49  /**
50     @brief This is the yat interface to GSL vector.
51
52     This is an interface class for vectors containing the const
53     interface. For mutable functionality see VectorMutable.
54  */
55  class VectorBase : public BasicVector<VectorBase>
56  {
57  public:
58    /**
59       value_type is double
60
61       \since New in yat 0.5
62     */
63    typedef double value_type;
64
65    /**
66       const_reference type is const double&
67
68       \since New in yat 0.5
69     */
70    typedef const double& const_reference;
71
72    /// \brief VectorBase::const_iterator
73    typedef StrideIterator<const double*> const_iterator;
74
75    /**
76       \brief Constructor.
77    */
78    VectorBase(const gsl_vector* v=NULL);
79
80    ///
81    /// The destructor.
82    ///
83    virtual ~VectorBase(void);
84
85    /**
86       \return read-only iterator to start of VectorBase
87     */
88    const_iterator begin(void) const;
89
90    /**
91       \return read-only iterator to end of VectorBase
92     */
93    const_iterator end(void) const;
94
95    /**
96       \brief Check whether VectorBases are equal within a user defined
97       precision, set by \a precision.
98
99       \return True if each element deviates less or equal than \a
100       d. If any VectorBase contain a NaN, false is always returned.
101
102       \see operator== and operator!=
103    */
104    bool equal(const VectorBase&, const double precision=0) const;
105
106    /**
107       \return A const pointer to the internal GSL vector,
108    */
109    const gsl_vector* gsl_vector_p(void) const;
110
111    /**
112       Check if the vector object is a view (sub-vector) to another
113       vector.
114
115       \return True if the object is a view, false othwerwise.
116     */
117    virtual bool isview(void) const=0;
118
119    /**
120       \return number of elements in the VectorBase.
121    */
122    size_t size(void) const;
123
124    /**
125       \brief Element access operator.
126
127       \return Const reference to element \a i.
128
129       \throw If GSL range checks are enabled in the underlying GSL
130       library a GSL_error exception is thrown if either index is out
131       of range.
132    */
133    const double& operator()(size_t i) const;
134
135    /**
136       \brief Comparison operator. Takes linear time.
137
138       Checks are performed with exact matching, i.e., rounding off
139       effects may destroy comparison. Use the equal function for
140       comparing elements within a user defined precision.
141
142       \return True if all elements are equal otherwise false.
143
144       \see equal(const VectorBase&, const double precision=0)
145    */
146    bool operator==(const VectorBase&) const;
147
148    /**
149       \brief Comparison operator. Takes linear time.
150
151       Checks are performed with exact matching, i.e., rounding off
152       effects may destroy comparison. Use the equal function for
153       comparing elements within a user defined precision.
154
155       \return False if all elements are equal otherwise true.
156
157       \see equal(const VectorBase&, const double precision=0)
158    */
159    bool operator!=(const VectorBase&) const;
160
161    ///
162    /// @return The dot product.
163    ///
164    double operator*(const VectorBase&) const;
165
166  protected:
167    /// pointer to underlying GSL vector
168    const gsl_vector* const_vec_;
169
170  private:
171    // copy assignment not allowed
172    const VectorBase& operator=(const VectorBase&);
173  };
174
175  /**
176     \brief Check if all elements of the VectorBase are zero.
177
178     \return True if all elements in the VectorBase is zero, false
179     othwerwise.
180
181     \relates VectorBase
182  */
183  bool isnull(const VectorBase&);
184
185  /**
186     \brief Get the maximum value of the VectorBase.
187
188     \return The maximum value of the VectorBase.
189
190     \relates VectorBase
191  */
192  double max(const VectorBase&);
193
194  /**
195     \brief Locate the maximum value in the VectorBase.
196
197     \return The index to the maximum value of the VectorBase.
198
199     \note Lower index has precedence.
200
201     \relates VectorBase
202  */
203  size_t max_index(const VectorBase&);
204
205  /**
206     \brief Get the minimum value of the VectorBase.
207
208     \return The minimum value of the VectorBase.
209
210     \relates VectorBase
211  */
212  double min(const VectorBase&);
213
214  /**
215     \brief Locate the minimum value in the VectorBase.
216
217     \return The index to the minimum value of the VectorBase.
218
219     \note Lower index has precedence.
220
221     \relates VectorBase
222  */
223  size_t min_index(const VectorBase&);
224
225  /**
226     \brief Create a VectorBase \a flag indicating NaN's in another VectorBase
227     \a templat.
228
229     The \a flag VectorBase is changed to contain 1's and 0's only. A 1
230     means that the corresponding element in the \a templat VectorBase is
231     valid and a zero means that the corresponding element is a NaN.
232
233     \note Space for vector \a flag is reallocated to fit the size of
234     VectorBase \a templat if sizes mismatch.
235
236     \return True if the \a templat VectorBase contains at least one NaN.
237
238     \relates VectorBase
239  */
240  bool nan(const VectorBase& templat, Vector& flag);
241
242  /**
243     \return L2 norm squared, or the dot product of \a v with itself
244
245     \relates VectorBase
246
247     \since New in yat 0.18
248   */
249  double norm2_squared(const VectorBase& v);
250
251
252  /**
253     Create a vector \a sort_index containing the indeces of
254     elements in a another VectorBase \a invec.  The elements of \a
255     sort_index give the index of the VectorBase element which would
256     have been stored in that position if the VectorBase had been sorted
257     in place. The first element of \a sort_index gives the index of the least
258     element in \a invec, and the last element of \a sort_index gives the
259     index of the greatest element in \a invec . The VectorBase \a invec
260     is not changed.
261
262     \relatesalso VectorBase
263  */
264  void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec);
265
266  /**
267      Similar to sort_index but creates a VectorBase with indices to
268      the \a k smallest elements in \a invec.
269
270     \relatesalso VectorBase
271  */
272  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
273                           const VectorBase& invec);
274
275  /**
276      Similar to sort_index but creates a VectorBase with indices to
277      the \a k largest elements in \a invec.
278
279     \relatesalso VectorBase
280  */
281  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
282                          const VectorBase& invec);
283
284  /**
285     \brief Calculate the sum of all VectorBase elements.
286
287     \return The sum.
288
289     \relates VectorBase
290  */
291  double sum(const VectorBase&);
292
293  /**
294     \brief The output operator for the VectorBase class.
295
296     Elements in VectorBase \a v are sent to ostream \a s and
297     separated with the fill character of stream \a s, s.fill(). If
298     you, for example, want to print the VectorBase \a v with the
299     elements separated by a ':', you do so by:
300     \verbatim
301     s << setfill(':') << v;
302     \endverbatim
303
304     \relates VectorBase
305  */
306  std::ostream& operator<<(std::ostream& s, const VectorBase& v);
307
308  /// \cond IGNORE_DOXYGEN
309
310  // Some convenience functions used in Vector and friends.
311  namespace detail {
312
313    /**
314       \brief Create a new gsl vector
315
316       Necessary memory for the new GSL vector is allocated and the
317       caller is responsible for freeing the allocated memory.
318
319       \return A pointer to created GSL vector.
320
321       \throw GSL_error if memory cannot be allocated for the new
322       GSL vector.
323    */
324    gsl_vector* create_gsl_vector(size_t n);
325
326    /**
327       same as above but sets all values to init
328     */
329    gsl_vector* create_gsl_vector(size_t n, double init);
330
331    /**
332       \brief Create a new copy of the internal GSL vector.
333
334       Necessary memory for the new GSL vector is allocated and the
335       caller is responsible for freeing the allocated memory.
336
337       \return A pointer to a copy of the internal GSL vector.
338
339       \throw GSL_error if memory cannot be allocated for the new
340       copy, or if dimensions mis-match.
341    */
342    gsl_vector* create_gsl_vector_copy(const gsl_vector*);
343
344
345    // return true if a and b overlap in memory
346    bool overlap(const gsl_vector* a, const gsl_vector* b);
347
348    /*
349      Is a stricter condition that overlap(2). Basically it returns
350      false if it is safe to write code like:
351      for (size_t i=0; i<a->size; ++i)
352        a[i] = foo(b[i])
353
354      in other words, it returns true if there exists i and j, i<j, such that
355      gsl_vector_const_ptr(a, j) == gsl_vector_const_ptr(b, i)
356
357      For speed reasons function might return true in some cases when
358      condition above is not true, when both a and b have stride>1 and
359      elements sit between each other.
360
361      Note, a and b must have same size.
362     */
363    bool serial_overlap(const gsl_vector* a, const gsl_vector* b);
364  }
365
366  /// \endcond
367
368}}} // of namespace utility, yat, and theplu
369
370#endif
Note: See TracBrowser for help on using the repository browser.