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

Generated on Tue Jul 13 09:30:52 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7