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

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

fixes #308

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