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

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

fixing Doxygen parsing

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