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

Last change on this file since 1547 was 1547, checked in by Peter, 13 years ago

refs #448

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