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

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

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