DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DenseLinAlgPack_DMatrixClass.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 //
42 // General matrix and matrix region (slice) classes
43 
44 #ifndef GEN_MATRIX_CLASS_H
45 #define GEN_MATRIX_CLASS_H
46 
47 #include "DenseLinAlgPack_DVectorClass.hpp"
48 #include "DenseLinAlgPack_DMatrixAssign.hpp"
49 
50 /* * @name {\bf Dense 2-D Rectangular Matrix Absractions}.
51  *
52  * The class DMatrix is a storage class for 2-D matrices while the class DMatrixSlice
53  * is used to represent rectangular regions of a DMatrix object.
54  */
55 // @{
56 // begin General Rectangular 2-D Matrices scope
57 
58 namespace DenseLinAlgPack {
59 
60 class DMatrix;
61 
62 /* * @name {\bf General Matrix Classes}. */
63 // @{
64 
65 // ////////////////////////////////////////////////////////////////////////////////////////////////////////
66 // GenMatrixClass
67 //
68 
70 /* * 2-D General Rectangular Matrix Subregion Class (Slice) (column major).
71  *
72  * This class is used to represent a rectangular matrix. It uses a BLAS-like
73  * slice of a raw C++ array. Objects of this class can represent
74  * an entire matrix or any rectangular subregion of a matrix.
75  */
76 class DMatrixSlice {
77 public:
78  /* * @name {\bf Nested Member Types (STL)}.
79  *
80  * These nested types give the types used in the interface to the class.
81  *
82  * \begin{description}
83  * <li>[#value_type#] - type being stored in the underlying valarray<>
84  * <li>[#size_type#] - type for the number rows and coluns
85  * <li>[#difference_type#] - type for the stride between elements (in a row or column)
86  * <li>[#reference#] - value_type&
87  * <li>[#const_reference#] - const value_type&
88  * \end{description}
89  */
90  // @{
91  // @}
92  typedef DenseLinAlgPack::value_type value_type;
93  typedef DenseLinAlgPack::size_type size_type;
94  typedef ptrdiff_t difference_type;
95  typedef value_type& reference;
96  typedef const value_type& const_reference;
97 
98  /* * @name {\bf Constructors}.
99  *
100  * These constructors are used by the other entities in the library
101  * to create DMatrixSlices. In general user need not use these
102  * constructors directly. Instead, the general user should use the
103  * subscript operators to create subregions of DMatrix and DMatrixSlice
104  * objects.
105  *
106  * The default copy constructor is used and is therefore not shown here.
107  */
108 
109  // @{
110 
112  /* * Construct an empty view.
113  *
114  * The client can then call bind(...) to bind the view.
115  */
116  DMatrixSlice();
117 
119  /* * Construct a veiw of a matrix from a raw C++ array.
120  *
121  * The DMatrixSlice constructed represents a 2-D matrix whos elements are stored
122  * in the array starting at #ptr#. This is how the BLAS represent general rectangular
123  * matrices.
124  * The class can be used to provide a non-constant view the elements (#DMatrix#)
125  * or a constant view (#const DMatrixSlice#). Here is an example of how to
126  * create a constant view:
127  *
128  \verbatim
129  const DMatrixSlice::size_type m = 4, n = 4;
130  const DMatrixSlice::value_type ptr[m*n] = { ... };
131  const GenMatrixslice mat(cosnt_cast<DMatrixSlice::value_type*>(ptr),m*n,m,m,n);
132  \endverbatim
133  *
134  * The #const_cast<...># such as in the example above is perfectly safe to use
135  * when constructing #const# veiw of #const# data. On the other hand casting
136  * away #const# and then non-#const# use is not safe in general since the original
137  * #const# data may reside in ROM for example. By using a non-#const# pointer in the
138  * constructor you avoid accidentally making a non-#const# view of #const# data.
139  *
140  * Preconditions: <ul>
141  * <li> #rows <= max_rows# (throw out_of_range)
142  * <li> #start + rows + max_rows * (cols - 1) <= v.size()# (throw out_of_range)
143  * </ul>
144  *
145  * @param ptr pointer to the storage of the elements of the matrix (column oriented).
146  * @param size total size of the storage pointed to by #ptr# (for size checking)
147  * @param max_rows number of rows in the full matrix this sub-matrix was taken from.
148  * This is equivalent to the leading dimension argument (LDA) in
149  * the BLAS.
150  * @param rows number of rows this matrix represents
151  * @param cols number of columns this matrix represents
152  */
153  DMatrixSlice( value_type* ptr, size_type size
154  , size_type max_rows, size_type rows, size_type cols );
155 
157  /* * Construct a submatrix region of another DMatrixSlice.
158  *
159  * This constructor simplifies the creation of subregions using the subscript
160  * operators.
161  *
162  * I and J must be bounded ranges (full_range() == false).
163  *
164  * Preconditions: <ul>
165  * <li> #I.full_range() == false# (throw out_of_range)
166  * <li> #J.full_range() == false# (throw out_of_range)
167  * <li> #I.ubound() <= gms.rows()# (throw out_of_range)
168  * <li> #J.ubound() <= gms.cols()# (throw out_of_range)
169  * </ul>
170  */
171  DMatrixSlice( DMatrixSlice& gms, const Range1D& I
172  , const Range1D& J );
173 
174  // @}
175 
177  /* * Set to the view of the input DMatrixSlice.
178  *
179  *
180  */
181  void bind( DMatrixSlice gms );
182 
183  /* * @name {\bf Dimensionality, Misc}. */
184  // @{
185 
187  size_type rows() const;
189  size_type cols() const;
190 
192  /* * Returns the degree of memory overlap of #this# and the DMatrixSlice object #gms#.
193  *
194  * @return
195  * \begin{description}
196  * <li>[NO_OVERLAP] There is no memory overlap between #this# and #gms#.
197  * <li>[SOME_OVERLAP] There is some memory locations that #this# and #gms# share.
198  * <li>[SAME_MEM] The DMatrixSlice objects #this# and #gms# share the exact same memory locations.
199  * \end{description}
200  */
201  EOverLap overlap(const DMatrixSlice& gms) const;
202 
203  // @}
204 
205  /* * @name {\bf Individual Element Access Subscripting (lvalue)}. */
206  // @{
207 
209  reference operator()(size_type i, size_type j);
211  const_reference operator()(size_type i, size_type j) const;
212 
213  // @}
214 
215  /* * @name {\bf Subregion Access (1-based)}.
216  *
217  * These member functions allow access to subregions of the DMatrixSlice object.
218  * The functions, row(i), col(j), and diag(k) return DVectorSlice objects while
219  * the subscripting operators opeator()(I,J) return DMatrixSlice objects for
220  * rectangular subregions.
221  */
222  // @{
223 
225  DVectorSlice row(size_type i);
227  const DVectorSlice row(size_type i) const;
229  DVectorSlice col(size_type j);
231  const DVectorSlice col(size_type j) const;
233  /* * Return DVectorSlice object representing a diagonal.
234  *
235  * Passing k == 0 returns the center diagonal. Values of k < 0 are the lower diagonals
236  * (k = -1, -2, ..., -#this->rows()# + 1). Values of k > 0 are the upper diagonals
237  * (k = 1, 2, ..., #this->cols()# - 1).
238  *
239  * Preconditions: <ul>
240  * <li> #[k < 0] k <= this->rows() + 1# (throw out_of_range)
241  * <li> #[k > 0] k <= this->cols() + 1# (throw out_of_range)
242  * </ul>
243  */
244  DVectorSlice diag(difference_type k = 0);
246  const DVectorSlice diag(difference_type k = 0) const;
248  /* * Extract a rectangular subregion containing rows I, and columns J.
249  *
250  * This operator function returns a DMatrixSlice that represents a
251  * rectangular region of this DMatrixSlice. This submatrix region
252  * represents the rows I.lbound() to I.ubound() and columns J.lbound()
253  * to J.lbound(). If I or J is unbounded (full_range() == true, constructed
254  * with Range1D()), the all of the rows or columns respectively will be
255  * selected. For example. To select all the rows and the first 5 columns of
256  * a matrix #A# you would use #A(Range1D(),Range1D(1,5))#.
257  *
258  * Preconditions: <ul>
259  * <li> #[I.full_range() == false] I.ubound() <= this->rows()# (throw out_of_range)
260  * <li> #[J.full_range() == false] J.ubound() <= this->cols()# (throw out_of_range)
261  * </ul>
262  */
263  DMatrixSlice operator()(const Range1D& I, const Range1D& J);
265  const DMatrixSlice operator()(const Range1D& I, const Range1D& J) const;
267  /* * Extract a rectangular subregion containing rows i1 to i2, and columns j1 to j2.
268  *
269  * This operator function returns a DMatrixSlice that represents a
270  * rectangular region of this DMatrixSlice. This submatrix region
271  * represents the rows i1 to 12 and colunms j1 to j2.
272  *
273  * Preconditions: <ul>
274  * <li> #i1 <= i2# (throw out_of_range)
275  * <li> #i2 <= this->rows()# (throw out_of_range)
276  * <li> #j1 <= j2# (throw out_of_range)
277  * <li> #j2 <= this->cols()# (throw out_of_range)
278  * </ul>
279  */
280  DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
281  , size_type j2);
283  const DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
284  , size_type j2) const;
287  return this;
288  }
290  const DMatrixSlice* operator&() const {
291  return this;
292  }
296  const DMatrixSlice& operator()() const;
297 
298  // @}
299 
300  /* * @name {\bf Assignment operators}. */
301  // @{
302 
304  /* * Sets all elements = alpha
305  *
306  * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to 1 x 1
307  * and the single element is set to alpha.
308  *
309  * Postcondtions: <ul>
310  * <li> #this->operator()(i,j) == alpha#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
311  * </ul>
312  */
313  DMatrixSlice& operator=(value_type alpha);
315  /* * Copies all of the elements of the DMatrixSlice, #rhs#, into the elements of #this#.
316  *
317  * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to
318  * the size of the rhs matrix.
319  *
320  * Precondtions: <ul>
321  * <li> #this->rows() == gms_rhs.rows()# (throw length_error)
322  * <li> #this->cols() == gms_rhs.cols()# (throw length_error)
323  * </ul>
324  *
325  * Postcondtions: <ul>
326  * <li> #this->operator()(i,j) == gms_rhs(i,j)#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
327  * </ul>
328  */
329  DMatrixSlice& operator=(const DMatrixSlice& gms_rhs);
330 
331  // @}
332 
333  /* * @name {\bf Raw data access}.
334  */
335  // @{
336 
338  size_type max_rows() const;
340  /* * Return pointer to the first element in the underlying array the jth
341  * col (j is 1-based here [1,cols]). If unsized col_ptr(1) returns zero if unsized.
342  */
343  value_type* col_ptr(size_type j);
345  const value_type* col_ptr(size_type j) const;
346 
347  // @}
348 
349 private:
350  value_type *ptr_; // contains the data for the matrix region
351  size_type max_rows_, // the number of rows in the full matrix
352  rows_, // the number of rows in this matrix region
353  cols_; // the number of cols in this matrix region
354 
355  // Assert the row subscript is in bounds (1-based), (throw std::out_of_range)
356  void validate_row_subscript(size_type i) const;
357  // Assert the column subscript is in bounds (1-based), (throw std::out_of_range)
358  void validate_col_subscript(size_type j) const;
359  // Assert that a constructed DMatrixSlice has a valid range, (throw std::out_of_range)
360  void validate_setup(size_type size) const;
361 
362  // Get a diagonal
363  DVectorSlice p_diag(difference_type k) const;
364 
365 }; // end class DMatrixSlice
366 
368 /* * 2-D General Rectangular Matrix (column major) Storage Class.
369  *
370  * This class provides the storage for 2-D rectangular matrices.
371  */
372 class DMatrix {
373 public:
374  /* * @name {\bf Nested Member Types (STL)}.
375  *
376  * These nested types give the types used in the interface to the class.
377  *
378  * \begin{description}
379  * <li>[#value_type#] type being stored in the underlying valarray<>
380  * <li>[#size_type#] type for the number rows and coluns
381  * <li>[#difference_type#] type for the stride between elements (in a row or column)
382  * <li>[#reference#] value_type&
383  * <li>[#const_reference#] const value_type&
384  * \end{description}
385  */
386  // @{
387  // @}
388 
389  typedef DenseLinAlgPack::value_type value_type;
390  typedef DenseLinAlgPack::size_type size_type;
391  typedef ptrdiff_t difference_type;
392  typedef value_type& reference;
393  typedef const value_type& const_reference;
394  typedef std::valarray<value_type> valarray;
395 
396  /* * @name {\bf Constructors}.
397  *
398  * The general user uses these constructors to create a matrix.
399  *
400  * The default constructor is used and is therefore not shown here.
401  */
402  // @{
403 
405  DMatrix();
407  explicit DMatrix(size_type rows, size_type cols);
409  /* * Construct rectangular matrix (rows x cols) with elements initialized to val.
410  *
411  * Postconditions: <ul>
412  * <li> #this->operator()(i,j) == val#, i = 1,2,...,#rows#, j = 1,2,...,#cols#
413  * </ul>
414  */
415  explicit DMatrix(value_type val, size_type rows, size_type cols);
417  /* * Construct rectangular matrix (rows x cols) initialized to elements of p (by column).
418  *
419  * Postconditions: <ul>
420  * <li> #this->operator()(i,j) == p[i-1 + rows * (j - 1)]#, i = 1,2,...,#rows#, j = 1,2,...,#cols#
421  * </ul>
422  */
423  explicit DMatrix(const value_type* p, size_type rows, size_type cols);
425  /* * Construct a matrix from the elements in another DMatrixSlice, #gms#.
426  *
427  * Postconditions: <ul>
428  * <li> #this->operator()(i,j) == gms(i,j)#, i = 1,2,...,#rows#, j = 1,2,...,#cols#
429  * </ul>
430  */
431  DMatrix(const DMatrixSlice& gms);
432 
433  // @}
434 
435  /* * @name {\bf Memory Management, Dimensionality, Misc}. */
436  // @{
437 
439  void resize(size_type rows, size_type cols, value_type val = value_type());
440 
442  void free();
443 
445  size_type rows() const;
446 
448  size_type cols() const;
449 
451  /* * Returns the degree of memory overlap of #this# and the DMatrixSlice object #gms#.
452  *
453  * @return
454  * \begin{description}
455  * <li>[NO_OVERLAP] There is no memory overlap between #this# and #gms#.
456  * <li>[SOME_OVERLAP] There is some memory locations that #this# and #gms# share.
457  * <li>[SAME_MEM] The DMatrixSlice objects #this# and #gms# share the exact same memory locations.
458  * \end{description}
459  */
460  EOverLap overlap(const DMatrixSlice& gms) const;
461 
462  // @}
463 
464  /* * @name {\bf Individual Element Access Subscripting (lvalue)}. */
465  // @{
466 
468  reference operator()(size_type i, size_type j);
469 
471  const_reference operator()(size_type i, size_type j) const;
472 
473  // @}
474 
475  /* * @name {\bf Subregion Access (1-based)}.
476  *
477  * These member functions allow access to subregions of the DMatrix object.
478  * The functions, row(i), col(j), and diag(k) return DVectorSlice objects while
479  * the subscripting operators opeator()(I,J) return DMatrixSlice objects for
480  * rectangular subregions.
481  */
482  // @{
483 
485  DVectorSlice row(size_type i);
486 
488  const DVectorSlice row(size_type i) const;
489 
491  DVectorSlice col(size_type j);
492 
494  const DVectorSlice col(size_type j) const;
495 
497  /* * Return DVectorSlice object representing a diagonal.
498  *
499  * Passing k == 0 returns the center diagonal. Values of k < 0 are the lower diagonals
500  * (k = -1, -2, ..., #this->rows()# - 1). Values of k > 0 are the upper diagonals
501  * (k = 1, 2, ..., #this->cols()# - 1).
502  *
503  * Preconditions: <ul>
504  * <li> #[k < 0] k <= this->rows() + 1# (throw out_of_range)
505  * <li> #[k > 0] k <= this->cols() + 1# (throw out_of_range)
506  * </ul>
507  */
508  DVectorSlice diag(difference_type k = 0);
509 
511  const DVectorSlice diag(difference_type k = 0) const;
512 
514  /* * Extract a rectangular subregion containing rows I, and columns J.
515  *
516  * This operator function returns a DMatrixSlice that represents a
517  * rectangular region of this DMatrixSlice. This submatrix region
518  * represents the rows I.lbound() to I.ubound() and columns J.lbound()
519  * to J.lbound(). If I or J is unbounded (full_range() == true, constructed
520  * with Range1D()), the all of the rows or columns respectively will be
521  * selected. For example. To select all the rows and the first 5 columns of
522  * a matrix #A# you would use #A(Range1D(),Range1D(1,5))#.
523  *
524  * Preconditions: <ul>
525  * <li> #[I.full_range() == false] I.ubound() <= this->rows()# (throw out_of_range)
526  * <li> #[J.full_range() == false] J.ubound() <= this->cols()# (throw out_of_range)
527  * </ul>
528  */
529  DMatrixSlice operator()(const Range1D& I, const Range1D& J);
530 
532  const DMatrixSlice operator()(const Range1D& I, const Range1D& J) const;
533 
535  /* * Extract a rectangular subregion containing rows i1 to i2, and columns j1 to j2.
536  *
537  * This operator function returns a DMatrixSlice that represents a
538  * rectangular region of this DMatrixSlice. This submatrix region
539  * represents the rows i1 to 12 and colunms j1 to j2.
540  *
541  * Preconditions: <ul>
542  * <li> #i1 <= i2# (throw out_of_range)
543  * <li> #i2 <= this->rows()# (throw out_of_range)
544  * <li> #j1 <= j2# (throw out_of_range)
545  * <li> #j2 <= this->cols()# (throw out_of_range)
546  * </ul>
547  */
548  DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
549  , size_type j2);
550 
552  const DMatrixSlice operator()(size_type i1, size_type i2, size_type j1
553  , size_type j2) const;
554 
557 
559  const DMatrixSlice operator()() const;
560 
561  // @}
562 
563  /* * @name {\bf Implicit conversion operators}.
564  *
565  * These functions allow for the implicit converstion from a DMatrix to a DMatrixSlice.
566  * This implicit converstion is important for the proper usage of much of the
567  * libraries functionality.
568  */
569  // @{
570 
572  operator DMatrixSlice();
574  operator const DMatrixSlice() const;
575 
576  // @}
577 
578  /* * @name {\bf Assignment Operators}. */
579  // @{
580 
582  /* * Sets all elements = alpha
583  *
584  * If the underlying valarray is unsized (#this->v().size() == 0#) the matrix is sized to 1 x 1
585  * and the single element is set to alpha.
586  *
587  * Postcondtions: <ul>
588  * <li> #this->operator()(i,j) == alpha#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
589  */
590  DMatrix& operator=(value_type rhs);
592  /* * Copies all of the elements of the DMatrixSlice, #rhs#, into the elements of #this#.
593  *
594  * If #this# is not the same size as gms_rhs the #this# is resized.
595  *
596  * Postcondtions: <ul>
597  * <li> #this->operator()(i,j) == gms_rhs(i,j)#, i = 1,2,...,#this->rows()#, j = 1,2,...,#this->cols()#
598  */
599  DMatrix& operator=(const DMatrixSlice& gms_rhs);
601  DMatrix& operator=(const DMatrix& rhs);
602 
603  // @}
604 
605  /* * @name {\bf Raw data access}.
606  */
607  // @{
608 
610  size_type max_rows() const;
612  /* * Return pointer to the first element in the underlying array the jth
613  * col (j is 1-based here [1,cols]). If unsized col_ptr(1) returns zero if unsized.
614  */
615  value_type* col_ptr(size_type j);
617  const value_type* col_ptr(size_type j) const;
618 
619  // @}
620 
621 private:
622  std::valarray<value_type> v_;
623  size_type rows_;
624 
625  // Assert the row subscript is in bounds (1-based), (throw std::out_of_range)
626  void validate_row_subscript(size_type i) const;
627  // Assert the column subscript is in bounds (1-based), (throw std::out_of_range)
628  void validate_col_subscript(size_type j) const;
629 
630  // Get a diagonal, (throw std::out_of_range)
631  DVectorSlice p_diag(difference_type k) const;
632 
633 }; // end class DMatrix
634 
635 // end General Matix Classes scope
636 // @}
637 
638 // ///////////////////////////////////////////////////////////////////////////////
639 // Non-member function declarations //
640 // ///////////////////////////////////////////////////////////////////////////////
641 
642 /* * @name {\bf DMatrix / DMatrixSlice Associated Non-Member Functions}. */
643 // @{
644 // begin non-member functions scope
645 
646 inline
648 /* * Explicit conversion function from DMatrix to DMatrixSlice.
649  *
650  * This is needed to allow a defered evaluation class (TCOL) to be evaluated using its
651  * implicit conversion operator temp_type() (which returns DMatrix for DMatrixSlice
652  * resulting expressions).
653  */
654 //DMatrixSlice EvaluateToDMatrixSlice(const DMatrix& gm)
655 //{ return DMatrixSlice(gm); }
656 
658 void assert_gms_sizes(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1, const DMatrixSlice& gms2
659  , BLAS_Cpp::Transp trans2);
660 
661 inline
663 void assert_gms_square(const DMatrixSlice& gms) {
664 #ifdef LINALGPACK_CHECK_SLICE_SETUP
665  if(gms.rows() != gms.cols())
666  throw std::length_error("Matrix must be square");
667 #endif
668 }
669 
670 inline
672 /* * Utility to check if a lhs matrix slice is the same size as a rhs matrix slice.
673  *
674  * A DMatrixSlice can not be resized since the rows_ property of the
675  * DMatrix it came from will not be updated. Allowing a DMatrixSlice
676  * to resize from unsized would require that the DMatrixSlice carry
677  * a reference to the DMatrix it was created from. If this is needed
678  * then it will be added.
679  */
680 void assert_gms_lhs(const DMatrixSlice& gms_lhs, size_type rows, size_type cols
681  , BLAS_Cpp::Transp trans_rhs = BLAS_Cpp::no_trans)
682 {
683  if(trans_rhs == BLAS_Cpp::trans) std::swap(rows,cols);
684  if(gms_lhs.rows() == rows && gms_lhs.cols() == cols) return; // same size
685  // not the same size so is an error
686  throw std::length_error("assert_gms_lhs(...): lhs DMatrixSlice dim does not match rhs dim");
687 }
688 
689 /* * @name Return rows or columns from a possiblly transposed DMatrix or DMatrixSlice. */
690 // @{
691 
692 inline
694 DVectorSlice row(DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type i) {
695  return (trans == BLAS_Cpp::no_trans) ? gms.row(i) : gms.col(i);
696 }
697 
698 inline
700 DVectorSlice col(DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type j) {
701  return (trans == BLAS_Cpp::no_trans) ? gms.col(j) : gms.row(j);
702 }
703 
704 inline
706 const DVectorSlice row(const DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type i) {
707  return (trans == BLAS_Cpp::no_trans) ? gms.row(i) : gms.col(i);
708 }
709 
710 inline
712 const DVectorSlice col(const DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type j) {
713  return (trans == BLAS_Cpp::no_trans) ? gms.col(j) : gms.row(j);
714 }
715 
716 inline
718 DVectorSlice row(DMatrix& gm, BLAS_Cpp::Transp trans, size_type i) {
719  return (trans == BLAS_Cpp::no_trans) ? gm.row(i) : gm.col(i);
720 }
721 
722 inline
724 DVectorSlice col(DMatrix& gm, BLAS_Cpp::Transp trans, size_type j) {
725  return (trans == BLAS_Cpp::no_trans) ? gm.col(j) : gm.row(j);
726 }
727 
728 inline
730 const DVectorSlice row(const DMatrix& gm, BLAS_Cpp::Transp trans, size_type i) {
731  return (trans == BLAS_Cpp::no_trans) ? gm.row(i) : gm.col(i);
732 }
733 
734 inline
736 const DVectorSlice col(const DMatrix& gm, BLAS_Cpp::Transp trans, size_type j) {
737  return (trans == BLAS_Cpp::no_trans) ? gm.col(j) : gm.row(j);
738 }
739 
740 // @}
741 
742 inline
744 void resize_gm_lhs(DMatrix* gm_rhs, size_type rows, size_type cols
745  , BLAS_Cpp::Transp trans_rhs)
746 {
747  if(trans_rhs == BLAS_Cpp::trans) std::swap(rows,cols);
748  gm_rhs->resize(rows,cols);
749 }
750 
751 // end non-member functions scope
752 // @}
753 
754 // end General Rectangular 2-D Matrices scope
755 // @}
756 
757 // ////////////////////////////////////////////////////////////////////////////////
758 // Inline definitions of computationally independent member function definitions //
759 // ////////////////////////////////////////////////////////////////////////////////
760 
761 // /////////////////////////////////////////////////////////////////////////////
762 // DMatrixSlice inline member function definitions
763 
764 // Private utilities
765 
766 #ifndef LINALGPACK_CHECK_RANGE
767 inline
768 void DMatrixSlice::validate_row_subscript(size_type i) const
769 {}
770 #endif
771 
772 #ifndef LINALGPACK_CHECK_RANGE
773 inline
774 void DMatrixSlice::validate_col_subscript(size_type j) const
775 {}
776 #endif
777 
778 #ifndef LINALGPACK_CHECK_SLICE_SETUP
779 inline
780 void DMatrixSlice::validate_setup(size_type size) const
781 {}
782 #endif
783 
784 // Constructors
785 
786 inline
787 DMatrixSlice::DMatrixSlice()
788  : ptr_(0), max_rows_(0), rows_(0), cols_(0)
789 {}
790 
791 inline
792 DMatrixSlice::DMatrixSlice( value_type* ptr, size_type size
793  , size_type max_rows, size_type rows, size_type cols )
794  : ptr_(ptr), max_rows_(max_rows), rows_(rows), cols_(cols)
795 {
796  validate_setup(size);
797 }
798 
799 inline
800 DMatrixSlice::DMatrixSlice( DMatrixSlice& gms, const Range1D& I
801  , const Range1D& J)
802  : ptr_( gms.col_ptr(1) + (I.lbound() - 1) + (J.lbound() - 1) * gms.max_rows() )
803  , max_rows_(gms.max_rows())
804  , rows_(I.size())
805  , cols_(J.size())
806 {
807  gms.validate_row_subscript(I.ubound());
808  gms.validate_col_subscript(J.ubound());
809 }
810 
811 inline
813  ptr_ = gms.ptr_;
814  max_rows_ = gms.max_rows_;
815  rows_ = gms.rows_;
816  cols_ = gms.cols_;
817 }
818 
819 // Size / Dimensionality
820 
821 inline
822 DMatrixSlice::size_type DMatrixSlice::rows() const {
823  return rows_;
824 }
825 
826 inline
827 DMatrixSlice::size_type DMatrixSlice::cols() const {
828  return cols_;
829 }
830 
831 // Misc
832 
833 // Element access
834 
835 inline
836 DMatrixSlice::reference DMatrixSlice::operator()(size_type i, size_type j)
837 {
838  validate_row_subscript(i);
839  validate_col_subscript(j);
840  return ptr_[(i-1) + (j-1) * max_rows_];
841 }
842 
843 inline
844 DMatrixSlice::const_reference DMatrixSlice::operator()(size_type i, size_type j) const
845 {
846  validate_row_subscript(i);
847  validate_col_subscript(j);
848  return ptr_[(i-1) + (j-1) * max_rows_];
849 }
850 
851 // Subregion access (validated by constructor for DMatrixSlice)
852 
853 inline
855  validate_row_subscript(i);
856  return DVectorSlice( ptr_ + (i-1), cols(), max_rows() );
857 }
858 
859 inline
860 const DVectorSlice DMatrixSlice::row(size_type i) const {
861  validate_row_subscript(i);
862  return DVectorSlice( const_cast<value_type*>(ptr_) + (i-1), cols(), max_rows() );
863 }
864 
865 inline
867  validate_col_subscript(j);
868  return DVectorSlice( ptr_ + (j-1)*max_rows(), rows(), 1 );
869 }
870 
871 inline
872 const DVectorSlice DMatrixSlice::col(size_type j) const {
873  validate_col_subscript(j);
874  return DVectorSlice( const_cast<value_type*>(ptr_) + (j-1)*max_rows(), rows(), 1 );
875 }
876 
877 inline
878 DVectorSlice DMatrixSlice::diag(difference_type k) {
879  return p_diag(k);
880 }
881 
882 inline
883 const DVectorSlice DMatrixSlice::diag(difference_type k) const {
884  return p_diag(k);
885 }
886 
887 inline
889  return DMatrixSlice(*this, RangePack::full_range(I, 1, rows()), RangePack::full_range(J,1,cols()));
890 }
891 
892 inline
893 const DMatrixSlice DMatrixSlice::operator()(const Range1D& I, const Range1D& J) const {
894  return DMatrixSlice( const_cast<DMatrixSlice&>(*this)
895  , RangePack::full_range(I, 1, rows()), RangePack::full_range(J,1,cols()) );
896 }
897 
898 inline
899 DMatrixSlice DMatrixSlice::operator()(size_type i1, size_type i2, size_type j1
900  , size_type j2)
901 {
902  return DMatrixSlice(*this, Range1D(i1,i2), Range1D(j1,j2));
903 }
904 
905 inline
906 const DMatrixSlice DMatrixSlice::operator()(size_type i1, size_type i2, size_type j1
907  , size_type j2) const
908 {
909  return DMatrixSlice( const_cast<DMatrixSlice&>(*this), Range1D(i1,i2)
910  , Range1D(j1,j2) );
911 }
912 
913 inline
915  return *this;
916 }
917 
918 inline
920  return *this;
921 }
922 
923 // Assignment operators
924 
925 inline
927  assign(this, alpha);
928  return *this;
929 }
930 
931 inline
933  assign(this, rhs, BLAS_Cpp::no_trans);
934  return *this;
935 }
936 
937 // Raw data access
938 
939 inline
940 DMatrixSlice::size_type DMatrixSlice::max_rows() const
941 { return max_rows_; }
942 
943 inline
944 DMatrixSlice::value_type* DMatrixSlice::col_ptr(size_type j) {
945  if( ptr_ )
946  validate_col_subscript(j);
947  return ptr_ + (j-1) * max_rows(); // will be 0 if not bound to a view.
948 }
949 
950 inline
951 const DMatrixSlice::value_type* DMatrixSlice::col_ptr(size_type j) const {
952  if( ptr_ )
953  validate_col_subscript(j);
954  return ptr_ + (j-1) * max_rows(); // will be 0 if not bound to a view.
955 }
956 
957 // ////////////////////////////////////////////////////////////////////////////////////////
958 // DMatrix inline member function definitions
959 
960 // Private utilities
961 
962 #ifndef LINALGPACK_CHECK_RANGE
963 inline
964 void DMatrix::validate_row_subscript(size_type i) const
965 {}
966 #endif
967 
968 #ifndef LINALGPACK_CHECK_RANGE
969 inline
970 void DMatrix::validate_col_subscript(size_type j) const
971 {}
972 #endif
973 
974 // constructors
975 
976 inline
977 DMatrix::DMatrix() : v_(), rows_(0)
978 {}
979 
980 inline
981 DMatrix::DMatrix(size_type rows, size_type cols)
982  : v_(rows*cols), rows_(rows)
983 {}
984 
985 inline
986 DMatrix::DMatrix(value_type val, size_type rows, size_type cols)
987  : v_(val,rows*cols), rows_(rows)
988 {}
989 
990 inline
991 DMatrix::DMatrix(const value_type* p, size_type rows, size_type cols)
992  : v_(rows*cols), rows_(rows)
993 {
994 // 6/7/00: valarray<> in libstdc++-2.90.7 has a bug in v_(p,size) so we do not
995 // use it. This is a hack until I can find the time to remove valarray all
996 // together.
997  std::copy( p, p + rows*cols, &v_[0] );
998 }
999 
1000 inline
1002  : v_(gms.rows() * gms.cols()), rows_(gms.rows())
1003 {
1004  assign(this, gms, BLAS_Cpp::no_trans);
1005 }
1006 
1007 // Memory management
1008 
1009 inline
1010 void DMatrix::resize(size_type rows, size_type cols, value_type val)
1011 {
1012  v_.resize(rows*cols,val);
1013  v_ = val;
1014  rows_ = rows;
1015 }
1016 
1017 inline
1019  v_.resize(0);
1020  rows_ = 0;
1021 }
1022 
1023 // Size / Dimensionality
1024 
1025 inline
1026 DMatrix::size_type DMatrix::rows() const {
1027  return rows_;
1028 }
1029 
1030 inline
1031 DMatrix::size_type DMatrix::cols() const {
1032  return rows_ > 0 ? v_.size() / rows_ : 0;
1033 }
1034 
1035 // Element access
1036 
1037 inline
1038 DMatrix::reference DMatrix::operator()(size_type i, size_type j)
1039 {
1040  validate_row_subscript(i); validate_col_subscript(j);
1041  return v_[(i-1) + (j-1) * rows_];
1042 }
1043 
1044 inline
1045 DMatrix::const_reference DMatrix::operator()(size_type i, size_type j) const
1046 {
1047  validate_row_subscript(i); validate_col_subscript(j);
1048  return (const_cast<std::valarray<value_type>&>(v_))[(i-1) + (j-1) * rows_];
1049 }
1050 
1051 // subregion access (range checked by constructors)
1052 
1053 inline
1055 {
1056  validate_row_subscript(i);
1057  return DVectorSlice( col_ptr(1) + (i-1), cols(), rows() );
1058 }
1059 
1060 inline
1061 const DVectorSlice DMatrix::row(size_type i) const
1062 {
1063  validate_row_subscript(i);
1064  return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + (i-1), cols(), rows() );
1065 }
1066 
1067 inline
1069 {
1070  validate_col_subscript(j);
1071  return DVectorSlice( col_ptr(1) + (j-1) * rows(), rows(), 1 );
1072 }
1073 
1074 inline
1075 const DVectorSlice DMatrix::col(size_type j) const
1076 {
1077  validate_col_subscript(j);
1078  return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + (j-1) * rows(), rows(), 1 ) ;
1079 }
1080 
1081 inline
1082 DVectorSlice DMatrix::diag(difference_type k)
1083 {
1084  return p_diag(k);
1085 }
1086 
1087 inline
1088 const DVectorSlice DMatrix::diag(difference_type k) const
1089 {
1090  return p_diag(k);
1091 }
1092 
1093 inline
1094 DMatrixSlice DMatrix::operator()(const Range1D& I, const Range1D& J)
1095 {
1096  Range1D Ix = RangePack::full_range(I,1,rows()), Jx = RangePack::full_range(J,1,cols());
1097  return DMatrixSlice( col_ptr(1) + (Ix.lbound() - 1) + (Jx.lbound() - 1) * rows()
1098  , max_rows() * cols(), max_rows(), Ix.size(), Jx.size() );
1099 }
1100 
1101 inline
1102 const DMatrixSlice DMatrix::operator()(const Range1D& I, const Range1D& J) const
1103 {
1104  Range1D Ix = RangePack::full_range(I,1,rows()), Jx = RangePack::full_range(J,1,cols());
1105  return DMatrixSlice( const_cast<value_type*>(col_ptr(1)) + (Ix.lbound() - 1) + (Jx.lbound() - 1) * rows()
1106  , max_rows() * cols(), max_rows(), Ix.size(), Jx.size() );
1107 }
1108 
1109 inline
1110 DMatrixSlice DMatrix::operator()(size_type i1, size_type i2, size_type j1
1111  , size_type j2)
1112 {
1113  return DMatrixSlice( col_ptr(1) + (i1 - 1) + (j1 - 1) * rows()
1114  , max_rows() * cols(), max_rows(), i2 - i1 + 1, j2 - j1 + 1 );
1115 }
1116 
1117 inline
1118 const DMatrixSlice DMatrix::operator()(size_type i1, size_type i2, size_type j1
1119  , size_type j2) const
1120 {
1121  return DMatrixSlice( const_cast<value_type*>(col_ptr(1)) + (i1 - 1) + (j1 - 1) * rows()
1122  , max_rows() * cols(), max_rows(), i2 - i1 + 1, j2 - j1 + 1 );
1123 }
1124 
1125 inline
1127 {
1128  return DMatrixSlice( col_ptr(1), max_rows() * cols(), max_rows(), rows(), cols() );
1129 }
1130 
1131 inline
1133 {
1134  return DMatrixSlice( const_cast<value_type*>(col_ptr(1)), max_rows() * cols(), max_rows()
1135  , rows(), cols() );
1136 }
1137 
1138 // Implicit conversion operators
1139 
1140 inline
1141 DMatrix::operator DMatrixSlice() {
1142  return (*this)();
1143 }
1144 
1145 inline
1146 DMatrix::operator const DMatrixSlice() const
1147 {
1148  return (*this)();
1149 }
1150 
1151 // Assignment operators
1152 
1153 inline
1154 DMatrix& DMatrix::operator=(value_type alpha)
1155 {
1156  assign(this, alpha);
1157  return *this;
1158 }
1159 
1160 inline
1162 {
1163  assign(this, rhs, BLAS_Cpp::no_trans);
1164  return *this;
1165 }
1166 
1167 inline
1169 {
1170  assign(this, rhs, BLAS_Cpp::no_trans);
1171  return *this;
1172 }
1173 
1174 // Raw data access
1175 
1176 inline
1177 DMatrix::size_type DMatrix::max_rows() const
1178 { return rows_; }
1179 
1180 inline
1181 DMatrix::value_type* DMatrix::col_ptr(size_type j)
1182 {
1183  if( v_.size() ) {
1184  validate_col_subscript(j);
1185  return &v_[ (j-1) * max_rows() ];
1186  }
1187  else {
1188  return 0;
1189  }
1190 }
1191 
1192 inline
1193 const DMatrix::value_type* DMatrix::col_ptr(size_type j) const
1194 {
1195  if( v_.size() ) {
1196  validate_col_subscript(j);
1197  return &const_cast<valarray&>(v_)[ (j-1) * max_rows() ];
1198  }
1199  else {
1200  return 0;
1201  }
1202 }
1203 
1204 } // end namespace DenseLinAlgPack
1205 
1206 #endif // GEN_MATRIX_CLASS_H
const DMatrixSlice * operator&() const
size_type cols() const
Return the number of columns.
DMatrixSlice * operator&()
Allow the address to be taken of an rvalue of this object.
Index size() const
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
DMatrixSlice & operator=(value_type alpha)
DMatrixSlice operator()()
Return a DMatrixSlice that represents this entire matrix.
size_type max_rows() const
Return the number of rows in the full matrix. Equivalent to BLAS LDA argument.
DVectorSlice row(size_type i)
Return DVectorSlice object representing the ith row (1-based; 1,2,..,#this->rows()#) ...
size_type rows() const
Return the number of rows.
size_t size_type
DVectorSlice col(size_type j)
Return DVectorSlice object representing the jth column (1-based; 1,2,..,#this->cols()#) ...
size_type max_rows() const
Return the number of rows in the full matrix. Equivalent to BLAS LDA argument.
void resize(size_type rows, size_type cols, value_type val=value_type())
Resize matrix to a (rows x cols) matrix and initializes any added elements by val.
void free()
frees memory and leaves a (0 x 0) matrix
DMatrix()
Construct a matrix with rows = cols = 0.
Index lbound() const
size_type rows() const
Return the number of rows.
DMatrixSlice & operator()()
Return reference of this. Included for iniformity with DMatrix.
DVectorSlice col(size_type j)
Return DVectorSlice object representing the jth column (1-based; 1,2,..,#this->cols()#, or throw std::out_of_range)
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
size_type cols() const
Return the number of columns.
DVectorSlice row(size_type i)
Return DVectorSlice object representing the ith row (1-based; 1,2,..,#this->rows()#, or throw std::out_of_range)