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

Generated on Tue Oct 20 12:51:46 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7