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

Last change on this file since 1890 was 1890, checked in by Peter, 12 years ago

fixes #454 - throw exception when detecting funky dimensions

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: Matrix.h 1890 2009-04-02 21:43:42Z 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, 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       If \a r is zero \a c must be zero and vice versa.
129
130       \throw GSL_error if memory allocation fails.
131    */
132    Matrix(const size_t& r, const size_t& c, double init_value=0);
133
134    /**
135       \brief The copy constructor.
136
137       \throw A GSL_error is indirectly thrown if memory allocation
138       fails.
139    */
140    Matrix(const Matrix&);
141
142    /**
143       \brief The istream constructor.
144
145       The std::istream will be interpreted as outlined here:
146
147       Missing values, i.e. empty elements, are treated as NaN values
148       (std::numeric_limits<double>::quiet_NaN() to be specific).
149
150       Matrix rows are separated with the new line character.
151
152       Column element separation has two modes depending on the value
153       of \a sep.
154
155       - If \a sep is the default '\\0' value then column elements are
156       separated with white space characters except the new line
157       character. Multiple sequential white space characters are
158       treated as one separator.
159
160       - Setting \a sep to something else than the default value will
161       change the behaviour to use the \a sep character as the
162       separator between column elements. Multiple sequential \a sep
163       characters will be treated as separating elements with missing
164       values.
165
166       End of input is the end of file marker and this treatment
167       cannot be redefined using the provided API.
168
169       \throw GSL_error if memory allocation fails, IO_error if
170       unexpected input is found in the input stream.
171    */
172    explicit Matrix(std::istream &, char sep='\0') 
173      throw(utility::IO_error, std::exception);
174
175    /**
176       \brief The destructor.
177    */
178    ~Matrix(void);
179
180    ///
181    /// Set all elements to \a value.
182    ///
183    void all(const double value);
184
185    /**
186       Iterator iterates along a row. When end of row is reached it
187       jumps to beginning of next row.
188
189       \return iterator pointing to upper-left element.
190     */
191    iterator begin(void);
192
193    /**
194       Iterator iterates along a row. When end of row is reached it
195       jumps to beginning of next row.
196
197       \return const_iterator pointing to upper-left element.
198     */
199    const_iterator begin(void) const;
200
201    /**
202       Iterator iterates along a column.
203
204       \return iterator pointing to first element of column \a i.
205     */
206    iterator begin_column(size_t i);
207
208    /**
209       Iterator iterates along a column.
210
211       \return const_iterator pointing to first element of column \a i.
212     */
213    const_iterator begin_column(size_t i) const;
214
215    /**
216       Iterator iterates along a row.
217
218       \return iterator pointing to first element of row \a i.
219     */
220    iterator begin_row(size_t i);
221
222    /**
223       Iterator iterates along a row.
224
225       \return const_iterator pointing to first element of row \a i.
226     */
227    const_iterator begin_row(size_t i) const;
228
229    /**
230       \return Vector view of column \a i
231     */
232    VectorView column_view(size_t i);
233
234    /**
235       \return const vector view of column \a i
236     */
237    const VectorConstView column_const_view(size_t) const;
238
239    ///
240    /// @return The number of columns in the matrix.
241    ///
242    size_t columns(void) const;
243
244    /**
245       Elementwise division of the elements of the calling matrix by
246       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
247       \forall i,j \f$. The result is stored into the calling matrix.
248
249       \throw GSL_error if dimensions mismatch.
250    */
251    void div(const Matrix& b);
252
253    /**
254       \return iterator pointing to end of matrix
255     */
256    iterator end(void);
257
258    /**
259       \return const_iterator pointing to end of matrix
260     */
261    const_iterator end(void) const;
262
263    /**
264       \return iterator pointing to end of column \a i
265     */
266    iterator end_column(size_t i);
267
268    /**
269       \return const_iterator pointing to end of column \a i
270     */
271    const_iterator end_column(size_t i) const;
272
273    /**
274       \return iterator pointing to end of row \a i
275     */
276    iterator end_row(size_t i);
277
278    /**
279       \return const_iterator pointing to end of row \a i
280     */
281    const_iterator end_row(size_t i) const;
282
283    /**
284       \brief Check whether matrices are equal within a user defined
285       precision, set by \a precision.
286
287       \return True if each element deviates less or equal than \a
288       d. If any matrix contain a NaN, false is always returned.
289
290       \see operator== and operator!=
291    */
292    bool equal(const Matrix&, const double precision=0) const;
293
294    ///
295    /// @return A const pointer to the internal GSL matrix.
296    ///
297    const gsl_matrix* gsl_matrix_p(void) const;
298
299    ///
300    /// @return A pointer to the internal GSL matrix.
301    ///
302    gsl_matrix* gsl_matrix_p(void);
303
304    /**
305       Multiply the elements of Matrix \a b with the elements of the
306       calling Matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
307       \f$. The result is stored into the calling Matrix.
308
309       \throw GSL_error if dimensions mismatch.
310    */
311    void mul(const Matrix& b);
312
313    /**
314       \brief Resize Matrix
315       
316       All elements are set to @a init_value.
317
318       If \a r is zero \a c must be zero and vice versa.
319
320       \note underlying GSL matrix is destroyed and views into this
321       Matrix becomes invalid.
322    */
323    void resize(size_t r, size_t c, double init_value=0);
324
325    ///
326    /// @return The number of rows in the matrix.
327    ///
328    size_t rows(void) const;
329
330    /**
331       \return Vector view of \a i
332     */
333    VectorView row_view(size_t);
334
335    /**
336       \return const vector view of \a i
337     */
338    const VectorConstView row_const_view(size_t) const;
339
340    /**
341       \brief Swap columns \a i and \a j.
342
343       \throw GSL_error if either index is out of bounds.
344    */
345    void swap_columns(const size_t i,const size_t j);
346
347    /**
348       \brief Swap row \a i and column \a j.
349
350       The Matrix must be square.
351
352       \throw GSL_error if either index is out of bounds, or if matrix
353       is not square.
354    */
355    void swap_rowcol(const size_t i,const size_t j);
356
357    /**
358       \brief Swap rows \a i and \a j.
359
360       \throw GSL_error if either index is out of bounds.
361    */
362    void swap_rows(const size_t i, const size_t j);
363
364    /**
365       \brief Transpose the matrix.
366
367       \throw GSL_error if memory allocation fails for the new
368       transposed matrix.
369    */
370    void transpose(void);
371
372    /**
373       \brief Element access operator.
374
375       \return Reference to the element position (\a row, \a column).
376
377       \throw If GSL range checks are enabled in the underlying GSL
378       library a GSL_error exception is thrown if either index is out
379       of range.
380    */
381    double& operator()(size_t row,size_t column);
382
383    /**
384       \brief Element access operator.
385
386       \return Const reference to the element position (\a row, \a
387       column).
388
389       \throw If GSL range checks are enabled in the underlying GSL
390       library a GSL_error exception is thrown if either index is out
391       of range.
392    */
393    const double& operator()(size_t row,size_t column) const;
394
395    /**
396       \brief Comparison operator. Takes squared time.
397
398       Checks are performed with exact matching, i.e., rounding off
399       effects may destroy comparison. Use the equal function for
400       comparing elements within a user defined precision.
401
402       \return True if all elements are equal otherwise false.
403
404       \see equal
405    */
406    bool operator==(const Matrix& other) const;
407
408    /**
409       \brief Comparison operator. Takes squared time.
410
411       Checks are performed with exact matching, i.e., rounding off
412       effects may destroy comparison. Use the equal function for
413       comparing elements within a user defined precision.
414
415       \return False if all elements are equal otherwise true.
416
417       \see equal
418    */
419    bool operator!=(const Matrix& other) const;
420
421    /**
422       \brief The assignment operator.
423
424       \return A const reference to the resulting Matrix.
425    */
426    const Matrix& operator=(const Matrix& other);
427
428    /**
429       \brief Add and assign operator.
430
431       Elementwise addition of the elements of Matrix \a b to the left
432       hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
433       \f$.
434
435       \return A const reference to the resulting Matrix.
436
437       \throw GSL_error if dimensions mismatch.
438    */
439    const Matrix& operator+=(const Matrix& b);
440
441    /**
442       \brief Add and assign operator
443
444       Add the scalar value \a d to the left hand side Matrix, \f$
445       a_{ij} = a_{ij} + d \; \forall i,j \f$.
446    */
447    const Matrix& operator+=(const double d);
448
449    /**
450       \brief Subtract and assign operator.
451
452       Elementwise subtraction of the elements of Matrix \a b to the
453       left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
454       i,j \f$.
455
456       \return A const reference to the resulting Matrix.
457
458       \throw GSL_error if dimensions mismatch.
459    */
460    const Matrix& operator-=(const Matrix&);
461
462    /**
463       \brief Subtract and assign operator
464
465       Subtract the scalar value \a d to the left hand side Matrix,
466       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
467    */
468    const Matrix& operator-=(const double d);
469
470    /**
471       \brief Multiply and assignment operator.
472
473       \return Const reference to the resulting Matrix.
474
475       \throw GSL_error if memory allocation fails.
476    */
477    const Matrix& operator*=(const Matrix&);
478
479    /**
480       \brief Multiply and assignment operator
481
482       Multiply the elements of the left hand side Matrix with a
483       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
484
485       \throw GSL_error if memory allocation fails.
486    */
487    const Matrix& operator*=(double d);
488
489  private:
490
491    /**
492       \brief Create a new copy of the internal GSL matrix.
493
494       Necessary memory for the new GSL matrix is allocated and the
495       caller is responsible for freeing the allocated memory.
496
497       \return A pointer to a copy of the internal GSL matrix.
498
499       \throw GSL_error if memory cannot be allocated for the new
500       copy, or if dimensions mismatch.
501    */
502    gsl_matrix* create_gsl_matrix_copy(void) const;
503
504    /**
505       \brief Clear all dynamically allocated memory.
506
507       Internal utility function.
508    */
509    void delete_allocated_memory(void);
510
511    // blas_result_ is used to temporarily store result in BLAS calls
512    // and when not NULL it should always have the same size as
513    // m_. Memory is not allocated for blas_result_ until it is
514    // needed.
515    gsl_matrix* blas_result_;
516    gsl_matrix* m_;
517  };
518
519  /**
520     \brief Check if all elements of the Matrix are zero.
521
522     \return True if all elements in the Matrix is zero, false
523     otherwise.
524
525     \relates Matrix
526  */
527  bool isnull(const Matrix&);
528
529  /**
530     \brief Get the maximum value of the Matrix.
531
532     \return The maximum value of the Matrix.
533
534     \relates Matrix
535  */
536  double max(const Matrix&);
537
538  /**
539     \brief Get the minimum value of the Matrix.
540
541     \return The minimum value of the Matrix.
542
543     \relates Matrix
544  */
545  double min(const Matrix&);
546
547  /**
548     \brief Locate the maximum and minimum element in the Matrix.
549
550     \return The indices to the element with the minimum and maximum
551     values of the Matrix, respectively.
552
553     \note Lower index has precedence (searching in row-major order).
554
555     \relates Matrix
556  */
557  void minmax_index(const Matrix&,
558                    std::pair<size_t,size_t>& min,
559                    std::pair<size_t,size_t>& max);
560
561  /**
562     \brief Create a Matrix \a flag indicating NaN's in another Matrix
563     \a templat.
564
565     The \a flag Matrix is changed to contain 1's and 0's only. A 1
566     means that the corresponding element in the \a templat Matrix is
567     valid and a zero means that the corresponding element is a NaN.
568
569     \note Space for Matrix \a flag is reallocated to fit the size of
570     Matrix \a templat if sizes mismatch.
571
572     \return True if the \a templat Matrix contains at least one NaN.
573
574     \relates Matrix
575  */
576  bool nan(const Matrix& templat, Matrix& flag);
577
578  /**
579     \brief Exchange all elements between the matrices by copying.
580
581     The two matrices must have the same size.
582
583     \throw GSL_error if sizes are not equal.
584
585     \relates Matrix
586  */
587  void swap(Matrix&, Matrix&);
588
589  /**
590     \brief The output operator for the Matrix class.
591
592     \relates Matrix
593  */
594  std::ostream& operator<< (std::ostream& s, const Matrix&);
595
596  /**
597     \brief Vector Matrix multiplication
598
599     \relates Matrix
600   */
601  Vector operator*(const Matrix&, const VectorBase&);
602
603  /**
604     \brief Matrix Vector multiplication
605
606     \relatesalso Matrix
607
608     \relatesalso VectorBase
609   */
610  Vector operator*(const VectorBase&, const Matrix&);
611
612}}} // of namespace utility, yat, and theplu
613
614#endif
Note: See TracBrowser for help on using the repository browser.