MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_DMatrixClass.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 // General matrix and matrix region (slice) classes
00043 
00044 #ifndef GEN_MATRIX_CLASS_H
00045 #define GEN_MATRIX_CLASS_H
00046 
00047 #include "DenseLinAlgPack_DVectorClass.hpp"
00048 #include "DenseLinAlgPack_DMatrixAssign.hpp"
00049 
00050 /* * @name {\bf Dense 2-D Rectangular Matrix Absractions}.
00051   *
00052   * The class DMatrix is a storage class for 2-D matrices while the class DMatrixSlice
00053   * is used to represent rectangular regions of a DMatrix object.
00054   */
00055 // @{
00056 //    begin General Rectangular 2-D Matrices scope
00057 
00058 namespace DenseLinAlgPack {
00059 
00060 class DMatrix;
00061 
00062 /* * @name {\bf General Matrix Classes}. */
00063 // @{
00064 
00065 // ////////////////////////////////////////////////////////////////////////////////////////////////////////
00066 // GenMatrixClass
00067 //
00068 
00070 /* * 2-D General Rectangular Matrix Subregion Class (Slice) (column major).
00071   *
00072   * This class is used to represent a rectangular matrix.  It uses a BLAS-like
00073   * slice of a raw C++ array.  Objects of this class can represent
00074   * an entire matrix or any rectangular subregion of a matrix.
00075   */
00076 class DMatrixSlice {
00077 public:
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 valarray<>      
00084     * <li>[#size_type#]       - type for the number rows and coluns
00085     * <li>[#difference_type#]   - type for the stride between elements (in a row or column)
00086     * <li>[#reference#]       - value_type&
00087     * <li>[#const_reference#]   - const value_type&
00088     * \end{description}
00089     */
00090   // @{
00091   // @}
00092   typedef DenseLinAlgPack::value_type         value_type;
00093   typedef DenseLinAlgPack::size_type          size_type;
00094   typedef ptrdiff_t               difference_type;
00095   typedef value_type&               reference;
00096   typedef const value_type&           const_reference;
00097 
00098   /* * @name {\bf Constructors}.
00099     *
00100     * These constructors are used by the other entities in the library
00101     * to create DMatrixSlices.  In general user need not use these
00102     * constructors directly.  Instead, the general user should use the 
00103     * subscript operators to create subregions of DMatrix and DMatrixSlice
00104     * objects.
00105     * 
00106     * The default copy constructor is used and is therefore not shown here.
00107     */
00108 
00109   // @{
00110 
00112   /* * Construct an empty view.
00113     *
00114     * The client can then call bind(...) to bind the view.
00115     */
00116   DMatrixSlice();
00117 
00119   /* * Construct a veiw of a matrix from a raw C++ array.
00120     *
00121     * The DMatrixSlice constructed represents a 2-D matrix whos elements are stored
00122     * in the array starting at #ptr#.  This is how the BLAS represent general rectangular
00123     * matrices.
00124     * The class can be used to provide a non-constant view the elements (#DMatrix#)
00125     * or a constant view (#const DMatrixSlice#).  Here is an example of how to
00126     * create a constant view:
00127     *
00128     \verbatim
00129     const DMatrixSlice::size_type m = 4, n = 4;
00130     const DMatrixSlice::value_type ptr[m*n] = { ... };
00131     const GenMatrixslice mat(cosnt_cast<DMatrixSlice::value_type*>(ptr),m*n,m,m,n);
00132     \endverbatim
00133     *
00134     * The #const_cast<...># such as in the example above is perfectly safe to use
00135     * when constructing  #const# veiw of #const# data.  On the other hand casting
00136     * away #const# and then non-#const# use is not safe in general since the original
00137     * #const# data may reside in ROM for example.  By using a non-#const# pointer in the
00138     * constructor you avoid accidentally making a non-#const# view of #const# data.
00139     *
00140     * Preconditions: <ul>
00141     *   <li> #rows <= max_rows# (throw out_of_range)
00142     *   <li> #start + rows + max_rows * (cols - 1) <= v.size()# (throw out_of_range)
00143     *   </ul>
00144     *
00145     * @param  ptr     pointer to the storage of the elements of the matrix (column oriented).
00146     * @param  size    total size of the storage pointed to by #ptr# (for size checking)
00147     * @param  max_rows  number of rows in the full matrix this sub-matrix was taken from.
00148     *           This is equivalent to the leading dimension argument (LDA) in
00149     *           the BLAS.
00150     * @param  rows    number of rows this matrix represents
00151     * @param  cols    number of columns this matrix represents
00152     */
00153   DMatrixSlice( value_type* ptr, size_type size
00154     , size_type max_rows, size_type rows, size_type cols );
00155 
00157   /* * Construct a submatrix region of another DMatrixSlice.
00158     *
00159     * This constructor simplifies the creation of subregions using the subscript
00160     * operators.
00161     *
00162     * I and J must be bounded ranges (full_range() == false).
00163     *
00164     * Preconditions: <ul>
00165     *   <li> #I.full_range() == false# (throw out_of_range)
00166     *   <li> #J.full_range() == false# (throw out_of_range)
00167     *   <li> #I.ubound() <= gms.rows()# (throw out_of_range)
00168     *   <li> #J.ubound() <= gms.cols()# (throw out_of_range)
00169     *   </ul>
00170     */
00171   DMatrixSlice( DMatrixSlice& gms, const Range1D& I
00172     , const Range1D& J );
00173 
00174   // @}
00175 
00177   /* * Set to the view of the input DMatrixSlice.
00178     *
00179     *
00180     */
00181   void bind( DMatrixSlice gms );
00182 
00183   /* * @name {\bf Dimensionality, Misc}. */
00184   // @{
00185 
00187   size_type   rows() const;
00189   size_type   cols() const;
00190 
00192   /* * Returns the degree of memory overlap of #this# and the DMatrixSlice object #gms#.
00193     *
00194     * @return 
00195     *   \begin{description}
00196     *   <li>[NO_OVERLAP]  There is no memory overlap between #this# and #gms#.
00197     *   <li>[SOME_OVERLAP]  There is some memory locations that #this# and #gms# share.
00198     *   <li>[SAME_MEM]    The DMatrixSlice objects #this# and #gms# share the exact same memory locations.
00199     *   \end{description}
00200     */
00201   EOverLap overlap(const DMatrixSlice& gms) const;
00202 
00203   // @}
00204 
00205   /* * @name {\bf Individual Element Access Subscripting (lvalue)}. */
00206   // @{
00207 
00209   reference   operator()(size_type i, size_type j);
00211   const_reference operator()(size_type i, size_type j) const;
00212 
00213   // @}
00214 
00215   /* * @name {\bf Subregion Access (1-based)}.
00216     *
00217     * These member functions allow access to subregions of the DMatrixSlice object.
00218     * The functions, row(i), col(j), and diag(k) return DVectorSlice objects while
00219     * the subscripting operators opeator()(I,J) return DMatrixSlice objects for
00220     * rectangular subregions.
00221     */
00222   // @{
00223 
00225   DVectorSlice      row(size_type i);
00227   const DVectorSlice  row(size_type i) const;
00229   DVectorSlice      col(size_type j);
00231   const DVectorSlice  col(size_type j) const;
00233   /* * Return DVectorSlice object representing a diagonal.
00234     *
00235     * Passing k == 0 returns the center diagonal.  Values of k < 0 are the lower diagonals
00236     * (k = -1, -2, ..., -#this->rows()# + 1).  Values of k > 0 are the upper diagonals
00237     * (k = 1, 2, ..., #this->cols()# - 1).
00238     *
00239     * Preconditions: <ul>
00240     *   <li> #[k < 0] k <= this->rows() + 1# (throw out_of_range)
00241     *   <li> #[k > 0] k <= this->cols() + 1# (throw out_of_range)
00242     *   </ul>
00243     */
00244   DVectorSlice      diag(difference_type k = 0);
00246   const DVectorSlice  diag(difference_type k = 0) const;
00248   /* * Extract a rectangular subregion containing rows I, and columns J.
00249     *
00250     * This operator function returns a DMatrixSlice that represents a
00251     * rectangular region of this DMatrixSlice.  This submatrix region
00252     * represents the rows I.lbound() to I.ubound() and columns J.lbound()
00253     * to J.lbound().  If I or J is unbounded (full_range() == true, constructed
00254     * with Range1D()), the all of the rows or columns respectively will be
00255     * selected.  For example. To select all the rows and the first 5 columns of
00256     * a matrix #A# you would use #A(Range1D(),Range1D(1,5))#.
00257     *
00258     * Preconditions: <ul>
00259     *   <li> #[I.full_range() == false] I.ubound() <= this->rows()# (throw out_of_range)
00260     *   <li> #[J.full_range() == false] J.ubound() <= this->cols()# (throw out_of_range)
00261     *   </ul>
00262     */
00263   DMatrixSlice operator()(const Range1D& I, const Range1D& J);
00265   const DMatrixSlice operator()(const Range1D& I, const Range1D& J) const;
00267   /* * Extract a rectangular subregion containing rows i1 to i2, and columns j1 to j2.
00268     *
00269     * This operator function returns a DMatrixSlice that represents a
00270     * rectangular region of this DMatrixSlice.  This submatrix region
00271     * represents the rows i1 to 12 and colunms j1 to j2.
00272     * 
00273     * Preconditions: <ul>
00274     *   <li> #i1 <= i2# (throw out_of_range)
00275     *   <li> #i2 <= this->rows()# (throw out_of_range)
00276     *   <li> #j1 <= j2# (throw out_of_range)
00277     *   <li> #j2 <= this->cols()# (throw out_of_range)
00278     *   </ul>
00279     */
00280   DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
00281     , size_type j2);
00283   const DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
00284     , size_type j2) const;
00286   DMatrixSlice* operator&() {
00287     return this;
00288   }
00290   const DMatrixSlice* operator&() const {
00291     return this;
00292   }
00294   DMatrixSlice& operator()();
00296   const DMatrixSlice& operator()() const;
00297 
00298   // @}
00299 
00300   /* * @name {\bf Assignment operators}. */
00301   // @{
00302 
00304   /* * Sets all elements = alpha
00305     *
00306     * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to 1 x 1
00307     * and the single element is set to alpha.
00308     *
00309     * Postcondtions: <ul>
00310     *   <li> #this->operator()(i,j) == alpha#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
00311     *   </ul>
00312     */
00313   DMatrixSlice& operator=(value_type alpha);
00315   /* *  Copies all of the elements of the DMatrixSlice, #rhs#, into the elements of #this#.
00316     *
00317     * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to
00318     * the size of the rhs matrix.
00319     *
00320     * Precondtions: <ul>
00321     *   <li> #this->rows() == gms_rhs.rows()# (throw length_error)
00322     *   <li> #this->cols() == gms_rhs.cols()# (throw length_error)
00323     *   </ul>
00324     *
00325     * Postcondtions: <ul>
00326     *   <li> #this->operator()(i,j) == gms_rhs(i,j)#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
00327     *   </ul>
00328     */
00329   DMatrixSlice& operator=(const DMatrixSlice& gms_rhs);
00330 
00331   // @}
00332 
00333   /* * @name {\bf Raw data access}.
00334     */
00335   // @{
00336 
00338   size_type     max_rows() const;
00340   /* * Return pointer to the first element in the underlying array the jth
00341      * col (j is 1-based here [1,cols]).  If unsized col_ptr(1) returns zero if unsized.
00342      */
00343   value_type*     col_ptr(size_type j);
00345   const value_type* col_ptr(size_type j) const;
00346 
00347   // @}
00348 
00349 private:
00350   value_type    *ptr_;    // contains the data for the matrix region
00351   size_type   max_rows_,  // the number of rows in the full matrix
00352           rows_,    // the number of rows in this matrix region
00353           cols_;    // the number of cols in this matrix region
00354 
00355   // Assert the row subscript is in bounds (1-based), (throw std::out_of_range)
00356   void validate_row_subscript(size_type i) const;
00357   // Assert the column subscript is in bounds (1-based), (throw std::out_of_range)
00358   void validate_col_subscript(size_type j) const;
00359   // Assert that a constructed DMatrixSlice has a valid range, (throw std::out_of_range)
00360   void validate_setup(size_type size) const;
00361   
00362   // Get a diagonal
00363   DVectorSlice p_diag(difference_type k) const;
00364 
00365 };  // end class DMatrixSlice
00366 
00368 /* * 2-D General Rectangular Matrix (column major) Storage Class.
00369   *
00370   * This class provides the storage for 2-D rectangular matrices.
00371   */
00372 class DMatrix {
00373 public:
00374   /* * @name {\bf Nested Member Types (STL)}.
00375     *
00376     * These nested types give the types used in the interface to the class.
00377     *
00378     * \begin{description}
00379     * <li>[#value_type#]        type being stored in the underlying valarray<>      
00380     * <li>[#size_type#]       type for the number rows and coluns
00381     * <li>[#difference_type#]   type for the stride between elements (in a row or column)
00382     * <li>[#reference#]       value_type&
00383     * <li>[#const_reference#]   const value_type&
00384     * \end{description}
00385     */
00386   // @{
00387   // @}
00388   
00389   typedef DenseLinAlgPack::value_type         value_type;
00390   typedef DenseLinAlgPack::size_type          size_type;
00391   typedef ptrdiff_t               difference_type;
00392   typedef value_type&               reference;
00393   typedef const value_type&           const_reference;
00394   typedef std::valarray<value_type>       valarray;
00395 
00396   /* * @name {\bf Constructors}.
00397     *
00398     * The general user uses these constructors to create a matrix.  
00399     *
00400     * The default constructor is used and is therefore not shown here.
00401     */
00402   // @{
00403 
00405   DMatrix();
00407   explicit DMatrix(size_type rows, size_type cols);
00409   /* * Construct rectangular matrix (rows x cols) with elements initialized to val.
00410     *
00411     * Postconditions: <ul>
00412     *   <li> #this->operator()(i,j) == val#, i = 1,2,...,#rows#, j = 1,2,...,#cols#  
00413     *   </ul>
00414     */
00415   explicit DMatrix(value_type val, size_type rows, size_type cols);
00417   /* * Construct rectangular matrix (rows x cols) initialized to elements of p (by column).
00418     *
00419     * Postconditions: <ul>
00420     *   <li> #this->operator()(i,j) == p[i-1 + rows * (j - 1)]#, i = 1,2,...,#rows#, j = 1,2,...,#cols#  
00421     *   </ul>
00422     */
00423   explicit DMatrix(const value_type* p, size_type rows, size_type cols);
00425   /* * Construct a matrix from the elements in another DMatrixSlice, #gms#.
00426     *
00427     * Postconditions: <ul>
00428     *   <li> #this->operator()(i,j) == gms(i,j)#, i = 1,2,...,#rows#, j = 1,2,...,#cols#  
00429     *   </ul>
00430     */
00431   DMatrix(const DMatrixSlice& gms);
00432 
00433   // @}
00434 
00435   /* * @name {\bf Memory Management, Dimensionality, Misc}. */
00436   // @{
00437   
00439   void resize(size_type rows, size_type cols, value_type val = value_type());
00440 
00442   void free();
00443   
00445   size_type   rows() const;
00446 
00448   size_type   cols() const;
00449 
00451   /* * Returns the degree of memory overlap of #this# and the DMatrixSlice object #gms#.
00452     *
00453     * @return 
00454     *   \begin{description}
00455     *   <li>[NO_OVERLAP]  There is no memory overlap between #this# and #gms#.
00456     *   <li>[SOME_OVERLAP]  There is some memory locations that #this# and #gms# share.
00457     *   <li>[SAME_MEM]    The DMatrixSlice objects #this# and #gms# share the exact same memory locations.
00458     *   \end{description}
00459     */
00460   EOverLap overlap(const DMatrixSlice& gms) const;
00461 
00462   // @}
00463 
00464   /* * @name {\bf Individual Element Access Subscripting (lvalue)}. */
00465   // @{
00466 
00468   reference   operator()(size_type i, size_type j);
00469 
00471   const_reference operator()(size_type i, size_type j) const;
00472 
00473   // @}
00474 
00475   /* * @name {\bf Subregion Access (1-based)}.
00476     *
00477     * These member functions allow access to subregions of the DMatrix object.
00478     * The functions, row(i), col(j), and diag(k) return DVectorSlice objects while
00479     * the subscripting operators opeator()(I,J) return DMatrixSlice objects for
00480     * rectangular subregions.
00481     */
00482   // @{
00483 
00485   DVectorSlice      row(size_type i);
00486 
00488   const DVectorSlice  row(size_type i) const;
00489 
00491   DVectorSlice      col(size_type j);
00492 
00494   const DVectorSlice  col(size_type j) const;
00495 
00497   /* * Return DVectorSlice object representing a diagonal.
00498     *
00499     * Passing k == 0 returns the center diagonal.  Values of k < 0 are the lower diagonals
00500     * (k = -1, -2, ..., #this->rows()# - 1).  Values of k > 0 are the upper diagonals
00501     * (k = 1, 2, ..., #this->cols()# - 1).
00502     *
00503     * Preconditions: <ul>
00504     *   <li> #[k < 0] k <= this->rows() + 1# (throw out_of_range)
00505     *   <li> #[k > 0] k <= this->cols() + 1# (throw out_of_range)
00506     *   </ul>
00507     */
00508   DVectorSlice      diag(difference_type k = 0);
00509 
00511   const DVectorSlice  diag(difference_type k = 0) const;
00512 
00514   /* * Extract a rectangular subregion containing rows I, and columns J.
00515     *
00516     * This operator function returns a DMatrixSlice that represents a
00517     * rectangular region of this DMatrixSlice.  This submatrix region
00518     * represents the rows I.lbound() to I.ubound() and columns J.lbound()
00519     * to J.lbound().  If I or J is unbounded (full_range() == true, constructed
00520     * with Range1D()), the all of the rows or columns respectively will be
00521     * selected.  For example. To select all the rows and the first 5 columns of
00522     * a matrix #A# you would use #A(Range1D(),Range1D(1,5))#.
00523     *
00524     * Preconditions: <ul>
00525     *   <li> #[I.full_range() == false] I.ubound() <= this->rows()# (throw out_of_range)
00526     *   <li> #[J.full_range() == false] J.ubound() <= this->cols()# (throw out_of_range)
00527     *   </ul>
00528     */
00529   DMatrixSlice operator()(const Range1D& I, const Range1D& J);
00530 
00532   const DMatrixSlice operator()(const Range1D& I, const Range1D& J) const;
00533 
00535   /* * Extract a rectangular subregion containing rows i1 to i2, and columns j1 to j2.
00536     *
00537     * This operator function returns a DMatrixSlice that represents a
00538     * rectangular region of this DMatrixSlice.  This submatrix region
00539     * represents the rows i1 to 12 and colunms j1 to j2.
00540     * 
00541     * Preconditions: <ul>
00542     *   <li> #i1 <= i2# (throw out_of_range)
00543     *   <li> #i2 <= this->rows()# (throw out_of_range)
00544     *   <li> #j1 <= j2# (throw out_of_range)
00545     *   <li> #j2 <= this->cols()# (throw out_of_range)
00546     *   </ul>
00547     */
00548   DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
00549     , size_type j2);
00550 
00552   const DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
00553     , size_type j2) const;
00554 
00556   DMatrixSlice operator()();
00557 
00559   const DMatrixSlice operator()() const;
00560 
00561   // @}
00562 
00563   /* * @name {\bf Implicit conversion operators}.
00564     *
00565     * These functions allow for the implicit converstion from a DMatrix to a DMatrixSlice.
00566     * This implicit converstion is important for the proper usage of much of the
00567     * libraries functionality.
00568     */
00569   // @{
00570 
00572   operator DMatrixSlice();
00574   operator const DMatrixSlice() const;
00575 
00576   // @}
00577 
00578   /* * @name {\bf Assignment Operators}. */
00579   // @{
00580   
00582   /* * Sets all elements = alpha
00583     *
00584     * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to 1 x 1
00585     * and the single element is set to alpha.
00586     *
00587     * Postcondtions: <ul>
00588     *   <li> #this->operator()(i,j) == alpha#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
00589     */
00590   DMatrix& operator=(value_type rhs);
00592   /* * Copies all of the elements of the DMatrixSlice, #rhs#, into the elements of #this#.
00593     *
00594     * If #this# is not the same size as gms_rhs the #this# is resized.
00595     *
00596     * Postcondtions: <ul>
00597     *   <li> #this->operator()(i,j) == gms_rhs(i,j)#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
00598     */
00599   DMatrix& operator=(const DMatrixSlice& gms_rhs);
00601   DMatrix& operator=(const DMatrix& rhs);
00602 
00603   // @}
00604 
00605   /* * @name {\bf Raw data access}.
00606     */
00607   // @{
00608 
00610   size_type     max_rows() const;
00612   /* * Return pointer to the first element in the underlying array the jth
00613      * col (j is 1-based here [1,cols]).  If unsized col_ptr(1) returns zero if unsized.
00614      */
00615   value_type*     col_ptr(size_type j);
00617   const value_type* col_ptr(size_type j) const;
00618 
00619   // @}
00620 
00621 private:
00622   std::valarray<value_type> v_;
00623   size_type         rows_;
00624 
00625   // Assert the row subscript is in bounds (1-based), (throw std::out_of_range)
00626   void validate_row_subscript(size_type i) const;
00627   // Assert the column subscript is in bounds (1-based), (throw std::out_of_range)
00628   void validate_col_subscript(size_type j) const;
00629 
00630   // Get a diagonal, (throw std::out_of_range)
00631   DVectorSlice p_diag(difference_type k) const;
00632 
00633 };  // end class DMatrix
00634 
00635 //    end General Matix Classes scope
00636 // @}
00637 
00638 // ///////////////////////////////////////////////////////////////////////////////
00639 // Non-member function declarations                       //
00640 // ///////////////////////////////////////////////////////////////////////////////
00641 
00642 /* * @name {\bf DMatrix / DMatrixSlice Associated Non-Member Functions}. */
00643 // @{
00644 //    begin non-member functions scope
00645 
00646 inline 
00648 /* * Explicit conversion function from DMatrix to DMatrixSlice.
00649   *
00650   * This is needed to allow a defered evaluation class (TCOL) to be evaluated using its
00651   * implicit conversion operator temp_type() (which returns DMatrix for DMatrixSlice
00652   * resulting expressions).
00653   */
00654 //DMatrixSlice EvaluateToDMatrixSlice(const DMatrix& gm)
00655 //{ return DMatrixSlice(gm); }
00656 
00658 void assert_gms_sizes(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1, const DMatrixSlice& gms2
00659   , BLAS_Cpp::Transp trans2);
00660 
00661 inline 
00663 void assert_gms_square(const DMatrixSlice& gms) {
00664 #ifdef LINALGPACK_CHECK_SLICE_SETUP
00665   if(gms.rows() != gms.cols())
00666     throw std::length_error("Matrix must be square");
00667 #endif
00668 } 
00669 
00670 inline 
00672 /* * Utility to check if a lhs matrix slice is the same size as a rhs matrix slice.
00673   *
00674   * A DMatrixSlice can not be resized since the rows_ property of the
00675   * DMatrix it came from will not be updated.  Allowing a DMatrixSlice
00676   * to resize from unsized would require that the DMatrixSlice carry
00677   * a reference to the DMatrix it was created from.  If this is needed
00678   * then it will be added.
00679   */
00680 void assert_gms_lhs(const DMatrixSlice& gms_lhs, size_type rows, size_type cols
00681   , BLAS_Cpp::Transp trans_rhs = BLAS_Cpp::no_trans)
00682 {
00683   if(trans_rhs == BLAS_Cpp::trans) std::swap(rows,cols);
00684   if(gms_lhs.rows() == rows && gms_lhs.cols() == cols) return; // same size
00685   // not the same size so is an error
00686   throw std::length_error("assert_gms_lhs(...):  lhs DMatrixSlice dim does not match rhs dim");
00687 }
00688 
00689 /* * @name Return rows or columns from a possiblly transposed DMatrix or DMatrixSlice. */
00690 // @{
00691 
00692 inline 
00694 DVectorSlice row(DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type i) {
00695   return (trans ==  BLAS_Cpp::no_trans) ? gms.row(i) : gms.col(i);
00696 } 
00697 
00698 inline 
00700 DVectorSlice col(DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type j) {
00701   return (trans ==  BLAS_Cpp::no_trans) ? gms.col(j) : gms.row(j);
00702 } 
00703 
00704 inline 
00706 const DVectorSlice row(const DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type i) {
00707   return (trans ==  BLAS_Cpp::no_trans) ? gms.row(i) : gms.col(i);
00708 } 
00709 
00710 inline 
00712 const DVectorSlice col(const DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type j) {
00713   return (trans ==  BLAS_Cpp::no_trans) ? gms.col(j) : gms.row(j);
00714 } 
00715 
00716 inline 
00718 DVectorSlice row(DMatrix& gm, BLAS_Cpp::Transp trans, size_type i) {
00719   return (trans ==  BLAS_Cpp::no_trans) ? gm.row(i) : gm.col(i);
00720 } 
00721 
00722 inline 
00724 DVectorSlice col(DMatrix& gm, BLAS_Cpp::Transp trans, size_type j) {
00725   return (trans ==  BLAS_Cpp::no_trans) ? gm.col(j) : gm.row(j);
00726 } 
00727 
00728 inline 
00730 const DVectorSlice row(const DMatrix& gm, BLAS_Cpp::Transp trans, size_type i) {
00731   return (trans ==  BLAS_Cpp::no_trans) ? gm.row(i) : gm.col(i);
00732 } 
00733 
00734 inline 
00736 const DVectorSlice col(const DMatrix& gm, BLAS_Cpp::Transp trans, size_type j) {
00737   return (trans ==  BLAS_Cpp::no_trans) ? gm.col(j) : gm.row(j);
00738 } 
00739 
00740 // @}
00741 
00742 inline 
00744 void resize_gm_lhs(DMatrix* gm_rhs, size_type rows, size_type cols
00745   , BLAS_Cpp::Transp trans_rhs)
00746 {
00747   if(trans_rhs == BLAS_Cpp::trans) std::swap(rows,cols);
00748   gm_rhs->resize(rows,cols);
00749 }
00750 
00751 //    end non-member functions scope
00752 // @}
00753 
00754 //    end General Rectangular 2-D Matrices scope
00755 // @}
00756 
00757 // ////////////////////////////////////////////////////////////////////////////////
00758 // Inline definitions of computationally independent member function definitions //
00759 // ////////////////////////////////////////////////////////////////////////////////
00760 
00761 // /////////////////////////////////////////////////////////////////////////////
00762 // DMatrixSlice inline member function definitions
00763 
00764 // Private utilities
00765 
00766 #ifndef LINALGPACK_CHECK_RANGE
00767 inline
00768 void DMatrixSlice::validate_row_subscript(size_type i) const
00769 {}
00770 #endif
00771 
00772 #ifndef LINALGPACK_CHECK_RANGE
00773 inline
00774 void DMatrixSlice::validate_col_subscript(size_type j) const
00775 {}
00776 #endif
00777 
00778 #ifndef LINALGPACK_CHECK_SLICE_SETUP
00779 inline
00780 void DMatrixSlice::validate_setup(size_type size) const
00781 {}
00782 #endif
00783 
00784 // Constructors
00785 
00786 inline
00787 DMatrixSlice::DMatrixSlice()
00788   : ptr_(0), max_rows_(0), rows_(0), cols_(0)
00789 {}
00790 
00791 inline
00792 DMatrixSlice::DMatrixSlice( value_type* ptr, size_type size
00793     , size_type max_rows, size_type rows, size_type cols )
00794   : ptr_(ptr), max_rows_(max_rows), rows_(rows), cols_(cols)
00795 { 
00796   validate_setup(size);
00797 }
00798 
00799 inline
00800 DMatrixSlice::DMatrixSlice( DMatrixSlice& gms, const Range1D& I
00801     , const Range1D& J)
00802   : ptr_( gms.col_ptr(1) + (I.lbound() - 1) + (J.lbound() - 1) * gms.max_rows() )
00803   , max_rows_(gms.max_rows())
00804   , rows_(I.size())
00805   , cols_(J.size())
00806 { 
00807   gms.validate_row_subscript(I.ubound());
00808   gms.validate_col_subscript(J.ubound());
00809 }
00810 
00811 inline
00812 void DMatrixSlice::bind(DMatrixSlice gms) {
00813   ptr_    = gms.ptr_;
00814   max_rows_ = gms.max_rows_;
00815   rows_   = gms.rows_;
00816   cols_   = gms.cols_;
00817 }
00818 
00819 // Size / Dimensionality
00820 
00821 inline
00822 DMatrixSlice::size_type DMatrixSlice::rows() const {
00823   return rows_;
00824 }
00825 
00826 inline
00827 DMatrixSlice::size_type DMatrixSlice::cols() const {
00828   return cols_;
00829 }
00830 
00831 // Misc
00832 
00833 // Element access
00834 
00835 inline
00836 DMatrixSlice::reference DMatrixSlice::operator()(size_type i, size_type j)
00837 { 
00838   validate_row_subscript(i);
00839   validate_col_subscript(j);
00840   return ptr_[(i-1) + (j-1) * max_rows_];
00841 }
00842 
00843 inline
00844 DMatrixSlice::const_reference DMatrixSlice::operator()(size_type i, size_type j) const
00845 {
00846   validate_row_subscript(i);
00847   validate_col_subscript(j);
00848   return ptr_[(i-1) + (j-1) * max_rows_];
00849 }
00850 
00851 // Subregion access (validated by constructor for DMatrixSlice)
00852 
00853 inline
00854 DVectorSlice  DMatrixSlice::row(size_type i) {
00855   validate_row_subscript(i);
00856   return DVectorSlice( ptr_ + (i-1), cols(), max_rows() );
00857 } 
00858 
00859 inline
00860 const DVectorSlice DMatrixSlice::row(size_type i) const {
00861   validate_row_subscript(i);
00862   return DVectorSlice( const_cast<value_type*>(ptr_) + (i-1), cols(), max_rows() );
00863 } 
00864 
00865 inline
00866 DVectorSlice  DMatrixSlice::col(size_type j) {
00867   validate_col_subscript(j);
00868   return DVectorSlice( ptr_ + (j-1)*max_rows(), rows(), 1 );
00869 } 
00870 
00871 inline
00872 const DVectorSlice DMatrixSlice::col(size_type j) const {
00873   validate_col_subscript(j);
00874   return DVectorSlice( const_cast<value_type*>(ptr_) + (j-1)*max_rows(), rows(), 1 );
00875 } 
00876 
00877 inline
00878 DVectorSlice DMatrixSlice::diag(difference_type k) {
00879   return p_diag(k);
00880 }
00881 
00882 inline
00883 const DVectorSlice DMatrixSlice::diag(difference_type k) const {
00884   return p_diag(k);
00885 }
00886 
00887 inline
00888 DMatrixSlice DMatrixSlice::operator()(const Range1D& I, const Range1D& J) {
00889   return DMatrixSlice(*this, RangePack::full_range(I, 1, rows()), RangePack::full_range(J,1,cols()));
00890 }
00891 
00892 inline
00893 const DMatrixSlice DMatrixSlice::operator()(const Range1D& I, const Range1D& J) const {
00894   return DMatrixSlice( const_cast<DMatrixSlice&>(*this)
00895     , RangePack::full_range(I, 1, rows()), RangePack::full_range(J,1,cols()) );
00896 }
00897 
00898 inline
00899 DMatrixSlice DMatrixSlice::operator()(size_type i1, size_type i2, size_type j1
00900   , size_type j2)
00901 {
00902   return DMatrixSlice(*this, Range1D(i1,i2), Range1D(j1,j2));
00903 }
00904 
00905 inline
00906 const DMatrixSlice DMatrixSlice::operator()(size_type i1, size_type i2, size_type j1
00907   , size_type j2) const
00908 {
00909   return DMatrixSlice( const_cast<DMatrixSlice&>(*this), Range1D(i1,i2)
00910     , Range1D(j1,j2) );
00911 }
00912 
00913 inline
00914 DMatrixSlice& DMatrixSlice::operator()() {
00915   return *this;
00916 }
00917 
00918 inline
00919 const DMatrixSlice& DMatrixSlice::operator()() const {
00920   return *this;
00921 }
00922 
00923 // Assignment operators
00924 
00925 inline
00926 DMatrixSlice& DMatrixSlice::operator=(value_type alpha) {
00927   assign(this, alpha);
00928   return *this;
00929 }
00930 
00931 inline
00932 DMatrixSlice& DMatrixSlice::operator=(const DMatrixSlice& rhs) {
00933   assign(this, rhs, BLAS_Cpp::no_trans);
00934   return *this;
00935 }
00936 
00937 // Raw data access
00938 
00939 inline
00940 DMatrixSlice::size_type DMatrixSlice::max_rows() const
00941 { return max_rows_; }
00942 
00943 inline
00944 DMatrixSlice::value_type* DMatrixSlice::col_ptr(size_type j) {
00945   if( ptr_ )
00946     validate_col_subscript(j);
00947   return ptr_ + (j-1) * max_rows(); // will be 0 if not bound to a view.
00948 }
00949 
00950 inline
00951 const DMatrixSlice::value_type* DMatrixSlice::col_ptr(size_type j) const {
00952   if( ptr_ )
00953     validate_col_subscript(j);
00954   return ptr_ + (j-1) * max_rows(); // will be 0 if not bound to a view.
00955 }
00956 
00957 // ////////////////////////////////////////////////////////////////////////////////////////
00958 // DMatrix inline member function definitions
00959 
00960 // Private utilities
00961 
00962 #ifndef LINALGPACK_CHECK_RANGE
00963 inline
00964 void DMatrix::validate_row_subscript(size_type i) const
00965 {}
00966 #endif
00967 
00968 #ifndef LINALGPACK_CHECK_RANGE
00969 inline
00970 void DMatrix::validate_col_subscript(size_type j) const
00971 {}
00972 #endif
00973 
00974 // constructors
00975 
00976 inline
00977 DMatrix::DMatrix() : v_(), rows_(0)
00978 {}
00979 
00980 inline
00981 DMatrix::DMatrix(size_type rows, size_type cols)
00982   : v_(rows*cols), rows_(rows)
00983 {}
00984 
00985 inline
00986 DMatrix::DMatrix(value_type val, size_type rows, size_type cols)
00987   : v_(val,rows*cols), rows_(rows)
00988 {}
00989 
00990 inline
00991 DMatrix::DMatrix(const value_type* p, size_type rows, size_type cols)
00992   : v_(rows*cols), rows_(rows)
00993 {
00994 // 6/7/00: valarray<> in libstdc++-2.90.7 has a bug in v_(p,size) so we do not
00995 // use it.  This is a hack until I can find the time to remove valarray all
00996 // together.
00997   std::copy( p, p + rows*cols, &v_[0] );
00998 }
00999 
01000 inline
01001 DMatrix::DMatrix(const DMatrixSlice& gms)
01002   : v_(gms.rows() * gms.cols()), rows_(gms.rows())
01003 { 
01004   assign(this, gms, BLAS_Cpp::no_trans);
01005 }
01006 
01007 // Memory management
01008 
01009 inline
01010 void DMatrix::resize(size_type rows, size_type cols, value_type val)
01011 {
01012   v_.resize(rows*cols,val);
01013   v_ = val;
01014   rows_ = rows;
01015 }
01016 
01017 inline
01018 void DMatrix::free() {
01019   v_.resize(0);
01020   rows_ = 0;
01021 }
01022 
01023 // Size / Dimensionality
01024 
01025 inline
01026 DMatrix::size_type DMatrix::rows() const {
01027   return rows_;
01028 }
01029 
01030 inline
01031 DMatrix::size_type DMatrix::cols() const {
01032   return rows_ > 0 ? v_.size() / rows_ : 0;
01033 }
01034 
01035 // Element access
01036 
01037 inline
01038 DMatrix::reference DMatrix::operator()(size_type i, size_type j)
01039 {  
01040   validate_row_subscript(i); validate_col_subscript(j);
01041   return v_[(i-1) + (j-1) * rows_];
01042 }
01043 
01044 inline
01045 DMatrix::const_reference DMatrix::operator()(size_type i, size_type j) const
01046 {
01047   validate_row_subscript(i); validate_col_subscript(j);
01048   return (const_cast<std::valarray<value_type>&>(v_))[(i-1) + (j-1) * rows_];
01049 }
01050 
01051 // subregion access (range checked by constructors)
01052 
01053 inline
01054 DVectorSlice DMatrix::row(size_type i)
01055 {
01056   validate_row_subscript(i);
01057   return DVectorSlice( col_ptr(1) + (i-1), cols(), rows() );
01058 } 
01059 
01060 inline
01061 const DVectorSlice DMatrix::row(size_type i) const
01062 {
01063   validate_row_subscript(i);
01064   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + (i-1), cols(), rows() );
01065 } 
01066 
01067 inline
01068 DVectorSlice  DMatrix::col(size_type j)
01069 {
01070   validate_col_subscript(j);
01071   return DVectorSlice( col_ptr(1) + (j-1) * rows(), rows(), 1 );
01072 } 
01073 
01074 inline
01075 const DVectorSlice DMatrix::col(size_type j) const
01076 {
01077   validate_col_subscript(j);
01078   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + (j-1) * rows(), rows(), 1 ) ;
01079 } 
01080 
01081 inline
01082 DVectorSlice DMatrix::diag(difference_type k)
01083 {
01084   return p_diag(k);
01085 } 
01086 
01087 inline
01088 const DVectorSlice DMatrix::diag(difference_type k) const
01089 {
01090   return p_diag(k);
01091 } 
01092 
01093 inline
01094 DMatrixSlice DMatrix::operator()(const Range1D& I, const Range1D& J)
01095 {
01096   Range1D Ix = RangePack::full_range(I,1,rows()), Jx = RangePack::full_range(J,1,cols());
01097   return DMatrixSlice( col_ptr(1) + (Ix.lbound() - 1) + (Jx.lbound() - 1) * rows()
01098     , max_rows() * cols(), max_rows(), Ix.size(), Jx.size() );
01099 }
01100 
01101 inline
01102 const DMatrixSlice DMatrix::operator()(const Range1D& I, const Range1D& J) const
01103 {
01104   Range1D Ix = RangePack::full_range(I,1,rows()), Jx = RangePack::full_range(J,1,cols());
01105   return DMatrixSlice( const_cast<value_type*>(col_ptr(1)) + (Ix.lbound() - 1) + (Jx.lbound() - 1) * rows()
01106     , max_rows() * cols(), max_rows(), Ix.size(), Jx.size() );
01107 }
01108 
01109 inline
01110 DMatrixSlice DMatrix::operator()(size_type i1, size_type i2, size_type j1
01111   , size_type j2)
01112 {
01113   return DMatrixSlice( col_ptr(1) + (i1 - 1) + (j1 - 1) * rows()
01114     , max_rows() * cols(), max_rows(), i2 - i1 + 1, j2 - j1 + 1 );
01115 }
01116 
01117 inline
01118 const DMatrixSlice DMatrix::operator()(size_type i1, size_type i2, size_type j1
01119   , size_type j2) const
01120 {
01121   return DMatrixSlice( const_cast<value_type*>(col_ptr(1)) + (i1 - 1) + (j1 - 1) * rows()
01122     , max_rows() * cols(), max_rows(), i2 - i1 + 1, j2 - j1 + 1 );
01123 }
01124 
01125 inline
01126 DMatrixSlice DMatrix::operator()()
01127 {
01128   return DMatrixSlice( col_ptr(1), max_rows() * cols(), max_rows(), rows(), cols() );
01129 }
01130 
01131 inline
01132 const DMatrixSlice DMatrix::operator()() const
01133 {
01134   return DMatrixSlice( const_cast<value_type*>(col_ptr(1)), max_rows() * cols(), max_rows()
01135     , rows(), cols() );
01136 }
01137 
01138 // Implicit conversion operators
01139 
01140 inline
01141 DMatrix::operator DMatrixSlice() {
01142   return (*this)();
01143 }
01144 
01145 inline
01146 DMatrix::operator const DMatrixSlice() const
01147 {
01148   return (*this)();
01149 }
01150 
01151 // Assignment operators
01152 
01153 inline
01154 DMatrix& DMatrix::operator=(value_type alpha)
01155 {
01156   assign(this, alpha);
01157   return *this;
01158 }
01159 
01160 inline
01161 DMatrix& DMatrix::operator=(const DMatrix& rhs)
01162 {
01163   assign(this, rhs, BLAS_Cpp::no_trans);
01164   return *this;
01165 }
01166 
01167 inline
01168 DMatrix& DMatrix::operator=(const DMatrixSlice& rhs)
01169 {
01170   assign(this, rhs, BLAS_Cpp::no_trans);
01171   return *this;
01172 }
01173 
01174 // Raw data access
01175 
01176 inline
01177 DMatrix::size_type DMatrix::max_rows() const
01178 { return rows_; }
01179 
01180 inline
01181 DMatrix::value_type* DMatrix::col_ptr(size_type j)
01182 {
01183   if( v_.size() ) {
01184     validate_col_subscript(j);
01185     return &v_[ (j-1) * max_rows() ];
01186   }
01187   else {
01188     return 0;
01189   }
01190 }
01191 
01192 inline
01193 const DMatrix::value_type* DMatrix::col_ptr(size_type j) const 
01194 {
01195   if( v_.size() ) {
01196     validate_col_subscript(j);
01197     return &const_cast<valarray&>(v_)[ (j-1) * max_rows() ];
01198   }
01199   else {
01200     return 0;
01201   }
01202 }
01203 
01204 } // end namespace DenseLinAlgPack
01205 
01206 #endif  // GEN_MATRIX_CLASS_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines