source: trunk/yat/utility/matrix.h @ 1064

Last change on this file since 1064 was 1064, checked in by Peter, 14 years ago

fixes #306 matrix iterators

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 1064 2008-02-10 22:06:52Z peter $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
10  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
11  Copyright (C) 2008 Peter Johansson
12
13  This file is part of the yat library, http://trac.thep.lu.se/yat
14
15  The yat library is free software; you can redistribute it and/or
16  modify it under the terms of the GNU General Public License as
17  published by the Free Software Foundation; either version 2 of the
18  License, or (at your option) any later version.
19
20  The yat library is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23  General Public License for more details.
24
25  You should have received a copy of the GNU General Public License
26  along with this program; if not, write to the Free Software
27  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
28  02111-1307, USA.
29*/
30
31#include "Exception.h"
32#include "StrideIterator.h"
33#include "vector.h"
34#include "VectorConstView.h"
35#include "VectorView.h"
36
37#include <gsl/gsl_matrix.h>
38#include <iostream>
39#include <utility>
40
41namespace theplu {
42namespace yat {
43namespace utility {
44
45  class VectorBase;
46
47  ///
48  /// @brief Interface to GSL matrix.
49  ///
50  /// For the time being 'double' is the only type supported.
51  ///
52  /// \par[File streams] Reading and writing vectors to file streams
53  /// are of course supported. These are implemented without using GSL
54  /// functionality, and thus binary read and write to streams are not
55  /// supported.
56  ///
57  /// \par[Matrix views] GSL matrix views are supported and these are
58  /// disguised as ordinary utility::matrix'es. A support function is
59  /// added, utility::matrix::isview(), that can be used to check if a
60  /// matrix object is a view. Note that view matricess do not own the
61  /// undelying data, and a view is not valid if the matrix owning the
62  /// data is deallocated.
63  ///
64  /// @note All GSL matrix related functions are not implement but
65  /// most functionality defined for GSL matrices can be achieved with
66  /// this interface class. Most notable GSL functionality not
67  /// supported are no binary file support and views on arrays,
68  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
69  /// superdiagonals.
70  ///
71  class matrix
72  {
73  public:
74    /**
75     */
76    typedef StrideIterator<double*> iterator;
77
78    /**
79     */
80    typedef StrideIterator<const double*> const_iterator;
81
82    /**
83       @brief The default constructor.
84   
85       This contructor does not initialize underlying (essential)
86       structures.
87    */
88    matrix(void);
89
90    /**
91       \brief Constructor allocating memory space for \a r times \a c
92       elements, and sets all elements to \a init_value.
93
94       \throw GSL_error if memory allocation fails.
95    */
96    matrix(const size_t& r, const size_t& c, double init_value=0);
97
98    /**
99       \brief The copy constructor.
100
101       \note If the object to be copied is a matrix view, the values
102       of the view will be copied, i.e. the view is not copied.
103
104       \throw A GSL_error is indirectly thrown if memory allocation
105       fails.
106    */
107    matrix(const matrix&);
108
109    /**
110       \brief The matrix view constructor.
111
112       Create a view of matrix \a m, with starting row \a offset_row,
113       starting column \a offset_column, row size \a n_row, and column
114       size \a n_column.
115
116       A matrix view can be used as any matrix with the difference
117       that changes made to the view will also change the object that
118       is viewed. Also, using the copy constructor will create a new
119       matrix object that is a copy of whatever is viewed. If a copy
120       of the view is needed then you should use this constructor to
121       obtain a copy.
122
123       \note If the object viewed by the view goes out of scope or is
124       deleted, the view becomes invalid and the result of further use
125       is undefined.
126
127       \throw GSL_error if a view cannot be set up.
128    */
129    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
130           size_t n_column);
131
132    /**
133       \brief The istream constructor.
134
135       Matrix rows are sepearated with the new line character, and
136       column elements in a row must be separated either with white
137       space characters except the new line (\\n) character (default),
138       or by the delimiter \a sep.  Rows, and columns, are read
139       consecutively and only rectangular matrices are
140       supported. Empty lines are ignored. End of input is at end of
141       file marker.
142
143       \throw GSL_error if memory allocation fails, IO_error if
144       unexpected input is found in the input stream.
145    */
146    explicit matrix(std::istream &, char sep='\0') 
147      throw(utility::IO_error, std::exception);
148
149    /**
150       \brief The destructor.
151    */
152    ~matrix(void);
153
154    ///
155    /// Set all elements to \a value.
156    ///
157    void all(const double value);
158
159    /**
160     */
161    iterator begin(void);
162
163    /**
164     */
165    const_iterator begin(void) const;
166
167    /**
168     */
169    iterator begin_column(size_t);
170
171    /**
172     */
173    const_iterator begin_column(size_t) const;
174
175    /**
176     */
177    iterator begin_row(size_t);
178
179    /**
180     */
181    const_iterator begin_row(size_t) const;
182
183    /**
184       \brief Make a copy of \a other.
185
186       This function will make a deep copy of \a other. Memory is
187       resized and view state is changed to same as in other.
188    */
189    const matrix& clone(const matrix& other);
190
191    /**
192     */
193    VectorView column_view(size_t);
194
195    /**
196     */
197    const VectorConstView column_const_view(size_t) const;
198
199    ///
200    /// @return The number of columns in the matrix.
201    ///
202    size_t columns(void) const;
203
204    /**
205       Elementwise division of the elemnts of the calling matrix by
206       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
207       \forall i,j \f$. The result is stored into the calling matrix.
208
209       \throw GSL_error if dimensions mis-match.
210    */
211    void div(const matrix& b);
212
213    /**
214     */
215    iterator end(void);
216
217    /**
218     */
219    const_iterator end(void) const;
220
221    /**
222     */
223    iterator end_column(size_t);
224
225    /**
226     */
227    const_iterator end_column(size_t) const;
228
229    /**
230     */
231    iterator end_row(size_t);
232
233    /**
234     */
235    const_iterator end_row(size_t) const;
236
237    /**
238       \brief Check whether matrices are equal within a user defined
239       precision, set by \a precision.
240
241       \return True if each element deviates less or equal than \a
242       d. If any matrix contain a NaN, false is always returned.
243
244       \see operator== and operator!=
245    */
246    bool equal(const matrix&, const double precision=0) const;
247
248    ///
249    /// @return A const pointer to the internal GSL matrix.
250    ///
251    const gsl_matrix* gsl_matrix_p(void) const;
252
253    ///
254    /// @return A pointer to the internal GSL matrix.
255    ///
256    gsl_matrix* gsl_matrix_p(void);
257
258    ///
259    /// @brief Check if the matrix object is a view (sub-matrix) to
260    /// another matrix.
261    ///
262    /// @return True if the object is a view, false othwerwise.
263    ///
264    bool isview(void) const;
265
266    /**
267       Multiply the elements of matrix \a b with the elements of the
268       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
269       \f$. The result is stored into the calling matrix.
270
271       \throw GSL_error if dimensions mis-match.
272    */
273    void mul(const matrix& b);
274
275    /**
276       \brief Resize matrix
277       
278       All elements are set to @a init_value.
279
280       \note underlying GSL matrix is destroyed and a view into this
281       matrix becomes invalid.
282    */
283    void resize(size_t, size_t, double init_value=0);
284
285    ///
286    /// @return The number of rows in the matrix.
287    ///
288    size_t rows(void) const;
289
290    /**
291     */
292    VectorView row_view(size_t);
293
294    /**
295     */
296    const VectorConstView row_const_view(size_t) const;
297
298    /**
299       \brief Swap columns \a i and \a j.
300
301       \throw GSL_error if either index is out of bounds.
302    */
303    void swap_columns(const size_t i,const size_t j);
304
305    /**
306       \brief Swap row \a i and column \a j.
307
308       The matrix must be square.
309
310       \throw GSL_error if either index is out of bounds, or if matrix
311       is not square.
312    */
313    void swap_rowcol(const size_t i,const size_t j);
314
315    /**
316       \brief Swap rows \a i and \a j.
317
318       \throw GSL_error if either index is out of bounds.
319    */
320    void swap_rows(const size_t i, const size_t j);
321
322    /**
323       \brief Transpose the matrix.
324
325       \throw GSL_error if memory allocation fails for the new
326       transposed matrix.
327    */
328    void transpose(void);
329
330    /**
331       \brief Element access operator.
332
333       \return Reference to the element position (\a row, \a column).
334
335       \throw If GSL range checks are enabled in the underlying GSL
336       library a GSL_error exception is thrown if either index is out
337       of range.
338    */
339    double& operator()(size_t row,size_t column);
340
341    /**
342       \brief Element access operator.
343
344       \return Const reference to the element position (\a row, \a
345       column).
346
347       \throw If GSL range checks are enabled in the underlying GSL
348       library a GSL_error exception is thrown if either index is out
349       of range.
350    */
351    const double& operator()(size_t row,size_t column) const;
352
353    /**
354       \brief Comparison operator. Takes squared time.
355
356       Checks are performed with exact matching, i.e., rounding off
357       effects may destroy comparison. Use the equal function for
358       comparing elements within a user defined precision.
359
360       \return True if all elements are equal otherwise false.
361
362       \see equal
363    */
364    bool operator==(const matrix& other) const;
365
366    /**
367       \brief Comparison operator. Takes squared time.
368
369       Checks are performed with exact matching, i.e., rounding off
370       effects may destroy comparison. Use the equal function for
371       comparing elements within a user defined precision.
372
373       \return False if all elements are equal otherwise true.
374
375       \see equal
376    */
377    bool operator!=(const matrix& other) const;
378
379    /**
380       \brief The assignment operator.
381
382       Dimensions of the matrices must match. If the LHS matrix is a
383       view, the underlying data will be changed.
384
385       \return A const reference to the resulting matrix.
386
387       \see void set(const matrix&).
388
389       \throw GSL_error if dimensions mis-match.
390    */
391    const matrix& operator=(const matrix& other);
392
393    /**
394       \brief Add and assign operator.
395
396       Elementwise addition of the elements of matrix \a b to the left
397       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
398       \f$.
399
400       \return A const reference to the resulting matrix.
401
402       \throw GSL_error if dimensions mis-match.
403    */
404    const matrix& operator+=(const matrix& b);
405
406    /**
407       \brief Add and assign operator
408
409       Add the scalar value \a d to the left hand side matrix, \f$
410       a_{ij} = a_{ij} + d \; \forall i,j \f$.
411    */
412    const matrix& operator+=(const double d);
413
414    /**
415       \brief Subtract and assign operator.
416
417       Elementwise subtraction of the elements of matrix \a b to the
418       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
419       i,j \f$.
420
421       \return A const reference to the resulting matrix.
422
423       \throw GSL_error if dimensions mis-match.
424    */
425    const matrix& operator-=(const matrix&);
426
427    /**
428       \brief Subtract and assign operator
429
430       Subtract the scalar value \a d to the left hand side matrix,
431       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
432    */
433    const matrix& operator-=(const double d);
434
435    /**
436       \brief Multiply and assigment operator.
437
438       \return Const reference to the resulting matrix.
439
440       \throw GSL_error if memory allocation fails.
441    */
442    const matrix& operator*=(const matrix&);
443
444    /**
445       \brief Multiply and assignment operator
446
447       Multiply the elements of the left hand side matrix with a
448       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
449
450       \throw GSL_error if memory allocation fails.
451    */
452    const matrix& operator*=(double d);
453
454  private:
455
456    /**
457       \brief Create a new copy of the internal GSL matrix.
458
459       Necessary memory for the new GSL matrix is allocated and the
460       caller is responsible for freeing the allocated memory.
461
462       \return A pointer to a copy of the internal GSL matrix.
463
464       \throw GSL_error if memory cannot be allocated for the new
465       copy, or if dimensions mis-match.
466    */
467    gsl_matrix* create_gsl_matrix_copy(void) const;
468
469    /**
470       \brief Clear all dynamically allocated memory.
471
472       Internal utility function.
473    */
474    void delete_allocated_memory(void);
475
476    // blas_result_ is used to temporarily store result in BLAS calls
477    // and when not NULL it should always have the same size as
478    // m_. Memory is not allocated for blas_result_ until it is
479    // needed.
480    gsl_matrix* blas_result_;
481    gsl_matrix* m_;
482    gsl_matrix_view* view_;
483    const gsl_matrix_const_view* view_const_;
484    // proxy_m_ is used to access the proper underlying gsl_matrix in
485    // all const member functions. It is not used by design for
486    // non-const vector functions and operators. This is to make sure
487    // that runtime errors occur if a const matrix is used in an
488    // inappropriate manner such as on left hand side in assignment
489    // (remember, m_ is null for const matrix views).
490    const gsl_matrix* proxy_m_;
491  };
492
493  /**
494     \brief Check if all elements of the matrix are zero.
495
496     \return True if all elements in the matrix is zero, false
497     othwerwise.
498  */
499  bool isnull(const matrix&);
500
501  /**
502     \brief Get the maximum value of the matrix.
503
504     \return The maximum value of the matrix.
505  */
506  double max(const matrix&);
507
508  /**
509     \brief Get the minimum value of the matrix.
510
511     \return The minimum value of the matrix.
512  */
513  double min(const matrix&);
514
515  /**
516     \brief Locate the maximum and minumum element in the matrix.
517
518     \return The indecies to the element with the minimum and maximum
519     values of the matrix, respectively.
520
521     \note Lower index has precedence (searching in row-major order).
522  */
523  void minmax_index(const matrix&,
524                    std::pair<size_t,size_t>& min,
525                    std::pair<size_t,size_t>& max);
526
527  /**
528     \brief Create a matrix \a flag indicating NaN's in another matrix
529     \a templat.
530
531     The \a flag matrix is changed to contain 1's and 0's only. A 1
532     means that the corresponding element in the \a templat matrix is
533     valid and a zero means that the corresponding element is a NaN.
534
535     \note Space for matrix \a flag is reallocated to fit the size of
536     matrix \a templat if sizes mismatch.
537
538     \return True if the \a templat matrix contains at least one NaN.
539  */
540  bool nan(const matrix& templat, matrix& flag);
541
542  /**
543     \brief Exchange all elements between the matrices by copying.
544
545     The two matrices must have the same size.
546
547     \throw GSL_error if either index is out of bounds.
548  */
549  void swap(matrix&, matrix&);
550
551  /**
552     \brief The output operator for the matrix class.
553  */
554  std::ostream& operator<< (std::ostream& s, const matrix&);
555
556  /**
557     \brief vector matrix multiplication
558   */
559  vector operator*(const matrix&, const VectorBase&);
560
561  /**
562     \brief matrix vector multiplication
563   */
564  vector operator*(const VectorBase&, const matrix&);
565
566}}} // of namespace utility, yat, and theplu
567
568#endif
Note: See TracBrowser for help on using the repository browser.