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

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

Addresses #436.

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