// Copyright (c) 1997-2000 // Utrecht University (The Netherlands), // ETH Zurich (Switzerland), // INRIA Sophia-Antipolis (France), // Max-Planck-Institute Saarbruecken (Germany), // and Tel-Aviv University (Israel). All rights reserved. // // This file is part of CGAL (www.cgal.org) // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Kernel_d/include/CGAL/Kernel_d/Vector__.h $ // $Id: Vector__.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Michael Seel #ifndef CGAL_VECTOR___H #define CGAL_VECTOR___H #include #include #include #include #include #include #include #include #include #include namespace CGAL { namespace Linear_Algebra { template class Vector_; template class Matrix_; /*{\Msubst # # Vector_#Vector Matrix_#Matrix }*/ /*{\Moptions print_title=yes}*/ /*{\Moptions outfile=Vector.man}*/ /*{\Xtext \headerline{Common Notation} The following data types use the concept of iterator ranges as an abstraction of tuples and sets. For an iterator range |[first,last)| we define |S = set [first,last)| as the ordered tuple $(|S[0]|,|S[1]|, \ldots |S[d-1]|)$ where $|S[i]| = |*| |++|^{(i)}|first|$ (the element obtained by forwarding the iterator by operator |++| $i$ times and then dereferencing it to get the value to which it points). We write |d = size [first,last)|. This extends the syntax of random access iterators to input iterators. If we index the tuple as above then we require that $|++|^{(d)}|first == last|$ (note that |last| points beyond the last element to be accepted).}*/ /*{\Manpage {Vector}{}{Vectors with NT Entries}{v}}*/ template class Vector_ { /*{\Mdefinition An instance of data type |Vector_| is a vector of variables of number type |NT|. Together with the type |Matrix_| it realizes the basic operations of linear algebra.}*/ public: /*{\Mtypes 5.5}*/ typedef NT_* pointer; typedef const NT_* const_pointer; typedef NT_ NT; /*{\Mtypemember the ring type of the components.}*/ typedef pointer iterator; /*{\Mtypemember the iterator type for accessing components.}*/ typedef const_pointer const_iterator; /*{\Mtypemember the const iterator type for accessing components.}*/ typedef AL_ allocator_type; /*{\Xtypemember the allocator type.}*/ protected: friend class Matrix_; NT* v_; int d_; allocator_type& allocator() { CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(allocator_type, MM); return MM; } inline void allocate_vec_space(NT*& vi, int di) { /* We use this procedure to allocate memory. We first get an appropriate piece of memory from the allocator and then initialize each cell by an inplace new. */ vi = allocator().allocate(di); NT* p = vi + di - 1; while (p >= vi) { new (p) NT(0); p--; } } inline void deallocate_vec_space(NT*& vi, int di) { /* We use this procedure to deallocate memory. We have to free it by the allocator scheme. We first call the destructor for type NT for each cell of the array and then return the piece of memory to the memory manager. */ NT* p = vi + di - 1; while (p >= vi) { std::allocator_traits::destroy(allocator(),p); p--; } //af: as proposed by sylvain allocator().deallocate(vi, di); vi = (NT*)0; } inline void check_dimensions(const Vector_& vec) const { CGAL_USE(vec); CGAL_assertion_msg((d_ == vec.d_), "Vector_::check_dimensions: object dimensions disagree."); } public: /*{\Mcreation v 3}*/ Vector_() : v_(0),d_(0) {} /*{\Mcreate creates an instance |\Mvar| of type |\Mname|.}*/ Vector_(int d) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to a vector of dimension $d$.}*/ { CGAL_assertion_msg( d >= 0 , "Vector_::constructor: negative dimension."); d_ = d; v_ = (NT*)0; if (d_ > 0){ allocate_vec_space(v_,d_); while (d--) v_[d] = NT(0); } } Vector_(int d, const NT& x) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to a vector of dimension $d$ with entries |x|.}*/ { CGAL_assertion_msg( d >= 0 , "Vector_::constructor: negative dimension."); d_ = d; v_ = (NT*)0; if (d_ > 0){ allocate_vec_space(v_,d_); while (d--) v_[d] = x; } } template Vector_(Forward_iterator first, Forward_iterator last) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|; |\Mvar| is initialized to the vector with entries |set [first,last)|. \require |Forward_iterator| has value type |NT|.}*/ { #if defined _MSC_VER && _MSC_VER == 1300 d_ = 0; Forward_iterator fit = first; while(fit++!=last) d_++; #else d_ = static_cast(std::distance(first, last)); #endif allocate_vec_space(v_,d_); iterator it = begin(); while (first != last) { *it = *first; ++it; ++first; } } Vector_(const Vector_& p) { d_ = p.d_; if (d_ > 0) allocate_vec_space(v_,d_); else v_ = (NT*)0; for(int i=0; i& operator=(const Vector_& vec) { if (&vec == this) return *this; int n = vec.d_; if (n != d_) { if (d_ > 0) deallocate_vec_space(v_,d_); d_=n; if (n > 0) allocate_vec_space(v_,n); else v_ = (NT*)0; } while (n--) v_[n] = vec.v_[n]; return *this; } ~Vector_() { if (d_ > 0) deallocate_vec_space(v_,d_); } /*{\Moperations 3 4}*/ int dimension() const { return d_; } /*{\Mop returns the dimension of |\Mvar|.}*/ bool is_zero() const /*{\Mop returns true iff |\Mvar| is the zero vector.}*/ { for(int i=0; i operator+(const Vector_& v1) const; /*{\Mbinop Addition. \precond\\ |v.dimension() == v1.dimension()|.}*/ Vector_ operator-(const Vector_& v1) const; /*{\Mbinop Subtraction. \precond\\ |v.dimension() = v1.dimension()|.}*/ NT operator*(const Vector_& v1) const; /*{\Mbinop Inner Product. \precond\\ |v.dimension() = v1.dimension()|.}*/ Vector_ compmul(const NT& r) const; Vector_ operator-() const; /*{\Munop Negation.}*/ Vector_& operator+=(const Vector_& v1); /*{\Mbinop Addition plus assignment. \precond\\ |v.dimension() == v1.dimension()|.}*/ Vector_& operator-=(const Vector_& v1); /*{\Mbinop Subtraction plus assignment. \precond\\ |v.dimension() == v1.dimension()|.}*/ Vector_& operator*=(const NT& s); /*{\Mbinop Scalar multiplication plus assignment.}*/ Vector_& operator/=(const NT& s); /*{\Mbinop Scalar division plus assignment.}*/ bool operator==(const Vector_& w) const; bool operator!=(const Vector_& w) const { return !(*this == w); } static int compare(const Vector_&, const Vector_&); }; template inline Vector_ operator*(const NT& r, const Vector_& v) /*{\Mbinopfunc Componentwise multiplication with number $r$.}*/ { return v.compmul(r); } template inline Vector_ operator*(const Vector_& v, const NT& r) /*{\Mbinopfunc Componentwise multiplication with number $r$.}*/ { return v.compmul(r); } template inline Vector_& Vector_:: operator+=(const Vector_& vec) { check_dimensions(vec); int n = d_; while (n--) v_[n] += vec.v_[n]; return *this; } template inline Vector_& Vector_:: operator-=(const Vector_& vec) { check_dimensions(vec); int n = d_; while (n--) v_[n] -= vec.v_[n]; return *this; } template inline Vector_& Vector_:: operator*=(const NT& s) { int n = d_; while (n--) v_[n] *= s; return *this; } template inline Vector_& Vector_:: operator/=(const NT& s) { int n = d_; while (n--) v_[n] /= s; return *this; } template inline Vector_ Vector_:: operator+(const Vector_& vec) const { check_dimensions(vec); int n = d_; Vector_ result(n); while (n--) result.v_[n] = v_[n]+vec.v_[n]; return result; } template inline Vector_ Vector_:: operator-(const Vector_& vec) const { check_dimensions(vec); int n = d_; Vector_ result(n); while (n--) result.v_[n] = v_[n]-vec.v_[n]; return result; } template inline Vector_ Vector_:: operator-() const // unary minus { int n = d_; Vector_ result(n); while (n--) result.v_[n] = -v_[n]; return result; } template inline Vector_ Vector_:: compmul(const NT& x) const { int n = d_; Vector_ result(n); while (n--) result.v_[n] = v_[n] * x; return result; } template inline NT_ Vector_:: operator*(const Vector_& vec) const { check_dimensions(vec); NT_ result=0; int n = d_; while (n--) result = result+v_[n]*vec.v_[n]; return result; } template inline bool Vector_:: operator==(const Vector_& vec) const { if (vec.d_ != d_) return false; int i = 0; while ((i int Vector_:: compare(const Vector_& v1, const Vector_& v2) { int i; v1.check_dimensions(v2); for(i=0; i < v1.dimension() && v1[i]==v2[i]; i++) {} if (i == v1.dimension()) return 0; return (v1[i] < v2[i]) ? -1 : 1; } template std::ostream& operator<<(std::ostream& os, const Vector_& v) /*{\Xbinopfunc writes |\Mvar| componentwise to the output stream $O$.}*/ { /* syntax: d x_0 x_1 ... x_d-1 */ int d = v.dimension(); switch (get_mode(os)) { case CGAL::IO::BINARY: CGAL::write( os, d); for ( int i = 0; i < d; ++i) CGAL::write( os, v[i]); break; case CGAL::IO::ASCII: os << d; for ( int i = 0; i < d; ++i) os << ' ' << v[i]; break; case CGAL::IO::PRETTY: os << "LA::Vector(" << d << " ["; for ( int i = 0; i < d; ++i) { if ( i > 0) os << ',' << ' '; os << v[i]; } os << "])"; break; } return os; } template std::istream& operator>>(std::istream& is, Vector_& v) /*{\Xbinopfunc reads |\Mvar| componentwise from the input stream $I$.}*/ { /* syntax: d x_0 x_1 ... x_d-1 */ int d; switch (get_mode(is)) { case CGAL::IO::ASCII : case CGAL::IO::BINARY : is >> d; v = Vector_(d); for ( int i = 0; i < d; ++i) { is >> v[i]; } break; default: std::cerr<<"\nStream must be in ascii or binary mode"<