04c1333727b9190df75e4f2093dac3c2f7e54c59
1 /*
2 * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;
3 * "views" mimic the semantic of C++ references: any operation performed
4 * on a "view" is actually performed on the "viewed object"
5 *
6 * Authors:
7 * Marco Cecchetti <mrcekets at gmail.com>
8 *
9 * Copyright 2008 authors
10 *
11 * This library is free software; you can redistribute it and/or
12 * modify it either under the terms of the GNU Lesser General Public
13 * License version 2.1 as published by the Free Software Foundation
14 * (the "LGPL") or, at your option, under the terms of the Mozilla
15 * Public License Version 1.1 (the "MPL"). If you do not alter this
16 * notice, a recipient may use your version of this file under either
17 * the MPL or the LGPL.
18 *
19 * You should have received a copy of the LGPL along with this library
20 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 * You should have received a copy of the MPL along with this library
23 * in the file COPYING-MPL-1.1
24 *
25 * The contents of this file are subject to the Mozilla Public License
26 * Version 1.1 (the "License"); you may not use this file except in
27 * compliance with the License. You may obtain a copy of the License at
28 * http://www.mozilla.org/MPL/
29 *
30 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
31 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
32 * the specific language governing rights and limitations.
33 */
38 #ifndef _NL_VECTOR_H_
39 #define _NL_VECTOR_H_
41 #include <cassert>
42 #include <algorithm> // for std::swap
43 #include <vector>
44 #include <sstream>
45 #include <string>
48 #include <gsl/gsl_vector.h>
49 #include <gsl/gsl_blas.h>
52 namespace Geom { namespace NL {
54 namespace detail
55 {
57 class BaseVectorImpl
58 {
59 public:
60 double const& operator[](size_t i) const
61 {
62 return *gsl_vector_const_ptr(m_vector, i);
63 }
65 const gsl_vector* get_gsl_vector() const
66 {
67 return m_vector;
68 }
69 bool is_zero() const
70 {
71 return gsl_vector_isnull(m_vector);
72 }
74 bool is_positive() const
75 {
76 for ( size_t i = 0; i < size(); ++i )
77 {
78 if ( (*this)[i] <= 0 ) return false;
79 }
80 return true;
81 }
83 bool is_negative() const
84 {
85 for ( size_t i = 0; i < size(); ++i )
86 {
87 if ( (*this)[i] >= 0 ) return false;
88 }
89 return true;
90 }
92 bool is_non_negative() const
93 {
94 for ( size_t i = 0; i < size(); ++i )
95 {
96 if ( (*this)[i] < 0 ) return false;
97 }
98 return true;
99 }
101 double max() const
102 {
103 return gsl_vector_max(m_vector);
104 }
106 double min() const
107 {
108 return gsl_vector_min(m_vector);
109 }
111 size_t max_index() const
112 {
113 return gsl_vector_max_index(m_vector);
114 }
116 size_t min_index() const
117 {
118 return gsl_vector_min_index(m_vector);
119 }
121 size_t size() const
122 {
123 return m_size;
124 }
126 std::string str() const;
128 virtual ~BaseVectorImpl()
129 {
130 }
132 protected:
133 size_t m_size;
134 gsl_vector* m_vector;
136 }; // end class BaseVectorImpl
139 inline
140 bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)
141 {
142 if (v1.size() != v2.size()) return false;
144 for (size_t i = 0; i < v1.size(); ++i)
145 {
146 if (v1[i] != v2[i]) return false;
147 }
148 return true;
149 }
151 template< class charT >
152 inline
153 std::basic_ostream<charT> &
154 operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)
155 {
156 if (_vector.size() == 0 ) return os;
157 os << "[" << _vector[0];
158 for (unsigned int i = 1; i < _vector.size(); ++i)
159 {
160 os << ", " << _vector[i];
161 }
162 os << "]";
163 return os;
164 }
166 inline
167 std::string BaseVectorImpl::str() const
168 {
169 std::ostringstream oss;
170 oss << (*this);
171 return oss.str();
172 }
174 inline
175 double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)
176 {
177 double result;
178 gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);
179 return result;
180 }
183 class VectorImpl : public BaseVectorImpl
184 {
185 public:
186 typedef BaseVectorImpl base_type;
188 public:
189 void set_all(double x)
190 {
191 gsl_vector_set_all(m_vector, x);
192 }
194 void set_basis(size_t i)
195 {
196 gsl_vector_set_basis(m_vector, i);
197 }
199 using base_type::operator[];
201 double & operator[](size_t i)
202 {
203 return *gsl_vector_ptr(m_vector, i);
204 }
206 using base_type::get_gsl_vector;
208 gsl_vector* get_gsl_vector()
209 {
210 return m_vector;
211 }
213 void swap_elements(size_t i, size_t j)
214 {
215 gsl_vector_swap_elements(m_vector, i, j);
216 }
218 void reverse()
219 {
220 gsl_vector_reverse(m_vector);
221 }
223 VectorImpl & scale(double x)
224 {
225 gsl_vector_scale(m_vector, x);
226 return (*this);
227 }
229 VectorImpl & translate(double x)
230 {
231 gsl_vector_add_constant(m_vector, x);
232 return (*this);
233 }
235 VectorImpl & operator+=(base_type const& _vector)
236 {
237 gsl_vector_add(m_vector, _vector.get_gsl_vector());
238 return (*this);
239 }
241 VectorImpl & operator-=(base_type const& _vector)
242 {
243 gsl_vector_sub(m_vector, _vector.get_gsl_vector());
244 return (*this);
245 }
247 }; // end class VectorImpl
249 } // end namespace detail
252 using detail::operator==;
253 using detail::operator<<;
255 class Vector : public detail::VectorImpl
256 {
257 public:
258 typedef detail::VectorImpl base_type;
260 public:
261 Vector(size_t n)
262 {
263 m_size = n;
264 m_vector = gsl_vector_alloc(n);
265 }
267 Vector(size_t n, double x)
268 {
269 m_size = n;
270 m_vector = gsl_vector_alloc(n);
271 gsl_vector_set_all(m_vector, x);
272 }
274 // create a vector with n elements all set to zero
275 // but the i-th that is set to 1
276 Vector(size_t n, size_t i)
277 {
278 m_size = n;
279 m_vector = gsl_vector_alloc(n);
280 gsl_vector_set_basis(m_vector, i);
281 }
283 Vector(Vector const& _vector)
284 : base_type()
285 {
286 m_size = _vector.size();
287 m_vector = gsl_vector_alloc(size());
288 gsl_vector_memcpy(m_vector, _vector.m_vector);
289 }
291 explicit
292 Vector(base_type::base_type const& _vector)
293 {
294 m_size = _vector.size();
295 m_vector = gsl_vector_alloc(size());
296 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
297 }
299 virtual ~Vector()
300 {
301 gsl_vector_free(m_vector);
302 }
305 Vector & operator=(Vector const& _vector)
306 {
307 assert( size() == _vector.size() );
308 gsl_vector_memcpy(m_vector, _vector.m_vector);
309 return (*this);
310 }
312 Vector & operator=(base_type::base_type const& _vector)
313 {
314 assert( size() == _vector.size() );
315 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
316 return (*this);
317 }
319 Vector & scale(double x)
320 {
321 return static_cast<Vector&>( base_type::scale(x) );
322 }
324 Vector & translate(double x)
325 {
326 return static_cast<Vector&>( base_type::translate(x) );
327 }
329 Vector & operator+=(base_type::base_type const& _vector)
330 {
331 return static_cast<Vector&>( base_type::operator+=(_vector) );
332 }
334 Vector & operator-=(base_type::base_type const& _vector)
335 {
336 return static_cast<Vector&>( base_type::operator-=(_vector) );
337 }
339 friend
340 void swap(Vector & v1, Vector & v2);
341 friend
342 void swap_any(Vector & v1, Vector & v2);
344 }; // end class Vector
347 // warning! these operations invalidate any view of the passed vector objects
348 inline
349 void swap(Vector & v1, Vector & v2)
350 {
351 assert( v1.size() == v2.size() );
352 std::swap(v1.m_vector, v2.m_vector);
353 }
355 inline
356 void swap_any(Vector & v1, Vector & v2)
357 {
358 std::swap(v1.m_vector, v2.m_vector);
359 std::swap(v1.m_size, v2.m_size);
360 }
363 class ConstVectorView : public detail::BaseVectorImpl
364 {
365 public:
366 typedef detail::BaseVectorImpl base_type;
368 public:
369 ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)
370 : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )
371 {
372 m_size = n;
373 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
374 }
376 ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)
377 : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )
378 {
379 m_size = n;
380 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
381 }
383 ConstVectorView(const double* _vector, size_t n, size_t offset = 0)
384 : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )
385 {
386 m_size = n;
387 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
388 }
390 ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)
391 : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )
392 {
393 m_size = n;
394 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
395 }
397 explicit
398 ConstVectorView(gsl_vector_const_view _gsl_vector_view)
399 : m_vector_view(_gsl_vector_view)
400 {
401 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
402 m_size = m_vector->size;
403 }
405 explicit
406 ConstVectorView(const std::vector<double>& _vector)
407 : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )
408 {
409 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
410 m_size = _vector.size();
411 }
413 ConstVectorView(const ConstVectorView & _vector)
414 : base_type(),
415 m_vector_view(_vector.m_vector_view)
416 {
417 m_size = _vector.size();
418 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
419 }
421 ConstVectorView(const base_type & _vector)
422 : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))
423 {
424 m_size = _vector.size();
425 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
426 }
428 private:
429 gsl_vector_const_view m_vector_view;
431 }; // end class ConstVectorView
436 class VectorView : public detail::VectorImpl
437 {
438 public:
439 typedef detail::VectorImpl base_type;
441 public:
442 VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)
443 {
444 m_size = n;
445 if (stride == 1)
446 {
447 m_vector_view
448 = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);
449 m_vector = &(m_vector_view.vector);
450 }
451 else
452 {
453 m_vector_view
454 = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);
455 m_vector = &(m_vector_view.vector);
456 }
457 }
459 VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)
460 {
461 m_size = n;
462 if (stride == 1)
463 {
464 m_vector_view
465 = gsl_vector_view_array(_vector + offset, n);
466 m_vector = &(m_vector_view.vector);
467 }
468 else
469 {
470 m_vector_view
471 = gsl_vector_view_array_with_stride(_vector + offset, stride, n);
472 m_vector = &(m_vector_view.vector);
473 }
475 }
477 VectorView(const VectorView & _vector)
478 : base_type()
479 {
480 m_size = _vector.size();
481 m_vector_view = _vector.m_vector_view;
482 m_vector = &(m_vector_view.vector);
483 }
485 VectorView(Vector & _vector)
486 {
487 m_size = _vector.size();
488 m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());
489 m_vector = &(m_vector_view.vector);
490 }
492 explicit
493 VectorView(gsl_vector_view _gsl_vector_view)
494 : m_vector_view(_gsl_vector_view)
495 {
496 m_vector = &(m_vector_view.vector);
497 m_size = m_vector->size;
498 }
500 explicit
501 VectorView(std::vector<double> & _vector)
502 {
503 m_size = _vector.size();
504 m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());
505 m_vector = &(m_vector_view.vector);
506 }
508 VectorView & operator=(VectorView const& _vector)
509 {
510 assert( size() == _vector.size() );
511 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
512 return (*this);
513 }
515 VectorView & operator=(base_type::base_type const& _vector)
516 {
517 assert( size() == _vector.size() );
518 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
519 return (*this);
520 }
522 VectorView & scale(double x)
523 {
524 return static_cast<VectorView&>( base_type::scale(x) );
525 }
527 VectorView & translate(double x)
528 {
529 return static_cast<VectorView&>( base_type::translate(x) );
530 }
532 VectorView & operator+=(base_type::base_type const& _vector)
533 {
534 return static_cast<VectorView&>( base_type::operator+=(_vector) );
535 }
537 VectorView & operator-=(base_type::base_type const& _vector)
538 {
539 return static_cast<VectorView&>( base_type::operator-=(_vector) );
540 }
542 friend
543 void swap_view(VectorView & v1, VectorView & v2);
545 private:
546 gsl_vector_view m_vector_view;
548 }; // end class VectorView
551 inline
552 void swap_view(VectorView & v1, VectorView & v2)
553 {
554 assert( v1.size() == v2.size() );
555 std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
556 }
559 } } // end namespaces
562 #endif /*_NL_VECTOR_H_*/
564 /*
565 Local Variables:
566 mode:c++
567 c-file-style:"stroustrup"
568 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
569 indent-tabs-mode:nil
570 fill-column:99
571 End:
572 */
573 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :