source: trunk/yat/utility/Matrix.h @ 1823

Last change on this file since 1823 was 1823, checked in by Jari Häkkinen, 13 years ago

Fixes #463. Improved documentation on istream constructor.

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