MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_DVectorClassTmpl.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef VECTOR_CLASS_TMPL_H
00043 #define VECTOR_CLASS_TMPL_H
00044 
00045 #include <vector>
00046 
00047 #include "DenseLinAlgPack_Types.hpp"
00048 #include "StrideIterPack_StrideIter.hpp"
00049 
00050 namespace DenseLinAlgPack{
00051 
00052 // ////////////////////////////////////////////////////////////////////////////////
00053 /* * @name {\bf Dense 1-D DVector Abstractions}.
00054   *
00055   * These are classes that abstract 1-D vectors. The class \Ref{DVector} is a storage class
00056   * for vectors while the class \Ref{DVectorSlice} is used to represent regions of vectors
00057   * , for rows, columns or diagonals of matrices (see \Ref{DMatrix}
00058   * and \Ref{DMatrixSlice}).
00059   */
00060 
00061 // @{
00062 
00063 /* * @name {\bf DVector Classes}. */
00064 // @{
00065   
00067 /* * Slice of a 1-D sequential C++ array treated as a vector.
00068   *
00069   * Objects of this class represent regions of vectors (continuous), rows of matrices
00070   * , columns of matrices or diagonals of matrices.
00071   * The underlying representation is of a continuous C++ array with non unit stride.
00072   * It uses the same convention that the BLAS use where a vector is represented as
00073   * the first element of
00074   * in an array, the stride between elements in that array and the number of elements.
00075   * Changes to elements through a DVectorSlice object result in changes to the elements
00076   * in the underlying value_type* data.
00077   *
00078   * DVectorSlice provides many STL compliant features such as typedef type members
00079   *, iterator returning functions
00080   * and the dim() function.  It also provides access to individual elements (lvalue)
00081   * through 0-based
00082   * and 1-based subscripting with operator[](i) and operator()(i) respectively.
00083   *  In addition and subregions can be created
00084   * by subscripting (operator()()) with an \Ref{Range1D} object or using lower (>0)
00085   * and upper bounds of the region.
00086   */
00087 template<class T>
00088 class VectorSliceTmpl {
00089 public:
00090   
00091   /* * @name {\bf Nested Member Types (STL)}.
00092     *
00093     * These nested types give the types used in the interface to the class.
00094     *
00095     * \begin{description}
00096     * <li>[#value_type#]        - type being stored in the underlying C++ array     
00097     * <li>[#size_type#]       - type used as an index and for the number of elements
00098     *                   in the vector
00099     * <li>[#difference_type#]   - type for the distance between elements and the stride
00100     * <li>[#iterator#]        - type for the forward non-constant iterator
00101     * <li>[#const_iterator#]      - type for the forward constant iterator (can't change elements)
00102     * <li>[#reverse_iterator#]    - type for the reverse non-constant iterator
00103     * <li>[#const_reverse_iterator#]  - type for the reverse constant iterator (can't change elements)
00104     * <li>[#reference#]       - type returned from subscripting, iterator deferencing etc.
00105     * <li>[#const_reference#]   - "" "" for const vector slice objects
00106     * \end{description}
00107     */
00108 
00109   // @{
00110   // @}
00111 
00112   typedef T                       value_type;
00113   typedef DenseLinAlgPack::size_type              size_type;
00114   typedef ptrdiff_t                   difference_type;
00115   typedef StrideIterPack::stride_iter<value_type*
00116     , value_type, value_type&, value_type*
00117     , difference_type>                  iterator;
00118   typedef StrideIterPack::stride_iter<const value_type*
00119     , value_type, const value_type&, const value_type*
00120     , difference_type>                  const_iterator;
00121 #if defined(_INTEL_CXX) || defined (_INTEL_CXX)
00122   typedef std::reverse_iterator<iterator, value_type
00123     , value_type&, value_type*, difference_type>    reverse_iterator;
00124   typedef std::reverse_iterator<const_iterator
00125     , value_type, const value_type&, const value_type*
00126     , difference_type>                  const_reverse_iterator;
00127 #else
00128   typedef std::reverse_iterator<iterator>         reverse_iterator;
00129   typedef std::reverse_iterator<const_iterator>     const_reverse_iterator;
00130 #endif
00131   typedef value_type&                   reference;
00132   typedef const value_type&               const_reference;
00133 
00134   /* * @name {\bf Constructors}.
00135     *
00136     * The user usually does not need not call any of these constructors
00137     * explicitly to create a vector slice.
00138     * These
00139     * constructors are used by the classes in the library to construct VectorSliceTmpl objects.
00140     * Instead, users create VectorSliceTmpl objects by indexing (\Ref{Range1D}) a \Ref{DVector}
00141     * , or \Ref{VectorSliceTmpl}
00142     * object or calling row(...), col(...) or diag(...) on a \Ref{DMatrix} or
00143     * \Ref{DMatrixSlice} object.
00144     * The default C++ copy constructor is used, and is therefore not show here.
00145     *
00146     * Constructors are also included for creating views of raw C++ arrays.
00147     * These constructors take non-#const# pointers.  However you can savely
00148     * create a #const# view of a #const# C++ array by using a constant cast.
00149     * For example:
00150     *
00151     \verbatim
00152     const VectorSliceTmpl<T>::size_type n = 5;
00153     const VectorSliceTmpl<T>::value_type ptr[n] = { 1.0, 2.0, 3.0, 4.0, 5,0 };
00154     const VectorSliceTmpl vec(const_cast<VectorSliceTmpl<T>::value_type*>(ptr),n);
00155     \endverbatim
00156     */
00157 
00158   // @{
00159 
00161   /* * Creates an empty view.
00162     *
00163     * You must use bind(...) to bind to a view to initialize after construction.
00164     */
00165   VectorSliceTmpl();
00167   /* * Creates a VectorSice object that represents a non-continous region of a raw C++ array.
00168     *
00169     * Of course the sequence of elements #ptr[stride * i]# for #i# = 0, 1, ..., #size#-1
00170     * must yield valid properly allocated regions of memory.  It is up the the user to insure
00171     * that they do.
00172     *
00173     * @param  ptr   Pointer to the first element in the raw C++ array
00174     * @param  size  number of elements in the vector slice
00175     * @param  stride  the distance (may be negative) between each successive element (default = 1)
00176     */
00177   VectorSliceTmpl( value_type* ptr, size_type size, difference_type stride = 1 );
00179   /* * Creates a VectorSliceTmpl object that represents a continous region of a raw C++ array.
00180     *
00181     * The VectorSliceTmpl Object represents the following elements of raw array:
00182     * 
00183     *      #ptr[rng.lbound()-1+i]#, for #i# = 0, 1, ..., #rng.ubound()-1#
00184     * 
00185     * Preconditions: <ul>
00186     *   <li> #rng.ubound() + 1 <= n# (throw std::out_of_range)
00187     *   </ul>
00188     *
00189     * @param  ptr   Pointer to the first element in the raw C++ array
00190     * @param  size  number of elements in the vector slice
00191     * @param  rng   index range (1-based) of the region being represented.  
00192     *         Here rng.full_range() can be true.
00193     */
00194   VectorSliceTmpl( value_type* ptr, size_type size, const Range1D& rng );
00196   /* * Create a VectorSliceTmpl that represents a continous region of the existing VectorSliceTmpl object, vs.
00197     *
00198     * The index, rng, is relative to the VectorSliceTmpl object, vs.
00199     * For example rng = [1,3] would create a VectorSliceTmpl object
00200     * representing the elements 2, 4 and 6.  The following
00201     * shows the elements represented by each of the participating objects.
00202     \verbatim
00203 
00204     vs =   [2, 4, 6, 8, 10]
00205     this = [2, 4, 6]
00206 
00207     \endverbatim
00208 
00209     * Preconditions: <ul>
00210     *   <li> rng.full_range() == false (throw #std::out_of_range#)
00211     *   <li> rng.dim() <= vs.dim() (throw #std::out_of_range#) 
00212     *   </ul>
00213     * 
00214     * @param  vs    VectorSliceTmpl object that this VectorSliceTmpl object is being created from
00215     * @param  rng   Range1D range of the vector slice being created.
00216     */
00217   VectorSliceTmpl( VectorSliceTmpl<value_type>& vs, const Range1D& rng );
00218 
00219   // @}
00220 
00222   void bind(VectorSliceTmpl<value_type> vs);
00223 
00224   /* * @name {\bf STL Iterator Access Functions}.
00225     *
00226     * These member functions return valid STL random access iterators to the elements in the
00227     * VectorSliceTmpl object.
00228     *
00229     * The forward iterators returned by begin() and end() iterator sequentialy from the first
00230     * element (same element as returned by operator()(1)) to the last
00231     * element (same element as returned by operator()(dim()).  This goes for reverse 
00232     * (stride() < 0) VectorSliceTmpl objects as well.  The reverse iterators returned by
00233     * rbegin() and rend() iterate in the reverse sequence.
00234     *
00235     * Warning! Beware of using iterations in a reverse vector slice (stride() < 0). 
00236     * In a reverse vector slice end() returns a slice iterator which is the current
00237     * element is one before the first allocated element.  Strictly speaking this is
00238     * not allowed so use iterators with reversed VectorSliceTmpl objects at own risk.
00239     */
00240 
00241   // @{
00242 
00244   iterator begin();
00246   iterator end();
00248   const_iterator begin() const;
00250   const_iterator end() const;
00252   reverse_iterator rbegin();
00254   reverse_iterator rend();
00256   const_reverse_iterator rbegin() const;
00258   const_reverse_iterator rend() const;
00259 
00260   // @}
00261 
00262   /* * @name {\bf Individual Element Access Subscripting (lvalue)}.
00263     *
00264     * These operator functions allow access (lvalue) to the individual elements
00265     * of the VectorSliceTmpl object.
00266     * 
00267     * The subscript i must be, 1 <= i <= this->dim(), for the 1-based element access
00268     * operators and, 0 <= i <= this->dim() - 1, for the 0-based element access operators.
00269     * If they are not then an #std::out_of_range# exception will be thrown.
00270     */
00271 
00272   // @{
00273 
00275   reference operator()(size_type i);
00277   const_reference operator()(size_type i) const;
00279   reference operator[](size_type i);
00281   const_reference operator[](size_type i) const;
00282 
00283   // @}
00284 
00285   /* * @name {\bf Subvector Access Operators}.
00286     *
00287     * These operator functions are used to create views of continous regions of the VectorSliceTmpl.
00288     * Each of them returns a VectorSliceTmpl object for the region.  Constant (const) VectorSliceTmpl objects
00289     * are returned for a const VectorSliceTmpl.  This means that the elements can not be changed
00290     * as should be the case.
00291     * 
00292     * Beware!  VC++ is returning non-const VectorSliceTmpl objects for the 
00293     * #VectorSliceTmpl operator()(...) const;# member functions and therefore a const \Ref{DVector} or
00294     * \Ref{VectorSliceTmpl} can be modifed my subsetting it.  Hopefully this problem will
00295     * be fixed in future versions of the compiler or I when will get another compiler.
00296     */
00297 
00298   // @{
00299 
00301   /* * Returns a VectorSliceTmpl object representing the entire vector slice.
00302     *
00303     * Included for uniformity with vector.
00304     */
00306   VectorSliceTmpl<value_type>* operator&() {
00307     return this;
00308   }
00310   const VectorSliceTmpl<value_type>* operator&() const {
00311     return this;
00312   }
00313   VectorSliceTmpl<value_type>& operator()();
00315   const VectorSliceTmpl<value_type>& operator()() const;
00317   VectorSliceTmpl<value_type> operator()(const Range1D& rng);
00319   /* * Returns a continous subregion of the VectorSliceTmpl object.
00320     *
00321     * The returned VectorSliceTmpl object represents the range of the rng argument.
00322     *
00323     * Preconditions: <ul>
00324     *   <li> #rng.ubound() - 1 <= this->dim()# (throw #out_of_range#)
00325     *   </ul>
00326     *
00327     * @param  rng   Indece range [lbound,ubound] of the region being returned.
00328     */
00329   const VectorSliceTmpl<value_type> operator()(const Range1D& rng) const;
00331   /* * Returns a VectorSliceTmpl object for the continous subregion [ubound, lbound].
00332     * 
00333     * Preconditions: <ul>
00334     *   <li> #lbound > 1# (throw out_of_range)
00335     *   <li> #lbound < ubound# (throw out_of_range)
00336     *   <li> #ubound <= this->dim()# (throw out_of_range)
00337     *   </ul>
00338     *
00339     * @param  rng   Range [lbound,ubound] of the region being returned.
00340     */
00341   VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound);
00343   const VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound) const;
00345   /* * Return a const VectorSliceTmpl object the reverse of this VectorSliceTmpl.
00346     *
00347     * In the reverse VectorSliceTmpl,
00348     * the first element becomes the last element and visa-versa.  For example, for 
00349     * #VectorSliceTmpl r = x.rev()#, #&x(1) == &z(z.dim())# and #&x(x.dim()) == &z(1)# are both true.
00350     * The iterators returned by \Ref{begin()} iterate from the first conceptual element to the last.
00351     */
00352   VectorSliceTmpl<value_type> rev();
00354   const VectorSliceTmpl<value_type> rev() const;
00355 
00356   // @}
00357   
00358   /* * @name {\bf Assignment operators}. */
00359 
00360   // @{
00361 
00363   /* * vs = alpha (Sets all the elements to the constant alpha).
00364     *
00365     * Preconditions: <ul>
00366     *   <li> #this->dim() > 0# (throw #std::length_error#)
00367     *   </ul>
00368     *
00369     * Postconditions: <ul>
00370     *   <li> #this->operator()(i) == alpha#, i = 1, 2, ... , #this->dim()#
00371     *   </ul>
00372     */
00373   VectorSliceTmpl<value_type>& operator=(value_type alpha);
00375   /* * vs = rhs (Copies the elements of rhs into the elements of this).
00376     *
00377     * Preconditions: <ul>
00378     *   <li> #this->dim() == rhs.dim()# (throw #out_of_range#)
00379     *   <li> #rhs.dim() > 0# (throw #out_of_range#)
00380     *   </ul>
00381     *
00382     * Postconditions: <ul>
00383     *   <li> #this->operator()(i) == rhs(i)#, i = 1, 2, ..., #this->dim()#
00384     *   </ul>
00385     */
00386   VectorSliceTmpl<value_type>& operator=(const VectorSliceTmpl<value_type>& rhs);
00387 
00388   // @}
00389 
00390   /* * @name {\bf Misc. Member Functions}. */
00391 
00392   // @{
00393 
00395   size_type dim() const;
00397   /* * Returns the degree of memory overlap of the two VectorSliceTmpl objects this and vs.
00398     *
00399     * @return 
00400     *   \begin{description}
00401     *   <li>[NO_OVERLAP]  There is no memory overlap between this and vs
00402     *   <li>[SOME_OVERLAP]  There is some memory locations that this and vs share
00403     *   <li>[SAME_MEM]    The VectorSliceTmpl objects this and vs share the exact same memory locations.
00404     *   \end{description}
00405     */
00406   EOverLap overlap(const VectorSliceTmpl<value_type>& vs) const;
00407 
00408   // @}
00409 
00410   /* * @name {\bf Raw data access}.
00411     *
00412     * Provides access to underlying raw data.
00413     */
00414 
00415   // @{
00416 
00418   value_type*     raw_ptr();
00420   const value_type* raw_ptr() const;
00422   value_type*     start_ptr();
00424   const value_type* start_ptr() const;
00426   difference_type   stride() const;
00427 
00428   // @}
00429 
00430 private:
00431   value_type          *ptr_;  // Pointer to first element in array.
00432   size_type         size_;  // # elements represented in v_
00433   difference_type       stride_;// # positions to skip between elements. Must be positive
00434   // Above the sequence represented is:
00435   //  ptr_[ i * stride_ ], for i = 0, ..., size_ - 1
00436 
00437 }; // end class VectorSliceTmpl<T>
00438 
00439 // /////////////////////////////////////////////////////////////////////////////////////////
00440 // DVector
00441 //
00442 
00444 /* * 1-D DVector Abstraction Storage Class.
00445   *
00446   * Holds the storage space for a 1-D vector of element type value_type.  The storage space class
00447   * used in a standard vector<> private member.  DVector provides much of the
00448   * same functionaliy of a VectorSliceTmpl object accept that DVector object can be resized at any time by
00449   * either explicitly calling #resize(...)# or to match an assignment to a rhs linear algebra expression.
00450   */
00451 template<class T>
00452 class VectorTmpl {
00453 public:
00454   /* * @name {\bf Nested Member Types (STL)}.
00455     *
00456     * These nested types give the types used in the interface to the class.
00457     *
00458     * \begin{description}
00459     * <li>[#value_type#]        - type being stored in the underlying vector<>      
00460     * <li>[#size_type#]       - type for the number of elements in the vector<>
00461     * <li>[#difference_type#]   - type for the distance between elements
00462     * <li>[#iterator#]        - type for the forward non-constant iterator
00463     * <li>[#const_iterator#]      - type for the forward constant iterator (can't change elements)
00464     * <li>[#reverse_iterator#]    - type for the reverse non-constant iterator
00465     * <li>[#const_reverse_iterator#]  - type for the reverse constant iterator (can't change elements)
00466     * <li>[#reference#]       - type returned from subscripting, iterator deferencing etc.
00467     * <li>[#const_reference#]   - "" "" for const vector slice objects
00468     * \end{description}
00469     */
00470 
00471   // @{
00472   // @}
00473 
00474   typedef T                   value_type;
00475   typedef DenseLinAlgPack::size_type          size_type;
00476   typedef ptrdiff_t               difference_type;
00477   typedef value_type*               iterator;
00478   typedef const value_type*           const_iterator;
00479 #if 0 /* defined(_INTEL_CXX) || defined(_WINDOWS) */
00480   typedef std::reverse_iterator<iterator, value_type
00481     , value_type&, value_type*, difference_type>    reverse_iterator;
00482   typedef std::reverse_iterator<const_iterator
00483     , value_type, const value_type&, const value_type*
00484     , difference_type>                  const_reverse_iterator;
00485 #else
00486   typedef std::reverse_iterator<iterator>         reverse_iterator;
00487   typedef std::reverse_iterator<const_iterator>     const_reverse_iterator;
00488 #endif
00489   typedef value_type&               reference;
00490   typedef const value_type&           const_reference;
00491   typedef std::vector<value_type>         valarray;
00492 
00493   /* * @name {\bf Constructors}.
00494     * 
00495     * These constructors allocate and may initialize the elements of a 1-D vector.
00496     * The default C++ copy constructor is used and is therefore not show here.
00497     */
00498 
00499   // @{
00500 
00502   VectorTmpl();
00504   VectorTmpl(size_type n);
00506   VectorTmpl(value_type val, size_type n);
00508   /* * Constructs a vector with n elements and intializes elements to those of an array.
00509     *
00510     * Postconditions: <ul>
00511     *   <li> #this->operator[](i) == p[i]#, i = 0, 1, ... n
00512     *   </ul>
00513     */  
00514   VectorTmpl(const value_type* p, size_type n);
00516   /* * Constructs a DVector object fron a VectorSliceTmpl object.
00517     *
00518     * Postconditions: <ul>
00519     *   <li> #this->dim() == vs.dim()#
00520     *   <li> #this->operator[](i) == vs[i]#, i = 0, 1, ... n
00521     *   </ul>
00522     */  
00523   VectorTmpl(const VectorSliceTmpl<value_type>& vs);
00524 
00525   // @}
00526 
00527   /* * @name {\bf Memory Management / Misc}. */
00528 
00529   // @{
00530 
00532   /* * Resize the vector to hold n elements.
00533     *
00534     * Any new elements added are initialized to val.
00535     *
00536     * Postconditions: <ul>
00537     *   <li> #this->dim() == n#
00538     *   </ul>
00539     */  
00540   void resize(size_type n, value_type val = value_type());
00542   /* * Free memory and resize DVector to this->dim() == 0.
00543     *
00544     * Postconditions: <ul>
00545     *   <li> #this->dim() == 0#
00546     *   </ul>
00547     */  
00548   void free();
00550   size_type dim() const;
00552   /* * Returns the degree of memory overlap of this and the VectorSliceTmpl object vs.
00553     *
00554     * @return 
00555     *   \begin{description}
00556     *   <li>[NO_OVERLAP]  There is no memory overlap between this and vs
00557     *   <li>[SOME_OVERLAP]  There is some memory locations that this and vs share
00558     *   <li>[SAME_MEM]    The VectorSliceTmpl objects this and vs share the exact same memory locations.
00559     *   \end{description}
00560     */
00561   EOverLap overlap(const VectorSliceTmpl<value_type>& vs) const;
00563   operator VectorSliceTmpl<value_type>();
00565   operator const VectorSliceTmpl<value_type>() const;
00566 
00567   // @}
00568   
00569 
00570   /* * @name {\bf STL Iterator Access functions}.
00571     *
00572     * The iterators returned are valid STL random access iterators.
00573     * The forward iterators returned iterate from the first element to the last element.
00574     * The reverse iterators returned iterate from the last element to the first element.
00575     */
00576 
00577   // @{
00578 
00580   iterator begin();
00582   iterator end();
00584   const_iterator begin() const;
00586   const_iterator end() const;
00588   reverse_iterator rbegin();
00590   reverse_iterator rend();
00592   const_reverse_iterator rbegin() const;
00594   const_reverse_iterator rend() const;
00595 
00596   // @}
00597 
00598   /* * @name {\bf Individual Element Access Subscripting (lvalue)}.
00599     *
00600     * These operator functions allow access (lvalue) to the individual elements
00601     * of the DVector object.
00602     * 
00603     * The subscript i must be, 1 <= i <= this->dim(), for the 1-based element access
00604     * operators and, 0 <= i <= this->dim() - 1, for the 0-based element access operators.
00605     * If they are not then an #std::out_of_range# exception will be thrown.
00606     */
00607 
00608   // @{
00609 
00611   reference operator()(size_type i);
00613   const_reference operator()(size_type i) const;
00615   reference operator[](size_type i);
00617   const_reference operator[](size_type i) const;
00618 
00619   // @}
00620 
00621   /* * @name {\bf Subvector Access Operators}.
00622     *
00623     * These operator functions are used to create views of continous regions of the DVector.
00624     * Each of them returns a VectorSliceTmpl object for the region.  Constant (const) VectorSliceTmpl objects
00625     * are returned for a const DVector.  This means that the elements can not be changed
00626     * as should be the case.
00627     * 
00628     * Beware!  VC++ is returning non-const VectorSliceTmpl objects for the 
00629     * #VectorSliceTmpl operator()(...) const;# member functions and therefore a const \Ref{DVector} or
00630     * \Ref{VectorSliceTmpl} can be modifed my subsetting it.  Hopefully this problem will
00631     * be fixed in future versions of the compiler or I when will get another compiler.
00632     */
00633 
00634   // @{
00635 
00637   /* * Returns a VectorSliceTmpl object representing the entire DVector. 
00638     * 
00639     * Call this member function to force a type conversion to VectorSliceTmpl.  Using the
00640     * VectorSliceTmpl of a DVector for algebraic expressions used with the TCOL allows a for simplier
00641     * implementaion of those operations by cutting down on the number combinations.  This is
00642     * especialy true for longer optimized expression.
00643     */
00644   VectorSliceTmpl<value_type> operator()();
00646   const VectorSliceTmpl<value_type> operator()() const;
00648   /* * Returns a continous subregion of the DVector object.
00649     *
00650     * The returned VectorSliceTmpl object represents the range of the rng argument.
00651     *
00652     * Preconditions: <ul>
00653     *   <li> #rng.ubound() - 1 <= this->dim()# (throw #out_of_range#)
00654     *   </ul>
00655     *
00656     * @param  rng   Indece range [lbound,ubound] of the region being returned.
00657     */
00658   VectorSliceTmpl<value_type> operator()(const Range1D& rng);
00660   const VectorSliceTmpl<value_type> operator()(const Range1D& rng) const;
00662   /* * Returns a VectorSliceTmpl object for the continous subregion [ubound, lbound].
00663     * 
00664     * Preconditions: <ul>
00665     *   <li> #lbound > 1# (throw #out_of_range#)
00666     *   <li> #lbound < ubound# (throw #out_of_range#)
00667     *   <li> #ubound <= this->dim()# (throw #out_of_range#)
00668     *   </ul>
00669     *
00670     * @param  rng   Range [lbound,ubound] of the region being taken.
00671     */
00672   VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound);
00674   const VectorSliceTmpl<value_type> operator()(size_type lbound, size_type ubound) const;
00676   /* * Return a VectorSliceTmpl object the reverse of this DVector.
00677     *
00678     * In the reverse VectorSliceTmpl,
00679     * the first element becomes the last element and visa-versa.  For example, for 
00680     * #VectorSliceTmpl r = x.rev()#, #&x(1) == &z(z.dim())# and #&x(x.dim()) == &z(1)# are both true.
00681     * The iterators returned by \Ref{begin()} iterate from the first conceptual element to the last.
00682     */
00683   VectorSliceTmpl<value_type> rev();
00685   const VectorSliceTmpl<value_type> rev() const;
00686 
00687   // @}
00688 
00689   /* * @name {\bf Assignment Operators}. */
00690 
00691   // @{
00692 
00694   /* * vs = alpha (Sets all the elements to the constant alpha).
00695     *
00696     * Preconditions: <ul>
00697     *   <li> #this->dim() > 0# (throw #std::length_error#)
00698     *   </ul>
00699     *
00700     * Postconditions: <ul>
00701     *   <li> #this->operator()(i) == alpha#, i = 1, 2, ... , #this->dim()#
00702     *   </ul>
00703     */
00704   VectorTmpl<value_type>& operator=(value_type alpha);
00706   /* * vs = rhs (Copies the elements of rhs into the elements of this).
00707     *
00708     * Preconditions: <ul>
00709     *   <li> #this->dim() == rhs.dim()# (throw #out_of_range#)
00710     *   <li> #rhs.dim() > 0# (throw #out_of_range#)
00711     *   </ul>
00712     *
00713     * Postconditions: <ul>
00714     *   <li> #this->operator()(i) == rhs(i)#, i = 1, 2, ..., #this->dim()#
00715     *   </ul>
00716     */
00717   VectorTmpl<value_type>& operator=(const VectorSliceTmpl<value_type>& rhs);
00719   /* * Needed to override the default assignment operator.
00720     */
00721   VectorTmpl<value_type>& operator=(const VectorTmpl<value_type>& rhs);
00722 
00723   // @}
00724 
00725   /* * @name {\bf Implementation Access}.
00726     *
00727     * Provides access to underlying raw data.
00728     */
00729 
00730   // @{
00731 
00733   value_type*     raw_ptr();
00735   const value_type* raw_ptr() const;
00737   value_type*     start_ptr();
00739   const value_type* start_ptr() const;
00741   difference_type   stride() const;
00742 
00743   // @}
00744   
00745 private:
00746   valarray v_;
00747 
00748 }; // end class VectorTmpl<T>
00749 
00750 //    end DVector Classes scope
00751 // @}
00752 
00753 // ///////////////////////////////////////////////////////////////////////////////
00754 // Non-member function declarations                       //
00755 // ///////////////////////////////////////////////////////////////////////////////
00756 
00757 
00758 /* * @name {\bf Non-Member Functions}. */
00759 
00760 // @{
00761 //    begin non-member functions scope
00762 
00763 //
00764 size_type vector_validate_sized(size_type size);
00765 //
00766 void vector_validate_range(size_type ubound, size_type max_ubound);
00767 //
00768 void vector_validate_subscript(size_type size, size_type i);
00770 /* * Utility for checking the sizes of two VectorSliceTmpl objects and throwing an exception
00771   * if the sizes are not the same.
00772   */
00773 void assert_vs_sizes(size_type size1, size_type size2);
00774 
00776 /* * Create a general vector slice.
00777   */
00778 template<class T>
00779 inline
00780 VectorSliceTmpl<T> gen_vs( VectorSliceTmpl<T>& vs, size_type start, size_type size, ptrdiff_t stride )
00781 {
00782   return VectorSliceTmpl<T>( vs.start_ptr() + vs.stride() * (start-1), size, vs.stride() * stride );
00783 }
00784 
00786 template<class T>
00787 inline
00788 const VectorSliceTmpl<T> gen_vs( const VectorSliceTmpl<T>& vs, size_type start, size_type size
00789   , ptrdiff_t stride )
00790 {
00791   return VectorSliceTmpl<T>( const_cast<typename VectorSliceTmpl<T>::value_type*>(vs.start_ptr()) + vs.stride() * (start-1)
00792     , size, vs.stride() * stride );
00793 }
00794 
00795 //    end non-member functions scope
00796 // @}
00797 
00798 //    end Vectors scope
00799 // @}
00800 
00801 } // end namespace DenseLinAlgPack
00802 
00803 // ////////////////////////////////////////////////////////////////////////////////
00804 // Inline definitions of member function definitions               //
00805 // ////////////////////////////////////////////////////////////////////////////////
00806 
00807 // ////////////////////////////////////////////////////////////////////////
00808 // Non-member functions / utilities
00809 
00810 #ifndef LINALGPACK_CHECK_SLICE_SETUP
00811 inline
00812 DenseLinAlgPack::size_type DenseLinAlgPack::vector_validate_sized(size_type size)
00813 {
00814   return size;
00815 }
00816 #endif
00817 
00818 #ifndef LINALGPACK_CHECK_RANGE
00819 inline
00820 void DenseLinAlgPack::vector_validate_range(size_type ubound, size_type max_ubound)
00821 {}
00822 #endif
00823 
00824 #ifndef LINALGPACK_CHECK_RANGE
00825 inline
00826 void DenseLinAlgPack::vector_validate_subscript(size_type size, size_type i)
00827 {}
00828 #endif
00829 
00830 #ifndef LINALGPACK_CHECK_RHS_SIZES
00831 inline
00832 void DenseLinAlgPack::assert_vs_sizes(size_type size1, size_type size2)
00833 {}
00834 #endif
00835 
00836 namespace DenseLinAlgPack {
00837 
00838 // /////////////////////////////////////////////////////////////////////////////
00839 // VectorSliceTmpl inline member function definitions
00840 
00841 // Constructors.  Use default copy constructor
00842 
00843 template<class T>
00844 inline
00845 VectorSliceTmpl<T>::VectorSliceTmpl()
00846   : ptr_(0)
00847   , size_(0)
00848   , stride_(0)
00849 {}
00850 
00851 template<class T>
00852 inline
00853 VectorSliceTmpl<T>::VectorSliceTmpl( value_type* ptr, size_type size, difference_type stride)
00854   : ptr_(ptr)
00855   , size_(size)
00856   , stride_(stride)
00857 {}
00858 
00859 template<class T>
00860 inline
00861 VectorSliceTmpl<T>::VectorSliceTmpl( value_type* ptr, size_type size, const Range1D& rng )
00862   : ptr_( ptr + rng.lbound() - 1 )
00863   , size_( rng.full_range() ? vector_validate_sized(size) : rng.size() )
00864   , stride_(1)
00865 {
00866   vector_validate_range( rng.full_range() ? size : rng.ubound(), size );
00867 }
00868 
00869 template<class T>
00870 inline
00871 VectorSliceTmpl<T>::VectorSliceTmpl( VectorSliceTmpl<T>& vs, const Range1D& rng )
00872   : ptr_( vs.start_ptr() + (rng.lbound() - 1) * vs.stride() )
00873   , size_( rng.full_range() ? vector_validate_sized(vs.dim()) : rng.size() )
00874   , stride_( vs.stride() )
00875 { vector_validate_range(  rng.full_range() ? vs.dim() : rng.ubound(), vs.dim() ); }
00876 
00877 template<class T>
00878 inline
00879 void VectorSliceTmpl<T>::bind(VectorSliceTmpl vs)
00880 {
00881   ptr_  = vs.ptr_;
00882   size_ = vs.size_;
00883   stride_ = vs.stride_;
00884 }
00885 
00886 // Iterator functions
00887 template<class T>
00888 inline
00889 typename VectorSliceTmpl<T>::iterator VectorSliceTmpl<T>::begin()
00890 { return iterator(start_ptr(), stride()); }
00891 
00892 template<class T>
00893 inline
00894 typename VectorSliceTmpl<T>::iterator VectorSliceTmpl<T>::end()
00895 { return iterator(start_ptr() + dim() * stride(), stride()); }
00896 
00897 template<class T>
00898 inline
00899 typename VectorSliceTmpl<T>::const_iterator VectorSliceTmpl<T>::begin() const
00900 { return const_iterator(start_ptr(), stride()); }
00901 
00902 template<class T>
00903 inline
00904 typename VectorSliceTmpl<T>::const_iterator VectorSliceTmpl<T>::end() const
00905 { return const_iterator(start_ptr() + dim() * stride(), stride()); }
00906 
00907 template<class T>
00908 inline
00909 typename VectorSliceTmpl<T>::reverse_iterator VectorSliceTmpl<T>::rbegin()
00910 { return reverse_iterator(end()); }
00911 
00912 template<class T>
00913 inline
00914 typename VectorSliceTmpl<T>::reverse_iterator VectorSliceTmpl<T>::rend()
00915 { return reverse_iterator(begin()); }
00916 
00917 template<class T>
00918 inline
00919 typename VectorSliceTmpl<T>::const_reverse_iterator VectorSliceTmpl<T>::rbegin() const
00920 { return const_reverse_iterator(end()); }
00921 
00922 template<class T>
00923 inline
00924 typename VectorSliceTmpl<T>::const_reverse_iterator VectorSliceTmpl<T>::rend() const
00925 { return const_reverse_iterator(begin()); }
00926 
00927 // Element access
00928 template<class T>
00929 inline
00930 typename VectorSliceTmpl<T>::reference VectorSliceTmpl<T>::operator()(size_type i) // 1 based
00931 {
00932   vector_validate_subscript(dim(),i);
00933   return ptr_[(i-1)*stride_];
00934 }
00935 
00936 template<class T>
00937 inline
00938 typename VectorSliceTmpl<T>::const_reference VectorSliceTmpl<T>::operator()(size_type i) const
00939 {
00940   vector_validate_subscript(dim(),i);
00941   return ptr_[(i-1)*stride_];
00942 }
00943 
00944 template<class T>
00945 inline
00946 typename VectorSliceTmpl<T>::reference VectorSliceTmpl<T>::operator[](size_type i) // 0 based   
00947 {
00948   vector_validate_subscript(dim(),i+1);
00949   return ptr_[(i)*stride_];
00950 }
00951 
00952 template<class T>
00953 inline
00954 typename VectorSliceTmpl<T>::const_reference VectorSliceTmpl<T>::operator[](size_type i) const
00955 {
00956   vector_validate_subscript(dim(),i+1);
00957   return ptr_[(i)*stride_];
00958 }
00959 
00960 // Subregion Access.  Let the constructors of VectorSliceTmpl validate the ranges
00961 template<class T>
00962 inline
00963 VectorSliceTmpl<T>& VectorSliceTmpl<T>::operator()() 
00964 { return *this; }
00965 
00966 template<class T>
00967 inline
00968 const VectorSliceTmpl<T>& VectorSliceTmpl<T>::operator()() const 
00969 { return *this; }
00970 
00971 template<class T>
00972 inline
00973 VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(const Range1D& rng) 
00974 { return VectorSliceTmpl(*this, RangePack::full_range(rng,1,dim())); }
00975 
00976 template<class T>
00977 inline
00978 const VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(const Range1D& rng) const
00979 { return VectorSliceTmpl(const_cast<VectorSliceTmpl<T>&>(*this), RangePack::full_range(rng,1,dim())); }
00980 
00981 template<class T>
00982 inline
00983 VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(size_type lbound, size_type ubound)
00984 { return VectorSliceTmpl(*this, Range1D(lbound, ubound)); }
00985 
00986 template<class T>
00987 inline
00988 const VectorSliceTmpl<T> VectorSliceTmpl<T>::operator()(size_type lbound, size_type ubound) const
00989 { return VectorSliceTmpl(const_cast<VectorSliceTmpl<T>&>(*this), Range1D(lbound, ubound)); }
00990 
00991 template<class T>
00992 inline
00993 VectorSliceTmpl<T> VectorSliceTmpl<T>::rev()
00994 { return VectorSliceTmpl( start_ptr() + stride() * (dim()-1), dim(), - stride() ); }
00995 
00996 template<class T>
00997 inline
00998 const VectorSliceTmpl<T> VectorSliceTmpl<T>::rev() const
00999 { return VectorSliceTmpl( const_cast<value_type*>(start_ptr()) + stride() * (dim()-1), dim(), - stride() ); }
01000 
01001 // Assignment Operators
01002 template<class T>
01003 inline
01004 VectorSliceTmpl<T>& VectorSliceTmpl<T>::operator=(value_type alpha)
01005 {
01006   std::fill(begin(),end(),alpha);
01007   return *this;
01008 }
01009 
01010 template<class T>
01011 inline
01012 VectorSliceTmpl<T>& VectorSliceTmpl<T>::operator=(const VectorSliceTmpl<T>& rhs) 
01013 {
01014   assert_vs_sizes(this->dim(),rhs.dim());
01015   std::copy(rhs.begin(),rhs.end(),begin());
01016   return *this;
01017 }
01018 
01019 // Misc. member functions
01020 
01021 template<class T>
01022 inline
01023 typename VectorSliceTmpl<T>::size_type VectorSliceTmpl<T>::dim() const
01024 { return size_; }
01025 
01026 // Raw pointer access
01027 
01028 template<class T>
01029 inline
01030 typename VectorSliceTmpl<T>::value_type*  VectorSliceTmpl<T>::raw_ptr()
01031 { return stride() > 0 ? start_ptr() : start_ptr() + stride() * (dim() - 1); }
01032 
01033 template<class T>
01034 inline
01035 const typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::raw_ptr() const
01036 { return stride() > 0 ? start_ptr() : start_ptr() + stride() * (dim() - 1); }
01037 
01038 template<class T>
01039 inline
01040 typename VectorSliceTmpl<T>::value_type*  VectorSliceTmpl<T>::start_ptr()
01041 { return ptr_; }
01042 
01043 template<class T>
01044 inline
01045 const typename VectorSliceTmpl<T>::value_type* VectorSliceTmpl<T>::start_ptr() const
01046 { return ptr_; }
01047 
01048 template<class T>
01049 inline
01050 typename VectorSliceTmpl<T>::difference_type VectorSliceTmpl<T>::stride() const
01051 { return stride_; }
01052 
01053 // /////////////////////////////////////////////////////////////////////////////
01054 // DVector inline member function definitions
01055 
01056 // Constructors
01057 template<class T>
01058 inline
01059 VectorTmpl<T>::VectorTmpl()
01060 {}  // used to shut satisfy compiler
01061 
01062 template<class T>
01063 inline
01064 VectorTmpl<T>::VectorTmpl(size_type n) 
01065   : v_(n)
01066 {}
01067 
01068 template<class T>
01069 inline
01070 VectorTmpl<T>::VectorTmpl(value_type val, size_type n) 
01071   : v_(n)
01072 {
01073   std::fill(begin(),end(),val);
01074 }
01075 
01076 template<class T>
01077 inline
01078 VectorTmpl<T>::VectorTmpl(const value_type* p, size_type n)
01079   : v_(n)
01080 {
01081   std::copy(p,p+n,begin());
01082 }
01083 
01084 template<class T>
01085 inline
01086 VectorTmpl<T>::VectorTmpl(const VectorSliceTmpl<T>& vs)
01087   : v_(vs.dim())
01088 {  
01089   std::copy(vs.begin(),vs.end(),begin());
01090 }
01091 
01092 // Memory management
01093 template<class T>
01094 inline
01095 void VectorTmpl<T>::resize(size_type n, value_type val)
01096 {
01097   v_.resize(n);
01098   std::fill(begin(),end(),val);
01099  }
01100 
01101 template<class T>
01102 inline
01103 void VectorTmpl<T>::free()
01104 {
01105   v_.resize(0);
01106 }
01107 
01108 // Size
01109 template<class T>
01110 inline
01111 typename VectorTmpl<T>::size_type VectorTmpl<T>::dim() const
01112 { return v_.size(); }
01113 
01114 // Iterator functions
01115 template<class T>
01116 inline
01117 typename VectorTmpl<T>::iterator VectorTmpl<T>::begin()
01118 { return start_ptr(); }
01119 
01120 template<class T>
01121 inline
01122 typename VectorTmpl<T>::iterator VectorTmpl<T>::end()
01123 { return start_ptr() + dim(); }
01124 
01125 template<class T>
01126 inline
01127 typename VectorTmpl<T>::const_iterator VectorTmpl<T>::begin() const
01128 { return start_ptr(); }
01129 
01130 template<class T>
01131 inline
01132 typename VectorTmpl<T>::const_iterator VectorTmpl<T>::end() const 
01133 { return start_ptr() + dim(); }
01134 
01135 template<class T>
01136 inline
01137 typename VectorTmpl<T>::reverse_iterator  VectorTmpl<T>::rbegin()
01138 { return reverse_iterator(end()); }
01139 
01140 template<class T>
01141 inline
01142 typename VectorTmpl<T>::reverse_iterator  VectorTmpl<T>::rend()
01143 { return reverse_iterator(begin()); }
01144 
01145 template<class T>
01146 inline
01147 typename VectorTmpl<T>::const_reverse_iterator VectorTmpl<T>::rbegin() const
01148 { return const_reverse_iterator(end()); }
01149 
01150 template<class T>
01151 inline
01152 typename VectorTmpl<T>::const_reverse_iterator VectorTmpl<T>::rend() const 
01153 { return const_reverse_iterator(begin()); }
01154 
01155 // Element access
01156 template<class T>
01157 inline
01158 typename VectorTmpl<T>::reference VectorTmpl<T>::operator()(size_type i)
01159 {
01160   vector_validate_subscript(dim(),i);
01161   return start_ptr()[i-1];
01162 }
01163 
01164 template<class T>
01165 inline
01166 typename VectorTmpl<T>::const_reference VectorTmpl<T>::operator()(size_type i) const
01167 {
01168   vector_validate_subscript(dim(),i);
01169   return start_ptr()[i-1];
01170 }
01171 
01172 template<class T>
01173 inline
01174 typename VectorTmpl<T>::reference VectorTmpl<T>::operator[](size_type i)
01175 {
01176   vector_validate_subscript(dim(),i+1);
01177   return start_ptr()[i];
01178 }
01179 
01180 template<class T>
01181 inline
01182 typename VectorTmpl<T>::const_reference VectorTmpl<T>::operator[](size_type i) const
01183 {
01184   vector_validate_subscript(dim(),i+1);
01185   return start_ptr()[i];
01186 }
01187 
01188 // Subregion Access.  Leave validation to the VectorSliceTmpl constructors.
01189 template<class T>
01190 inline
01191 VectorSliceTmpl<T> VectorTmpl<T>::operator()() 
01192 { return VectorSliceTmpl<T>(start_ptr(),dim()); }
01193 
01194 template<class T>
01195 inline
01196 const VectorSliceTmpl<T> VectorTmpl<T>::operator()() const
01197 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()),dim()); }
01198 
01199 template<class T>
01200 inline
01201 VectorSliceTmpl<T> VectorTmpl<T>::operator()(const Range1D& rng) 
01202 { return VectorSliceTmpl<T>(start_ptr(),dim(),rng); }
01203 
01204 template<class T>
01205 inline
01206 const VectorSliceTmpl<T> VectorTmpl<T>::operator()(const Range1D& rng) const
01207 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()),dim(),rng); }
01208 
01209 template<class T>
01210 inline
01211 VectorSliceTmpl<T> VectorTmpl<T>::operator()(size_type lbound, size_type ubound)
01212 { return VectorSliceTmpl<T>(start_ptr(), dim(), Range1D(lbound, ubound)); }
01213 
01214 template<class T>
01215 inline
01216 const VectorSliceTmpl<T> VectorTmpl<T>::operator()(size_type lbound, size_type ubound) const
01217 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()), dim(), Range1D(lbound, ubound)); }
01218 
01219 template<class T>
01220 inline
01221 VectorSliceTmpl<T> VectorTmpl<T>::rev()
01222 { return VectorSliceTmpl<T>( start_ptr() + dim() - 1, dim(), -1 ); }
01223 
01224 template<class T>
01225 inline
01226 const VectorSliceTmpl<T> VectorTmpl<T>::rev() const
01227 { return VectorSliceTmpl<T>( const_cast<value_type*>(start_ptr()) + dim() - 1, dim(), -1 ); }
01228 
01229 // Conversion operators
01230 template<class T>
01231 inline
01232 VectorTmpl<T>::operator VectorSliceTmpl<T>()
01233 { return VectorSliceTmpl<T>(start_ptr(), dim()); }
01234 
01235 template<class T>
01236 inline
01237 VectorTmpl<T>::operator const VectorSliceTmpl<T>() const
01238 { return VectorSliceTmpl<T>(const_cast<value_type*>(start_ptr()), dim()); }
01239 
01240 // Assignment Operators
01241 template<class T>
01242 inline
01243 VectorTmpl<T>& VectorTmpl<T>::operator=(value_type alpha) 
01244 {
01245   if(!dim()) resize(1);
01246   std::fill(begin(),end(),alpha);
01247   return *this;
01248 }
01249 
01250 template<class T>
01251 inline
01252 VectorTmpl<T>& VectorTmpl<T>::operator=(const VectorTmpl<T>& rhs) 
01253 {
01254   resize(rhs.dim());
01255   std::copy(rhs.begin(),rhs.end(),begin());
01256   return *this;
01257 }
01258 
01259 template<class T>
01260 inline
01261 VectorTmpl<T>& VectorTmpl<T>::operator=(const VectorSliceTmpl<T>& rhs) 
01262 {
01263   resize(rhs.dim());
01264   std::copy(rhs.begin(),rhs.end(),begin());
01265   return *this;
01266 }
01267 
01268 // Raw pointer access
01269 
01270 template<class T>
01271 inline
01272 typename VectorTmpl<T>::value_type* VectorTmpl<T>::raw_ptr()
01273 { return start_ptr(); }
01274 
01275 template<class T>
01276 inline
01277 const typename VectorTmpl<T>::value_type* VectorTmpl<T>::raw_ptr() const
01278 { return start_ptr(); }
01279 
01280 template<class T>
01281 inline
01282 typename VectorTmpl<T>::value_type* VectorTmpl<T>::start_ptr()
01283 { return dim() ? &(v_)[0] : 0; }
01284 
01285 template<class T>
01286 inline
01287 const typename VectorTmpl<T>::value_type* VectorTmpl<T>::start_ptr() const
01288 { return &const_cast<valarray&>((v_))[0]; }
01289 
01290 template<class T>
01291 inline
01292 typename VectorTmpl<T>::difference_type VectorTmpl<T>::stride() const
01293 { return 1; }
01294 
01295 // //////////////////////////////////////////////////
01296 // Non-inlined members
01297 
01298 // Assume that the vector slices are the rows, cols or diag of a 2-D matrix.
01299 template<class T>
01300 EOverLap VectorSliceTmpl<T>::overlap(const VectorSliceTmpl<value_type>& vs) const {
01301 
01302   const typename VectorSliceTmpl<T>::value_type
01303     *raw_ptr1 = ( stride() > 0 ? start_ptr() : start_ptr() + (dim()-1)*stride() ),
01304     *raw_ptr2 = ( vs.stride() > 0 ? vs.start_ptr() : vs.start_ptr() + (vs.dim()-1)*vs.stride() );
01305   typename VectorSliceTmpl<T>::size_type
01306     size1 = dim(),
01307     size2 = vs.dim();
01308   typename VectorSliceTmpl<T>::difference_type
01309     stride1 = std::abs(stride()),
01310     stride2 = std::abs(vs.stride());
01311 
01312   // Establish a frame of reference where raw_ptr1 < raw_ptr2
01313   if(raw_ptr1 > raw_ptr2) {
01314     std::swap(raw_ptr1,raw_ptr2);
01315     std::swap(stride1,stride2);
01316     std::swap(size1,size2);
01317   }
01318 
01319   if( raw_ptr1 + stride1 * (size1 - 1) < raw_ptr2 ) {
01320     return NO_OVERLAP; // can't be any overlap
01321   }
01322 
01323   typename VectorSliceTmpl<T>::size_type
01324     start1 = 0,
01325     start2 = raw_ptr2 - raw_ptr1;
01326 
01327   if(start1 == start2 && stride1 == stride2 && size1 == size2)
01328     return SAME_MEM;
01329 //  else if(start1 == start2)
01330 //    return SOME_OVERLAP;  // First elements are the same
01331 //  else if(stride1 + (size1 - 1) * stride1 == stride2 + (size2 - 1) * stride2)
01332 //    return SOME_OVERLAP;  // Last elements are the same
01333   else if(stride1 == stride2) {
01334     if(!((start2 - start1) % stride1))
01335       return SOME_OVERLAP;
01336     else
01337       return NO_OVERLAP; // a different row, col or diag of a matrix
01338   }
01339   else {
01340     if(stride1 == 1 || stride2 == 1) {
01341       // One of them is a column vector.
01342       // Make vs1 the column vector.
01343       bool switch_them = (stride2 == 1);
01344       if(switch_them) {
01345         std::swap(start1,start2);
01346         std::swap(stride1,stride2);
01347         std::swap(size1,size2);
01348       }
01349 
01350       // Determine if the other vector could be row vector
01351       // or must be a diag vector.  If using stride2 makes
01352       // the first and last elements of the column vector
01353       // on different rows, then the other vector must be a diagonal.
01354       // col_first = start1/stride2, col_last = (start1+size1-1)/stride2
01355       // if(col_last - col_first > 0) then vs2 must be a diagonal vector
01356       // with max_rows = stride2 - 1.
01357       size_t max_rows = (start1+size1-1)/stride2 - start1/stride2 > 0 ? stride2 - 1 : stride2;
01358 
01359       // find the index (0-based) of vs2 that intersects this column.
01360       size_t vs2_col_i = start1/max_rows - start2/max_rows;
01361       
01362       // See if the first element of the column is above the element in vs2
01363       // and the last element of the column is below the element.  If it is
01364       // then we conclude that there is an itersection.
01365       size_t vs2_col_rng = start2 + vs2_col_i * stride2;
01366       if(start1 <= vs2_col_rng && vs2_col_rng <= start1+size1-1)
01367         return SOME_OVERLAP;
01368       else
01369         return NO_OVERLAP;
01370     }
01371     // They are not the same and nether is a column vector so one is a row vector
01372     // and the other is a diagonal vector.
01373     // Nether is a column vector so choose as vs1 the row vector (the smaller stride).
01374     bool switch_them = stride2 < stride1;
01375     if(switch_them) {
01376       std::swap(start1,start2);
01377       std::swap(stride1,stride2);
01378       std::swap(size1,size2);
01379     }
01380 
01381     size_t max_rows = stride1;
01382     // Determine the first and last columns (0-based) in the original
01383     // matrix where there vs1 and vs2 intersect.
01384     size_t  sec_first_col = (start1 > start2) ? start1/max_rows : start2/max_rows,
01385         last1 = start1 + (size1 - 1) * stride1,
01386         last2 = start2 + (size2 - 1) * stride2,
01387         sec_last_col = (last1 < last2) ? last1/max_rows : last2/max_rows;
01388     // Determine the vector indexes (0-based) of vs1 and vs2 for the start and end
01389     // in this region
01390     size_t  vs1_first_col = start1 / max_rows,
01391         vs2_first_col = start2 / max_rows;
01392     
01393     // Determine the indexes in the valarray of the two vectors for the two ends
01394     size_t  vs1_first_col_i = sec_first_col - vs1_first_col,
01395         vs1_last_col_i = sec_last_col - vs1_first_col,
01396         vs2_first_col_i = sec_first_col - vs2_first_col,
01397         vs2_last_col_i = sec_last_col - vs2_first_col;
01398 
01399     // Compare the indexes in the valarray at the two ends.  If they cross then
01400     // there must be an element of overlap.  Uses equivalent of the intermediate
01401     // value therorm.
01402     // Must cast to an int that can hold a negative value
01403     ptrdiff_t diff1 = (start1 + vs1_first_col_i * stride1) 
01404             - static_cast<ptrdiff_t>((start2 + vs2_first_col_i * stride2)),
01405           diff2 = (start1 + vs1_last_col_i * stride1) 
01406             - static_cast<ptrdiff_t>((start2 + vs2_last_col_i * stride2));
01407     if(diff1 * diff2 > 0 )
01408       return NO_OVERLAP;    // they do not cross
01409     else
01410       return SOME_OVERLAP;  // they share an element
01411   }
01412 }
01413 
01414 template<class T>
01415 EOverLap VectorTmpl<T>::overlap(const VectorSliceTmpl<value_type>& vs) const {
01416 
01417   const typename VectorSliceTmpl<T>::value_type
01418     *raw_ptr1 = ( stride() > 0 ? start_ptr() : start_ptr() + (dim()-1)*stride() ),
01419     *raw_ptr2 = ( vs.stride() > 0 ? vs.start_ptr() : vs.start_ptr() + (vs.dim()-1)*vs.stride() );
01420   typename VectorSliceTmpl<T>::size_type
01421     size1 = dim(),
01422     size2 = vs.dim();
01423 
01424   if( raw_ptr1 <= raw_ptr2 && raw_ptr2 + size2 <= raw_ptr1 + size1 ) {
01425     if( raw_ptr1 == raw_ptr2 && size1 == size2 && 1 == vs.stride() )
01426       return SAME_MEM;
01427     return SOME_OVERLAP;
01428   }
01429   return NO_OVERLAP;
01430 }
01431 
01432 } // end namespace DenseLinAlgPack
01433 
01434 #endif  // end VECTOR_CLASS_TMPL_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines