DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DenseLinAlgPack_DVectorClassTmpl.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef VECTOR_CLASS_TMPL_H
43 #define VECTOR_CLASS_TMPL_H
44 
45 #include <vector>
46 
47 #include "DenseLinAlgPack_Types.hpp"
48 #include "StrideIterPack_StrideIter.hpp"
49 
50 namespace DenseLinAlgPack{
51 
52 // ////////////////////////////////////////////////////////////////////////////////
53 /* * @name {\bf Dense 1-D DVector Abstractions}.
54  *
55  * These are classes that abstract 1-D vectors. The class \Ref{DVector} is a storage class
56  * for vectors while the class \Ref{DVectorSlice} is used to represent regions of vectors
57  * , for rows, columns or diagonals of matrices (see \Ref{DMatrix}
58  * and \Ref{DMatrixSlice}).
59  */
60 
61 // @{
62 
63 /* * @name {\bf DVector Classes}. */
64 // @{
65 
67 /* * Slice of a 1-D sequential C++ array treated as a vector.
68  *
69  * Objects of this class represent regions of vectors (continuous), rows of matrices
70  * , columns of matrices or diagonals of matrices.
71  * The underlying representation is of a continuous C++ array with non unit stride.
72  * It uses the same convention that the BLAS use where a vector is represented as
73  * the first element of
74  * in an array, the stride between elements in that array and the number of elements.
75  * Changes to elements through a DVectorSlice object result in changes to the elements
76  * in the underlying value_type* data.
77  *
78  * DVectorSlice provides many STL compliant features such as typedef type members
79  *, iterator returning functions
80  * and the dim() function. It also provides access to individual elements (lvalue)
81  * through 0-based
82  * and 1-based subscripting with operator[](i) and operator()(i) respectively.
83  * In addition and subregions can be created
84  * by subscripting (operator()()) with an \Ref{Range1D} object or using lower (>0)
85  * and upper bounds of the region.
86  */
87 template<class T>
89 public:
90 
91  /* * @name {\bf Nested Member Types (STL)}.
92  *
93  * These nested types give the types used in the interface to the class.
94  *
95  * \begin{description}
96  * <li>[#value_type#] - type being stored in the underlying C++ array
97  * <li>[#size_type#] - type used as an index and for the number of elements
98  * in the vector
99  * <li>[#difference_type#] - type for the distance between elements and the stride
100  * <li>[#iterator#] - type for the forward non-constant iterator
101  * <li>[#const_iterator#] - type for the forward constant iterator (can't change elements)
102  * <li>[#reverse_iterator#] - type for the reverse non-constant iterator
103  * <li>[#const_reverse_iterator#] - type for the reverse constant iterator (can't change elements)
104  * <li>[#reference#] - type returned from subscripting, iterator deferencing etc.
105  * <li>[#const_reference#] - "" "" for const vector slice objects
106  * \end{description}
107  */
108 
109  // @{
110  // @}
111 
112  typedef T value_type;
113  typedef DenseLinAlgPack::size_type size_type;
114  typedef ptrdiff_t difference_type;
115  typedef StrideIterPack::stride_iter<value_type*
116  , value_type, value_type&, value_type*
117  , difference_type> iterator;
118  typedef StrideIterPack::stride_iter<const value_type*
119  , value_type, const value_type&, const value_type*
120  , difference_type> const_iterator;
121 #if defined(_INTEL_CXX) || defined (_INTEL_CXX)
122  typedef std::reverse_iterator<iterator, value_type
123  , value_type&, value_type*, difference_type> reverse_iterator;
124  typedef std::reverse_iterator<const_iterator
125  , value_type, const value_type&, const value_type*
126  , difference_type> const_reverse_iterator;
127 #else
128  typedef std::reverse_iterator<iterator> reverse_iterator;
129  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
130 #endif
131  typedef value_type& reference;
132  typedef const value_type& const_reference;
133 
134  /* * @name {\bf Constructors}.
135  *
136  * The user usually does not need not call any of these constructors
137  * explicitly to create a vector slice.
138  * These
139  * constructors are used by the classes in the library to construct VectorSliceTmpl objects.
140  * Instead, users create VectorSliceTmpl objects by indexing (\Ref{Range1D}) a \Ref{DVector}
141  * , or \Ref{VectorSliceTmpl}
142  * object or calling row(...), col(...) or diag(...) on a \Ref{DMatrix} or
143  * \Ref{DMatrixSlice} object.
144  * The default C++ copy constructor is used, and is therefore not show here.
145  *
146  * Constructors are also included for creating views of raw C++ arrays.
147  * These constructors take non-#const# pointers. However you can savely
148  * create a #const# view of a #const# C++ array by using a constant cast.
149  * For example:
150  *
151  \verbatim
152  const VectorSliceTmpl<T>::size_type n = 5;
153  const VectorSliceTmpl<T>::value_type ptr[n] = { 1.0, 2.0, 3.0, 4.0, 5,0 };
154  const VectorSliceTmpl vec(const_cast<VectorSliceTmpl<T>::value_type*>(ptr),n);
155  \endverbatim
156  */
157 
158  // @{
159 
161  /* * Creates an empty view.
162  *
163  * You must use bind(...) to bind to a view to initialize after construction.
164  */
165  VectorSliceTmpl();
167  /* * Creates a VectorSice object that represents a non-continous region of a raw C++ array.
168  *
169  * Of course the sequence of elements #ptr[stride * i]# for #i# = 0, 1, ..., #size#-1
170  * must yield valid properly allocated regions of memory. It is up the the user to insure
171  * that they do.
172  *
173  * @param ptr Pointer to the first element in the raw C++ array
174  * @param size number of elements in the vector slice
175  * @param stride the distance (may be negative) between each successive element (default = 1)
176  */
177  VectorSliceTmpl( value_type* ptr, size_type size, difference_type stride = 1 );
179  /* * Creates a VectorSliceTmpl object that represents a continous region of a raw C++ array.
180  *
181  * The VectorSliceTmpl Object represents the following elements of raw array:
182  *
183  * #ptr[rng.lbound()-1+i]#, for #i# = 0, 1, ..., #rng.ubound()-1#
184  *
185  * Preconditions: <ul>
186  * <li> #rng.ubound() + 1 <= n# (throw std::out_of_range)
187  * </ul>
188  *
189  * @param ptr Pointer to the first element in the raw C++ array
190  * @param size number of elements in the vector slice
191  * @param rng index range (1-based) of the region being represented.
192  * Here rng.full_range() can be true.
193  */
194  VectorSliceTmpl( value_type* ptr, size_type size, const Range1D& rng );
196  /* * Create a VectorSliceTmpl that represents a continous region of the existing VectorSliceTmpl object, vs.
197  *
198  * The index, rng, is relative to the VectorSliceTmpl object, vs.
199  * For example rng = [1,3] would create a VectorSliceTmpl object
200  * representing the elements 2, 4 and 6. The following
201  * shows the elements represented by each of the participating objects.
202  \verbatim
203 
204  vs = [2, 4, 6, 8, 10]
205  this = [2, 4, 6]
206 
207  \endverbatim
208 
209  * Preconditions: <ul>
210  * <li> rng.full_range() == false (throw #std::out_of_range#)
211  * <li> rng.dim() <= vs.dim() (throw #std::out_of_range#)
212  * </ul>
213  *
214  * @param vs VectorSliceTmpl object that this VectorSliceTmpl object is being created from
215  * @param rng Range1D range of the vector slice being created.
216  */
218 
219  // @}
220 
223 
224  /* * @name {\bf STL Iterator Access Functions}.
225  *
226  * These member functions return valid STL random access iterators to the elements in the
227  * VectorSliceTmpl object.
228  *
229  * The forward iterators returned by begin() and end() iterator sequentialy from the first
230  * element (same element as returned by operator()(1)) to the last
231  * element (same element as returned by operator()(dim()). This goes for reverse
232  * (stride() < 0) VectorSliceTmpl objects as well. The reverse iterators returned by
233  * rbegin() and rend() iterate in the reverse sequence.
234  *
235  * Warning! Beware of using iterations in a reverse vector slice (stride() < 0).
236  * In a reverse vector slice end() returns a slice iterator which is the current
237  * element is one before the first allocated element. Strictly speaking this is
238  * not allowed so use iterators with reversed VectorSliceTmpl objects at own risk.
239  */
240 
241  // @{
242 
244  iterator begin();
246  iterator end();
248  const_iterator begin() const;
250  const_iterator end() const;
252  reverse_iterator rbegin();
254  reverse_iterator rend();
256  const_reverse_iterator rbegin() const;
258  const_reverse_iterator rend() const;
259 
260  // @}
261 
262  /* * @name {\bf Individual Element Access Subscripting (lvalue)}.
263  *
264  * These operator functions allow access (lvalue) to the individual elements
265  * of the VectorSliceTmpl object.
266  *
267  * The subscript i must be, 1 <= i <= this->dim(), for the 1-based element access
268  * operators and, 0 <= i <= this->dim() - 1, for the 0-based element access operators.
269  * If they are not then an #std::out_of_range# exception will be thrown.
270  */
271 
272  // @{
273 
275  reference operator()(size_type i);
277  const_reference operator()(size_type i) const;
279  reference operator[](size_type i);
281  const_reference operator[](size_type i) const;
282 
283  // @}
284 
285  /* * @name {\bf Subvector Access Operators}.
286  *
287  * These operator functions are used to create views of continous regions of the VectorSliceTmpl.
288  * Each of them returns a VectorSliceTmpl object for the region. Constant (const) VectorSliceTmpl objects
289  * are returned for a const VectorSliceTmpl. This means that the elements can not be changed
290  * as should be the case.
291  *
292  * Beware! VC++ is returning non-const VectorSliceTmpl objects for the
293  * #VectorSliceTmpl operator()(...) const;# member functions and therefore a const \Ref{DVector} or
294  * \Ref{VectorSliceTmpl} can be modifed my subsetting it. Hopefully this problem will
295  * be fixed in future versions of the compiler or I when will get another compiler.
296  */
297 
298  // @{
299 
301  /* * Returns a VectorSliceTmpl object representing the entire vector slice.
302  *
303  * Included for uniformity with vector.
304  */
307  return this;
308  }
311  return this;
312  }
313  VectorSliceTmpl<value_type>& operator()();
315  const VectorSliceTmpl<value_type>& operator()() const;
317  VectorSliceTmpl<value_type> operator()(const Range1D& rng);
319  /* * Returns a continous subregion of the VectorSliceTmpl object.
320  *
321  * The returned VectorSliceTmpl object represents the range of the rng argument.
322  *
323  * Preconditions: <ul>
324  * <li> #rng.ubound() - 1 <= this->dim()# (throw #out_of_range#)
325  * </ul>
326  *
327  * @param rng Indece range [lbound,ubound] of the region being returned.
328  */
329  const VectorSliceTmpl<value_type> operator()(const Range1D& rng) const;
331  /* * Returns a VectorSliceTmpl object for the continous subregion [ubound, lbound].
332  *
333  * Preconditions: <ul>
334  * <li> #lbound > 1# (throw out_of_range)
335  * <li> #lbound < ubound# (throw out_of_range)
336  * <li> #ubound <= this->dim()# (throw out_of_range)
337  * </ul>
338  *
339  * @param rng Range [lbound,ubound] of the region being returned.
340  */
341  VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound);
343  const VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound) const;
345  /* * Return a const VectorSliceTmpl object the reverse of this VectorSliceTmpl.
346  *
347  * In the reverse VectorSliceTmpl,
348  * the first element becomes the last element and visa-versa. For example, for
349  * #VectorSliceTmpl r = x.rev()#, #&x(1) == &z(z.dim())# and #&x(x.dim()) == &z(1)# are both true.
350  * The iterators returned by \Ref{begin()} iterate from the first conceptual element to the last.
351  */
354  const VectorSliceTmpl<value_type> rev() const;
355 
356  // @}
357 
358  /* * @name {\bf Assignment operators}. */
359 
360  // @{
361 
363  /* * vs = alpha (Sets all the elements to the constant alpha).
364  *
365  * Preconditions: <ul>
366  * <li> #this->dim() > 0# (throw #std::length_error#)
367  * </ul>
368  *
369  * Postconditions: <ul>
370  * <li> #this->operator()(i) == alpha#, i = 1, 2, ... , #this->dim()#
371  * </ul>
372  */
373  VectorSliceTmpl<value_type>& operator=(value_type alpha);
375  /* * vs = rhs (Copies the elements of rhs into the elements of this).
376  *
377  * Preconditions: <ul>
378  * <li> #this->dim() == rhs.dim()# (throw #out_of_range#)
379  * <li> #rhs.dim() > 0# (throw #out_of_range#)
380  * </ul>
381  *
382  * Postconditions: <ul>
383  * <li> #this->operator()(i) == rhs(i)#, i = 1, 2, ..., #this->dim()#
384  * </ul>
385  */
387 
388  // @}
389 
390  /* * @name {\bf Misc. Member Functions}. */
391 
392  // @{
393 
395  size_type dim() const;
397  /* * Returns the degree of memory overlap of the two VectorSliceTmpl objects this and vs.
398  *
399  * @return
400  * \begin{description}
401  * <li>[NO_OVERLAP] There is no memory overlap between this and vs
402  * <li>[SOME_OVERLAP] There is some memory locations that this and vs share
403  * <li>[SAME_MEM] The VectorSliceTmpl objects this and vs share the exact same memory locations.
404  * \end{description}
405  */
406  EOverLap overlap(const VectorSliceTmpl<value_type>& vs) const;
407 
408  // @}
409 
410  /* * @name {\bf Raw data access}.
411  *
412  * Provides access to underlying raw data.
413  */
414 
415  // @{
416 
418  value_type* raw_ptr();
420  const value_type* raw_ptr() const;
422  value_type* start_ptr();
424  const value_type* start_ptr() const;
426  difference_type stride() const;
427 
428  // @}
429 
430 private:
431  value_type *ptr_; // Pointer to first element in array.
432  size_type size_; // # elements represented in v_
433  difference_type stride_;// # positions to skip between elements. Must be positive
434  // Above the sequence represented is:
435  // ptr_[ i * stride_ ], for i = 0, ..., size_ - 1
436 
437 }; // end class VectorSliceTmpl<T>
438 
439 // /////////////////////////////////////////////////////////////////////////////////////////
440 // DVector
441 //
442 
444 /* * 1-D DVector Abstraction Storage Class.
445  *
446  * Holds the storage space for a 1-D vector of element type value_type. The storage space class
447  * used in a standard vector<> private member. DVector provides much of the
448  * same functionaliy of a VectorSliceTmpl object accept that DVector object can be resized at any time by
449  * either explicitly calling #resize(...)# or to match an assignment to a rhs linear algebra expression.
450  */
451 template<class T>
452 class VectorTmpl {
453 public:
454  /* * @name {\bf Nested Member Types (STL)}.
455  *
456  * These nested types give the types used in the interface to the class.
457  *
458  * \begin{description}
459  * <li>[#value_type#] - type being stored in the underlying vector<>
460  * <li>[#size_type#] - type for the number of elements in the vector<>
461  * <li>[#difference_type#] - type for the distance between elements
462  * <li>[#iterator#] - type for the forward non-constant iterator
463  * <li>[#const_iterator#] - type for the forward constant iterator (can't change elements)
464  * <li>[#reverse_iterator#] - type for the reverse non-constant iterator
465  * <li>[#const_reverse_iterator#] - type for the reverse constant iterator (can't change elements)
466  * <li>[#reference#] - type returned from subscripting, iterator deferencing etc.
467  * <li>[#const_reference#] - "" "" for const vector slice objects
468  * \end{description}
469  */
470 
471  // @{
472  // @}
473 
474  typedef T value_type;
475  typedef DenseLinAlgPack::size_type size_type;
476  typedef ptrdiff_t difference_type;
477  typedef value_type* iterator;
478  typedef const value_type* const_iterator;
479 #if 0 /* defined(_INTEL_CXX) || defined(_WINDOWS) */
480  typedef std::reverse_iterator<iterator, value_type
481  , value_type&, value_type*, difference_type> reverse_iterator;
482  typedef std::reverse_iterator<const_iterator
483  , value_type, const value_type&, const value_type*
484  , difference_type> const_reverse_iterator;
485 #else
486  typedef std::reverse_iterator<iterator> reverse_iterator;
487  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
488 #endif
489  typedef value_type& reference;
490  typedef const value_type& const_reference;
491  typedef std::vector<value_type> valarray;
492 
493  /* * @name {\bf Constructors}.
494  *
495  * These constructors allocate and may initialize the elements of a 1-D vector.
496  * The default C++ copy constructor is used and is therefore not show here.
497  */
498 
499  // @{
500 
502  VectorTmpl();
504  VectorTmpl(size_type n);
506  VectorTmpl(value_type val, size_type n);
508  /* * Constructs a vector with n elements and intializes elements to those of an array.
509  *
510  * Postconditions: <ul>
511  * <li> #this->operator[](i) == p[i]#, i = 0, 1, ... n
512  * </ul>
513  */
514  VectorTmpl(const value_type* p, size_type n);
516  /* * Constructs a DVector object fron a VectorSliceTmpl object.
517  *
518  * Postconditions: <ul>
519  * <li> #this->dim() == vs.dim()#
520  * <li> #this->operator[](i) == vs[i]#, i = 0, 1, ... n
521  * </ul>
522  */
524 
525  // @}
526 
527  /* * @name {\bf Memory Management / Misc}. */
528 
529  // @{
530 
532  /* * Resize the vector to hold n elements.
533  *
534  * Any new elements added are initialized to val.
535  *
536  * Postconditions: <ul>
537  * <li> #this->dim() == n#
538  * </ul>
539  */
540  void resize(size_type n, value_type val = value_type());
542  /* * Free memory and resize DVector to this->dim() == 0.
543  *
544  * Postconditions: <ul>
545  * <li> #this->dim() == 0#
546  * </ul>
547  */
548  void free();
550  size_type dim() const;
552  /* * Returns the degree of memory overlap of this and the VectorSliceTmpl object vs.
553  *
554  * @return
555  * \begin{description}
556  * <li>[NO_OVERLAP] There is no memory overlap between this and vs
557  * <li>[SOME_OVERLAP] There is some memory locations that this and vs share
558  * <li>[SAME_MEM] The VectorSliceTmpl objects this and vs share the exact same memory locations.
559  * \end{description}
560  */
561  EOverLap overlap(const VectorSliceTmpl<value_type>& vs) const;
563  operator VectorSliceTmpl<value_type>();
565  operator const VectorSliceTmpl<value_type>() const;
566 
567  // @}
568 
569 
570  /* * @name {\bf STL Iterator Access functions}.
571  *
572  * The iterators returned are valid STL random access iterators.
573  * The forward iterators returned iterate from the first element to the last element.
574  * The reverse iterators returned iterate from the last element to the first element.
575  */
576 
577  // @{
578 
580  iterator begin();
582  iterator end();
584  const_iterator begin() const;
586  const_iterator end() const;
588  reverse_iterator rbegin();
590  reverse_iterator rend();
592  const_reverse_iterator rbegin() const;
594  const_reverse_iterator rend() const;
595 
596  // @}
597 
598  /* * @name {\bf Individual Element Access Subscripting (lvalue)}.
599  *
600  * These operator functions allow access (lvalue) to the individual elements
601  * of the DVector object.
602  *
603  * The subscript i must be, 1 <= i <= this->dim(), for the 1-based element access
604  * operators and, 0 <= i <= this->dim() - 1, for the 0-based element access operators.
605  * If they are not then an #std::out_of_range# exception will be thrown.
606  */
607 
608  // @{
609 
611  reference operator()(size_type i);
613  const_reference operator()(size_type i) const;
615  reference operator[](size_type i);
617  const_reference operator[](size_type i) const;
618 
619  // @}
620 
621  /* * @name {\bf Subvector Access Operators}.
622  *
623  * These operator functions are used to create views of continous regions of the DVector.
624  * Each of them returns a VectorSliceTmpl object for the region. Constant (const) VectorSliceTmpl objects
625  * are returned for a const DVector. This means that the elements can not be changed
626  * as should be the case.
627  *
628  * Beware! VC++ is returning non-const VectorSliceTmpl objects for the
629  * #VectorSliceTmpl operator()(...) const;# member functions and therefore a const \Ref{DVector} or
630  * \Ref{VectorSliceTmpl} can be modifed my subsetting it. Hopefully this problem will
631  * be fixed in future versions of the compiler or I when will get another compiler.
632  */
633 
634  // @{
635 
637  /* * Returns a VectorSliceTmpl object representing the entire DVector.
638  *
639  * Call this member function to force a type conversion to VectorSliceTmpl. Using the
640  * VectorSliceTmpl of a DVector for algebraic expressions used with the TCOL allows a for simplier
641  * implementaion of those operations by cutting down on the number combinations. This is
642  * especialy true for longer optimized expression.
643  */
644  VectorSliceTmpl<value_type> operator()();
646  const VectorSliceTmpl<value_type> operator()() const;
648  /* * Returns a continous subregion of the DVector object.
649  *
650  * The returned VectorSliceTmpl object represents the range of the rng argument.
651  *
652  * Preconditions: <ul>
653  * <li> #rng.ubound() - 1 <= this->dim()# (throw #out_of_range#)
654  * </ul>
655  *
656  * @param rng Indece range [lbound,ubound] of the region being returned.
657  */
658  VectorSliceTmpl<value_type> operator()(const Range1D& rng);
660  const VectorSliceTmpl<value_type> operator()(const Range1D& rng) const;
662  /* * Returns a VectorSliceTmpl object for the continous subregion [ubound, lbound].
663  *
664  * Preconditions: <ul>
665  * <li> #lbound > 1# (throw #out_of_range#)
666  * <li> #lbound < ubound# (throw #out_of_range#)
667  * <li> #ubound <= this->dim()# (throw #out_of_range#)
668  * </ul>
669  *
670  * @param rng Range [lbound,ubound] of the region being taken.
671  */
672  VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound);
674  const VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound) const;
676  /* * Return a VectorSliceTmpl object the reverse of this DVector.
677  *
678  * In the reverse VectorSliceTmpl,
679  * the first element becomes the last element and visa-versa. For example, for
680  * #VectorSliceTmpl r = x.rev()#, #&x(1) == &z(z.dim())# and #&x(x.dim()) == &z(1)# are both true.
681  * The iterators returned by \Ref{begin()} iterate from the first conceptual element to the last.
682  */
685  const VectorSliceTmpl<value_type> rev() const;
686 
687  // @}
688 
689  /* * @name {\bf Assignment Operators}. */
690 
691  // @{
692 
694  /* * vs = alpha (Sets all the elements to the constant alpha).
695  *
696  * Preconditions: <ul>
697  * <li> #this->dim() > 0# (throw #std::length_error#)
698  * </ul>
699  *
700  * Postconditions: <ul>
701  * <li> #this->operator()(i) == alpha#, i = 1, 2, ... , #this->dim()#
702  * </ul>
703  */
704  VectorTmpl<value_type>& operator=(value_type alpha);
706  /* * vs = rhs (Copies the elements of rhs into the elements of this).
707  *
708  * Preconditions: <ul>
709  * <li> #this->dim() == rhs.dim()# (throw #out_of_range#)
710  * <li> #rhs.dim() > 0# (throw #out_of_range#)
711  * </ul>
712  *
713  * Postconditions: <ul>
714  * <li> #this->operator()(i) == rhs(i)#, i = 1, 2, ..., #this->dim()#
715  * </ul>
716  */
719  /* * Needed to override the default assignment operator.
720  */
722 
723  // @}
724 
725  /* * @name {\bf Implementation Access}.
726  *
727  * Provides access to underlying raw data.
728  */
729 
730  // @{
731 
733  value_type* raw_ptr();
735  const value_type* raw_ptr() const;
737  value_type* start_ptr();
739  const value_type* start_ptr() const;
741  difference_type stride() const;
742 
743  // @}
744 
745 private:
746  valarray v_;
747 
748 }; // end class VectorTmpl<T>
749 
750 // end DVector Classes scope
751 // @}
752 
753 // ///////////////////////////////////////////////////////////////////////////////
754 // Non-member function declarations //
755 // ///////////////////////////////////////////////////////////////////////////////
756 
757 
758 /* * @name {\bf Non-Member Functions}. */
759 
760 // @{
761 // begin non-member functions scope
762 
763 //
764 size_type vector_validate_sized(size_type size);
765 //
766 void vector_validate_range(size_type ubound, size_type max_ubound);
767 //
768 void vector_validate_subscript(size_type size, size_type i);
770 /* * Utility for checking the sizes of two VectorSliceTmpl objects and throwing an exception
771  * if the sizes are not the same.
772  */
773 void assert_vs_sizes(size_type size1, size_type size2);
774 
776 /* * Create a general vector slice.
777  */
778 template<class T>
779 inline
780 VectorSliceTmpl<T> gen_vs( VectorSliceTmpl<T>& vs, size_type start, size_type size, ptrdiff_t stride )
781 {
782  return VectorSliceTmpl<T>( vs.start_ptr() + vs.stride() * (start-1), size, vs.stride() * stride );
783 }
784 
786 template<class T>
787 inline
788 const VectorSliceTmpl<T> gen_vs( const VectorSliceTmpl<T>& vs, size_type start, size_type size
789  , ptrdiff_t stride )
790 {
791  return VectorSliceTmpl<T>( const_cast<typename VectorSliceTmpl<T>::value_type*>(vs.start_ptr()) + vs.stride() * (start-1)
792  , size, vs.stride() * stride );
793 }
794 
795 // end non-member functions scope
796 // @}
797 
798 // end Vectors scope
799 // @}
800 
801 } // end namespace DenseLinAlgPack
802 
803 // ////////////////////////////////////////////////////////////////////////////////
804 // Inline definitions of member function definitions //
805 // ////////////////////////////////////////////////////////////////////////////////
806 
807 // ////////////////////////////////////////////////////////////////////////
808 // Non-member functions / utilities
809 
810 #ifndef LINALGPACK_CHECK_SLICE_SETUP
811 inline
812 DenseLinAlgPack::size_type DenseLinAlgPack::vector_validate_sized(size_type size)
813 {
814  return size;
815 }
816 #endif
817 
818 #ifndef LINALGPACK_CHECK_RANGE
819 inline
820 void DenseLinAlgPack::vector_validate_range(size_type ubound, size_type max_ubound)
821 {}
822 #endif
823 
824 #ifndef LINALGPACK_CHECK_RANGE
825 inline
826 void DenseLinAlgPack::vector_validate_subscript(size_type size, size_type i)
827 {}
828 #endif
829 
830 #ifndef LINALGPACK_CHECK_RHS_SIZES
831 inline
832 void DenseLinAlgPack::assert_vs_sizes(size_type size1, size_type size2)
833 {}
834 #endif
835 
836 namespace DenseLinAlgPack {
837 
838 // /////////////////////////////////////////////////////////////////////////////
839 // VectorSliceTmpl inline member function definitions
840 
841 // Constructors. Use default copy constructor
842 
843 template<class T>
844 inline
846  : ptr_(0)
847  , size_(0)
848  , stride_(0)
849 {}
850 
851 template<class T>
852 inline
853 VectorSliceTmpl<T>::VectorSliceTmpl( value_type* ptr, size_type size, difference_type stride)
854  : ptr_(ptr)
855  , size_(size)
856  , stride_(stride)
857 {}
858 
859 template<class T>
860 inline
861 VectorSliceTmpl<T>::VectorSliceTmpl( value_type* ptr, size_type size, const Range1D& rng )
862  : ptr_( ptr + rng.lbound() - 1 )
863  , size_( rng.full_range() ? vector_validate_sized(size) : rng.size() )
864  , stride_(1)
865 {
866  vector_validate_range( rng.full_range() ? size : rng.ubound(), size );
867 }
868 
869 template<class T>
870 inline
872  : ptr_( vs.start_ptr() + (rng.lbound() - 1) * vs.stride() )
873  , size_( rng.full_range() ? vector_validate_sized(vs.dim()) : rng.size() )
874  , stride_( vs.stride() )
875 { vector_validate_range( rng.full_range() ? vs.dim() : rng.ubound(), vs.dim() ); }
876 
877 template<class T>
878 inline
880 {
881  ptr_ = vs.ptr_;
882  size_ = vs.size_;
883  stride_ = vs.stride_;
884 }
885 
886 // Iterator functions
887 template<class T>
888 inline
890 { return iterator(start_ptr(), stride()); }
891 
892 template<class T>
893 inline
895 { return iterator(start_ptr() + dim() * stride(), stride()); }
896 
897 template<class T>
898 inline
900 { return const_iterator(start_ptr(), stride()); }
901 
902 template<class T>
903 inline
905 { return const_iterator(start_ptr() + dim() * stride(), stride()); }
906 
907 template<class T>
908 inline
909 typename VectorSliceTmpl<T>::reverse_iterator VectorSliceTmpl<T>::rbegin()
910 { return reverse_iterator(end()); }
911 
912 template<class T>
913 inline
914 typename VectorSliceTmpl<T>::reverse_iterator VectorSliceTmpl<T>::rend()
915 { return reverse_iterator(begin()); }
916 
917 template<class T>
918 inline
919 typename VectorSliceTmpl<T>::const_reverse_iterator VectorSliceTmpl<T>::rbegin() const
920 { return const_reverse_iterator(end()); }
921 
922 template<class T>
923 inline
924 typename VectorSliceTmpl<T>::const_reverse_iterator VectorSliceTmpl<T>::rend() const
925 { return const_reverse_iterator(begin()); }
926 
927 // Element access
928 template<class T>
929 inline
930 typename VectorSliceTmpl<T>::reference VectorSliceTmpl<T>::operator()(size_type i) // 1 based
931 {
932  vector_validate_subscript(dim(),i);
933  return ptr_[(i-1)*stride_];
934 }
935 
936 template<class T>
937 inline
938 typename VectorSliceTmpl<T>::const_reference VectorSliceTmpl<T>::operator()(size_type i) const
939 {
940  vector_validate_subscript(dim(),i);
941  return ptr_[(i-1)*stride_];
942 }
943 
944 template<class T>
945 inline
946 typename VectorSliceTmpl<T>::reference VectorSliceTmpl<T>::operator[](size_type i) // 0 based
947 {
948  vector_validate_subscript(dim(),i+1);
949  return ptr_[(i)*stride_];
950 }
951 
952 template<class T>
953 inline
954 typename VectorSliceTmpl<T>::const_reference VectorSliceTmpl<T>::operator[](size_type i) const
955 {
956  vector_validate_subscript(dim(),i+1);
957  return ptr_[(i)*stride_];
958 }
959 
960 // Subregion Access. Let the constructors of VectorSliceTmpl validate the ranges
961 template<class T>
962 inline
964 { return *this; }
965 
966 template<class T>
967 inline
969 { return *this; }
970 
971 template<class T>
972 inline
974 { return VectorSliceTmpl(*this, RangePack::full_range(rng,1,dim())); }
975 
976 template<class T>
977 inline
978 const VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(const Range1D& rng) const
979 { return VectorSliceTmpl(const_cast<VectorSliceTmpl<T>&>(*this), RangePack::full_range(rng,1,dim())); }
980 
981 template<class T>
982 inline
983 VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(size_type lbound, size_type ubound)
984 { return VectorSliceTmpl(*this, Range1D(lbound, ubound)); }
985 
986 template<class T>
987 inline
988 const VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(size_type lbound, size_type ubound) const
989 { return VectorSliceTmpl(const_cast<VectorSliceTmpl<T>&>(*this), Range1D(lbound, ubound)); }
990 
991 template<class T>
992 inline
994 { return VectorSliceTmpl( start_ptr() + stride() * (dim()-1), dim(), - stride() ); }
995 
996 template<class T>
997 inline
999 { return VectorSliceTmpl( const_cast<value_type*>(start_ptr()) + stride() * (dim()-1), dim(), - stride() ); }
1000 
1001 // Assignment Operators
1002 template<class T>
1003 inline
1005 {
1006  std::fill(begin(),end(),alpha);
1007  return *this;
1008 }
1009 
1010 template<class T>
1011 inline
1013 {
1014  assert_vs_sizes(this->dim(),rhs.dim());
1015  std::copy(rhs.begin(),rhs.end(),begin());
1016  return *this;
1017 }
1018 
1019 // Misc. member functions
1020 
1021 template<class T>
1022 inline
1023 typename VectorSliceTmpl<T>::size_type VectorSliceTmpl<T>::dim() const
1024 { return size_; }
1025 
1026 // Raw pointer access
1027 
1028 template<class T>
1029 inline
1030 typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::raw_ptr()
1031 { return stride() > 0 ? start_ptr() : start_ptr() + stride() * (dim() - 1); }
1032 
1033 template<class T>
1034 inline
1035 const typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::raw_ptr() const
1036 { return stride() > 0 ? start_ptr() : start_ptr() + stride() * (dim() - 1); }
1037 
1038 template<class T>
1039 inline
1040 typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::start_ptr()
1041 { return ptr_; }
1042 
1043 template<class T>
1044 inline
1045 const typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::start_ptr() const
1046 { return ptr_; }
1047 
1048 template<class T>
1049 inline
1050 typename VectorSliceTmpl<T>::difference_type VectorSliceTmpl<T>::stride() const
1051 { return stride_; }
1052 
1053 // /////////////////////////////////////////////////////////////////////////////
1054 // DVector inline member function definitions
1055 
1056 // Constructors
1057 template<class T>
1058 inline
1060 {} // used to shut satisfy compiler
1061 
1062 template<class T>
1063 inline
1065  : v_(n)
1066 {}
1067 
1068 template<class T>
1069 inline
1070 VectorTmpl<T>::VectorTmpl(value_type val, size_type n)
1071  : v_(n)
1072 {
1073  std::fill(begin(),end(),val);
1074 }
1075 
1076 template<class T>
1077 inline
1078 VectorTmpl<T>::VectorTmpl(const value_type* p, size_type n)
1079  : v_(n)
1080 {
1081  std::copy(p,p+n,begin());
1082 }
1083 
1084 template<class T>
1085 inline
1087  : v_(vs.dim())
1088 {
1089  std::copy(vs.begin(),vs.end(),begin());
1090 }
1091 
1092 // Memory management
1093 template<class T>
1094 inline
1095 void VectorTmpl<T>::resize(size_type n, value_type val)
1096 {
1097  v_.resize(n);
1098  std::fill(begin(),end(),val);
1099  }
1100 
1101 template<class T>
1102 inline
1104 {
1105  v_.resize(0);
1106 }
1107 
1108 // Size
1109 template<class T>
1110 inline
1111 typename VectorTmpl<T>::size_type VectorTmpl<T>::dim() const
1112 { return v_.size(); }
1113 
1114 // Iterator functions
1115 template<class T>
1116 inline
1117 typename VectorTmpl<T>::iterator VectorTmpl<T>::begin()
1118 { return start_ptr(); }
1119 
1120 template<class T>
1121 inline
1122 typename VectorTmpl<T>::iterator VectorTmpl<T>::end()
1123 { return start_ptr() + dim(); }
1124 
1125 template<class T>
1126 inline
1127 typename VectorTmpl<T>::const_iterator VectorTmpl<T>::begin() const
1128 { return start_ptr(); }
1129 
1130 template<class T>
1131 inline
1132 typename VectorTmpl<T>::const_iterator VectorTmpl<T>::end() const
1133 { return start_ptr() + dim(); }
1134 
1135 template<class T>
1136 inline
1137 typename VectorTmpl<T>::reverse_iterator VectorTmpl<T>::rbegin()
1138 { return reverse_iterator(end()); }
1139 
1140 template<class T>
1141 inline
1142 typename VectorTmpl<T>::reverse_iterator VectorTmpl<T>::rend()
1143 { return reverse_iterator(begin()); }
1144 
1145 template<class T>
1146 inline
1147 typename VectorTmpl<T>::const_reverse_iterator VectorTmpl<T>::rbegin() const
1148 { return const_reverse_iterator(end()); }
1149 
1150 template<class T>
1151 inline
1152 typename VectorTmpl<T>::const_reverse_iterator VectorTmpl<T>::rend() const
1153 { return const_reverse_iterator(begin()); }
1154 
1155 // Element access
1156 template<class T>
1157 inline
1158 typename VectorTmpl<T>::reference VectorTmpl<T>::operator()(size_type i)
1159 {
1160  vector_validate_subscript(dim(),i);
1161  return start_ptr()[i-1];
1162 }
1163 
1164 template<class T>
1165 inline
1166 typename VectorTmpl<T>::const_reference VectorTmpl<T>::operator()(size_type i) const
1167 {
1168  vector_validate_subscript(dim(),i);
1169  return start_ptr()[i-1];
1170 }
1171 
1172 template<class T>
1173 inline
1174 typename VectorTmpl<T>::reference VectorTmpl<T>::operator[](size_type i)
1175 {
1176  vector_validate_subscript(dim(),i+1);
1177  return start_ptr()[i];
1178 }
1179 
1180 template<class T>
1181 inline
1182 typename VectorTmpl<T>::const_reference VectorTmpl<T>::operator[](size_type i) const
1183 {
1184  vector_validate_subscript(dim(),i+1);
1185  return start_ptr()[i];
1186 }
1187 
1188 // Subregion Access. Leave validation to the VectorSliceTmpl constructors.
1189 template<class T>
1190 inline
1192 { return VectorSliceTmpl<T>(start_ptr(),dim()); }
1193 
1194 template<class T>
1195 inline
1197 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()),dim()); }
1198 
1199 template<class T>
1200 inline
1202 { return VectorSliceTmpl<T>(start_ptr(),dim(),rng); }
1203 
1204 template<class T>
1205 inline
1207 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()),dim(),rng); }
1208 
1209 template<class T>
1210 inline
1211 VectorSliceTmpl<T> VectorTmpl<T>::operator()(size_type lbound, size_type ubound)
1212 { return VectorSliceTmpl<T>(start_ptr(), dim(), Range1D(lbound, ubound)); }
1213 
1214 template<class T>
1215 inline
1216 const VectorSliceTmpl<T> VectorTmpl<T>::operator()(size_type lbound, size_type ubound) const
1217 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()), dim(), Range1D(lbound, ubound)); }
1218 
1219 template<class T>
1220 inline
1222 { return VectorSliceTmpl<T>( start_ptr() + dim() - 1, dim(), -1 ); }
1223 
1224 template<class T>
1225 inline
1227 { return VectorSliceTmpl<T>( const_cast<value_type*>(start_ptr()) + dim() - 1, dim(), -1 ); }
1228 
1229 // Conversion operators
1230 template<class T>
1231 inline
1233 { return VectorSliceTmpl<T>(start_ptr(), dim()); }
1234 
1235 template<class T>
1236 inline
1237 VectorTmpl<T>::operator const VectorSliceTmpl<T>() const
1238 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()), dim()); }
1239 
1240 // Assignment Operators
1241 template<class T>
1242 inline
1244 {
1245  if(!dim()) resize(1);
1246  std::fill(begin(),end(),alpha);
1247  return *this;
1248 }
1249 
1250 template<class T>
1251 inline
1253 {
1254  resize(rhs.dim());
1255  std::copy(rhs.begin(),rhs.end(),begin());
1256  return *this;
1257 }
1258 
1259 template<class T>
1260 inline
1262 {
1263  resize(rhs.dim());
1264  std::copy(rhs.begin(),rhs.end(),begin());
1265  return *this;
1266 }
1267 
1268 // Raw pointer access
1269 
1270 template<class T>
1271 inline
1272 typename VectorTmpl<T>::value_type* VectorTmpl<T>::raw_ptr()
1273 { return start_ptr(); }
1274 
1275 template<class T>
1276 inline
1277 const typename VectorTmpl<T>::value_type* VectorTmpl<T>::raw_ptr() const
1278 { return start_ptr(); }
1279 
1280 template<class T>
1281 inline
1282 typename VectorTmpl<T>::value_type* VectorTmpl<T>::start_ptr()
1283 { return dim() ? &(v_)[0] : 0; }
1284 
1285 template<class T>
1286 inline
1287 const typename VectorTmpl<T>::value_type* VectorTmpl<T>::start_ptr() const
1288 { return &const_cast<valarray&>((v_))[0]; }
1289 
1290 template<class T>
1291 inline
1292 typename VectorTmpl<T>::difference_type VectorTmpl<T>::stride() const
1293 { return 1; }
1294 
1295 // //////////////////////////////////////////////////
1296 // Non-inlined members
1297 
1298 // Assume that the vector slices are the rows, cols or diag of a 2-D matrix.
1299 template<class T>
1300 EOverLap VectorSliceTmpl<T>::overlap(const VectorSliceTmpl<value_type>& vs) const {
1301 
1302  const typename VectorSliceTmpl<T>::value_type
1303  *raw_ptr1 = ( stride() > 0 ? start_ptr() : start_ptr() + (dim()-1)*stride() ),
1304  *raw_ptr2 = ( vs.stride() > 0 ? vs.start_ptr() : vs.start_ptr() + (vs.dim()-1)*vs.stride() );
1305  typename VectorSliceTmpl<T>::size_type
1306  size1 = dim(),
1307  size2 = vs.dim();
1308  typename VectorSliceTmpl<T>::difference_type
1309  stride1 = std::abs(stride()),
1310  stride2 = std::abs(vs.stride());
1311 
1312  // Establish a frame of reference where raw_ptr1 < raw_ptr2
1313  if(raw_ptr1 > raw_ptr2) {
1314  std::swap(raw_ptr1,raw_ptr2);
1315  std::swap(stride1,stride2);
1316  std::swap(size1,size2);
1317  }
1318 
1319  if( raw_ptr1 + stride1 * (size1 - 1) < raw_ptr2 ) {
1320  return NO_OVERLAP; // can't be any overlap
1321  }
1322 
1324  start1 = 0,
1325  start2 = raw_ptr2 - raw_ptr1;
1326 
1327  if(start1 == start2 && stride1 == stride2 && size1 == size2)
1328  return SAME_MEM;
1329 // else if(start1 == start2)
1330 // return SOME_OVERLAP; // First elements are the same
1331 // else if(stride1 + (size1 - 1) * stride1 == stride2 + (size2 - 1) * stride2)
1332 // return SOME_OVERLAP; // Last elements are the same
1333  else if(stride1 == stride2) {
1334  if(!((start2 - start1) % stride1))
1335  return SOME_OVERLAP;
1336  else
1337  return NO_OVERLAP; // a different row, col or diag of a matrix
1338  }
1339  else {
1340  if(stride1 == 1 || stride2 == 1) {
1341  // One of them is a column vector.
1342  // Make vs1 the column vector.
1343  bool switch_them = (stride2 == 1);
1344  if(switch_them) {
1345  std::swap(start1,start2);
1346  std::swap(stride1,stride2);
1347  std::swap(size1,size2);
1348  }
1349 
1350  // Determine if the other vector could be row vector
1351  // or must be a diag vector. If using stride2 makes
1352  // the first and last elements of the column vector
1353  // on different rows, then the other vector must be a diagonal.
1354  // col_first = start1/stride2, col_last = (start1+size1-1)/stride2
1355  // if(col_last - col_first > 0) then vs2 must be a diagonal vector
1356  // with max_rows = stride2 - 1.
1357  size_t max_rows = (start1+size1-1)/stride2 - start1/stride2 > 0 ? stride2 - 1 : stride2;
1358 
1359  // find the index (0-based) of vs2 that intersects this column.
1360  size_t vs2_col_i = start1/max_rows - start2/max_rows;
1361 
1362  // See if the first element of the column is above the element in vs2
1363  // and the last element of the column is below the element. If it is
1364  // then we conclude that there is an itersection.
1365  size_t vs2_col_rng = start2 + vs2_col_i * stride2;
1366  if(start1 <= vs2_col_rng && vs2_col_rng <= start1+size1-1)
1367  return SOME_OVERLAP;
1368  else
1369  return NO_OVERLAP;
1370  }
1371  // They are not the same and nether is a column vector so one is a row vector
1372  // and the other is a diagonal vector.
1373  // Nether is a column vector so choose as vs1 the row vector (the smaller stride).
1374  bool switch_them = stride2 < stride1;
1375  if(switch_them) {
1376  std::swap(start1,start2);
1377  std::swap(stride1,stride2);
1378  std::swap(size1,size2);
1379  }
1380 
1381  size_t max_rows = stride1;
1382  // Determine the first and last columns (0-based) in the original
1383  // matrix where there vs1 and vs2 intersect.
1384  size_t sec_first_col = (start1 > start2) ? start1/max_rows : start2/max_rows,
1385  last1 = start1 + (size1 - 1) * stride1,
1386  last2 = start2 + (size2 - 1) * stride2,
1387  sec_last_col = (last1 < last2) ? last1/max_rows : last2/max_rows;
1388  // Determine the vector indexes (0-based) of vs1 and vs2 for the start and end
1389  // in this region
1390  size_t vs1_first_col = start1 / max_rows,
1391  vs2_first_col = start2 / max_rows;
1392 
1393  // Determine the indexes in the valarray of the two vectors for the two ends
1394  size_t vs1_first_col_i = sec_first_col - vs1_first_col,
1395  vs1_last_col_i = sec_last_col - vs1_first_col,
1396  vs2_first_col_i = sec_first_col - vs2_first_col,
1397  vs2_last_col_i = sec_last_col - vs2_first_col;
1398 
1399  // Compare the indexes in the valarray at the two ends. If they cross then
1400  // there must be an element of overlap. Uses equivalent of the intermediate
1401  // value therorm.
1402  // Must cast to an int that can hold a negative value
1403  ptrdiff_t diff1 = (start1 + vs1_first_col_i * stride1)
1404  - static_cast<ptrdiff_t>((start2 + vs2_first_col_i * stride2)),
1405  diff2 = (start1 + vs1_last_col_i * stride1)
1406  - static_cast<ptrdiff_t>((start2 + vs2_last_col_i * stride2));
1407  if(diff1 * diff2 > 0 )
1408  return NO_OVERLAP; // they do not cross
1409  else
1410  return SOME_OVERLAP; // they share an element
1411  }
1412 }
1413 
1414 template<class T>
1415 EOverLap VectorTmpl<T>::overlap(const VectorSliceTmpl<value_type>& vs) const {
1416 
1417  const typename VectorSliceTmpl<T>::value_type
1418  *raw_ptr1 = ( stride() > 0 ? start_ptr() : start_ptr() + (dim()-1)*stride() ),
1419  *raw_ptr2 = ( vs.stride() > 0 ? vs.start_ptr() : vs.start_ptr() + (vs.dim()-1)*vs.stride() );
1421  size1 = dim(),
1422  size2 = vs.dim();
1423 
1424  if( raw_ptr1 <= raw_ptr2 && raw_ptr2 + size2 <= raw_ptr1 + size1 ) {
1425  if( raw_ptr1 == raw_ptr2 && size1 == size2 && 1 == vs.stride() )
1426  return SAME_MEM;
1427  return SOME_OVERLAP;
1428  }
1429  return NO_OVERLAP;
1430 }
1431 
1432 } // end namespace DenseLinAlgPack
1433 
1434 #endif // end VECTOR_CLASS_TMPL_H
size_type dim() const
Returns the number of elements of the VectorSliceTmpl.
value_type * start_ptr()
Return a pointer to the conceptual first element in the underlying array.
Index ubound() const
value_type * start_ptr()
Return a pointer to the conceptual first element in the underlying array.
void resize(size_type n, value_type val=value_type())
VectorTmpl()
Constructs a vector with 0 elements (this->dim()==0).
VectorSliceTmpl< value_type > * operator&()
Allow the address to be taken of an rvalue of this object.
value_type * raw_ptr()
Return a pointer to the address of the first memory location of underlying array. ...
VectorSliceTmpl< value_type > & operator=(value_type alpha)
reference operator()(size_type i)
1-based element access (lvalue)
size_t size_type
value_type * raw_ptr()
Return a pointer to the address of the first memory location of underlying array. ...
bool full_range() const
size_type dim() const
Returns the number of elements of the DVector.
difference_type stride() const
Return the distance (+,-) (in units of elements) between adjacent elements in the underlying array...
void bind(VectorSliceTmpl< value_type > vs)
Bind to the view of another VectorSliceTmpl.
const VectorSliceTmpl< value_type > * operator&() const
difference_type stride() const
Return the distance (+,-) (in units of elements) between adjacent elements in the underlying array...
reference operator[](size_type i)
1-based element access (lvalue)
reference operator()(size_type i)
1-based element access (lvalue)
reference operator[](size_type i)
1-based element access (lvalue)
VectorTmpl< value_type > & operator=(value_type alpha)