source: trunk/yat/utility/matrix.h @ 788

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

Improved comparison operator and equal function doxygen doc. Added for self comparison in equal.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 788 2007-03-10 19:59:51Z jari $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
10  Copyright (C) 2006 Jari Häkkinen, Peter Johansson
11  Copyright (C) 2007 Jari Häkkinen
12
13  This file is part of the yat library, http://lev.thep.lu.se/trac/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 "vector.h"
32#include "Exception.h"
33
34#include <gsl/gsl_matrix.h>
35#include <iostream>
36#include <utility>
37
38namespace theplu {
39namespace yat {
40namespace utility {
41
42  class vector;
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  /// \par[Matrix views] GSL matrix views are supported and these are
55  /// disguised as ordinary utility::matrix'es. A support function is
56  /// added, utility::matrix::isview(), that can be used to check if a
57  /// matrix object is a view. Note that view matricess do not own the
58  /// undelying data, and a view is not valid if the matrix owning the
59  /// data is deallocated.
60  ///
61  /// @note All GSL matrix related functions are not implement but
62  /// most functionality defined for GSL matrices can be achieved with
63  /// this interface class. Most notable GSL functionality not
64  /// supported are no binary file support and views on arrays,
65  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
66  /// superdiagonals. If there is an interest from the user community,
67  /// the omitted functionality could be included.
68  ///
69  class matrix
70  {
71  public:
72
73    /**
74       @brief The default constructor.
75   
76       This contructor does not initialize underlying (essential)
77       structures.
78    */
79    matrix(void);
80
81    /**
82       \brief Constructor allocating memory space for \a r times \a c
83       elements, and sets all elements to \a init_value.
84
85       \throw GSL_error if memory allocation fails.
86    */
87    matrix(const size_t& r, const size_t& c, double init_value=0);
88
89    /**
90       \brief The copy constructor.
91
92       \note If the object to be copied is a matrix view, the values
93       of the view will be copied, i.e. the view is not copied.
94
95       \throw A GSL_error is indirectly thrown if memory allocation
96       fails.
97    */
98    matrix(const matrix&);
99
100    /**
101       \brief The matrix view constructor.
102
103       Create a view of matrix \a m, with starting row \a offset_row,
104       starting column \a offset_column, row size \a n_row, and column
105       size \a n_column.
106
107       A matrix view can be used as any matrix with the difference
108       that changes made to the view will also change the object that
109       is viewed. Also, using the copy constructor will create a new
110       matrix object that is a copy of whatever is viewed. If a copy
111       of the view is needed then you should use this constructor to
112       obtain a copy.
113
114       \note If the object viewed by the view goes out of scope or is
115       deleted, the view becomes invalid and the result of further use
116       is undefined.
117
118       \throw GSL_error if a view cannot be set up.
119    */
120    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
121           size_t n_column);
122
123    /**
124       \brief The istream constructor.
125
126       Matrix rows are sepearated with the new line character, and
127       column elements in a row must be separated either with white
128       space characters except the new line (\\n) character (default),
129       or by the delimiter \a sep.  Rows, and columns, are read
130       consecutively and only rectangular matrices are
131       supported. Empty lines are ignored. End of input is at end of
132       file marker.
133
134       \throw GSL_error if memory allocation fails, IO_error if
135       unexpected input is found in the input stream.
136    */
137    explicit matrix(std::istream &, char sep='\0') 
138      throw(utility::IO_error, std::exception);
139
140    /**
141       \brief The destructor.
142    */
143    ~matrix(void);
144
145    ///
146    /// @return The number of columns in the matrix.
147    ///
148    size_t columns(void) const;
149
150    /**
151       Elementwise division of the elemnts of the calling matrix by
152       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
153       \forall i,j \f$. The result is stored into the calling matrix.
154
155       \throw GSL_error if dimensions mis-match.
156    */
157    void div_elements(const matrix& b);
158
159    /**
160       \brief Check whether matrices are equal within a user defined
161       precision, set by \a precision.
162
163       \return True if each element deviates less or equal than \a
164       d. If any matrix contain a NaN, false is always returned.
165
166       \see operator== and operator!=
167    */
168    bool equal(const matrix&, const double precision=0) const;
169
170    ///
171    /// @return A const pointer to the internal GSL matrix.
172    ///
173    const gsl_matrix* gsl_matrix_p(void) const;
174
175    ///
176    /// @return A pointer to the internal GSL matrix.
177    ///
178    gsl_matrix* gsl_matrix_p(void);
179
180    ///
181    /// @brief Check if the matrix object is a view (sub-matrix) to
182    /// another matrix.
183    ///
184    /// @return True if the object is a view, false othwerwise.
185    ///
186    bool isview(void) const;
187
188    /**
189       Multiply the elements of matrix \a b with the elements of the
190       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
191       \f$. The result is stored into the calling matrix.
192
193       \throw GSL_error if dimensions mis-match.
194    */
195    void mul_elements(const matrix& b);
196
197    ///
198    /// @return The number of rows in the matrix.
199    ///
200    size_t rows(void) const;
201
202    /**
203       \brief Set element values to values in \a mat.
204
205       This function is needed for assignment of viewed elements.
206
207       \see const matrix& operator=(const matrix&)
208
209       \throw GSL_error if dimensions mis-match.
210    */
211    void set(const matrix& mat);
212
213    ///
214    /// Set all elements to \a value.
215    ///
216    void set_all(const double value);
217
218    /**
219       \brief Set \a column values to values in \a vec.
220
221       \note No check on size is done.
222
223       \throw GSL_error if index is out of range or mis-match in
224       sizes.
225    */
226    void set_column(const size_t column, const vector& vec);
227
228    /**
229       \brief Set \a row values to values in \a vec.
230
231       \note No check on size is done.
232
233       \throw GSL_error if index is out of range or mis-match in
234       sizes.
235    */
236    void set_row(const size_t row, const vector& vec);
237
238    /**
239       \brief Swap columns \a i and \a j.
240
241       \throw GSL_error if either index is out of bounds.
242    */
243    void swap_columns(const size_t i,const size_t j);
244
245    /**
246       \brief Swap row \a i and column \a j.
247
248       The matrix must be square.
249
250       \throw GSL_error if either index is out of bounds, or if matrix
251       is not square.
252    */
253    void swap_rowcol(const size_t i,const size_t j);
254
255    /**
256       \brief Swap rows \a i and \a j.
257
258       \throw GSL_error if either index is out of bounds.
259    */
260    void swap_rows(const size_t i, const size_t j);
261
262    /**
263       \brief Transpose the matrix.
264
265       \throw GSL_error if memory allocation fails for the new
266       transposed matrix.
267    */
268    void transpose(void);
269
270    /**
271       \brief Element access operator.
272
273       \return Reference to the element position (\a row, \a column).
274
275       \throw If GSL range checks are enabled in the underlying GSL
276       library a GSL_error exception is thrown if either index is out
277       of range.
278    */
279    double& operator()(size_t row,size_t column);
280
281    /**
282       \brief Element access operator.
283
284       \return Const reference to the element position (\a row, \a
285       column).
286
287       \throw If GSL range checks are enabled in the underlying GSL
288       library a GSL_error exception is thrown if either index is out
289       of range.
290    */
291    const double& operator()(size_t row,size_t column) const;
292
293    /**
294       \brief Comparison operator. Takes squared time.
295
296       Checks are performed with exact matching, i.e., rounding off
297       effects may destroy comparison. Use the equal function for
298       comparing elements within a user defined precision.
299
300       \return True if all elements are equal otherwise false.
301
302       \see equal
303    */
304    bool operator==(const matrix& other) const;
305
306    /**
307       \brief Comparison operator. Takes squared time.
308
309       Checks are performed with exact matching, i.e., rounding off
310       effects may destroy comparison. Use the equal function for
311       comparing elements within a user defined precision.
312
313       \return False if all elements are equal otherwise true.
314
315       \see equal
316    */
317    bool operator!=(const matrix& other) const;
318
319    /**
320       \brief The assignment operator.
321
322       There is no requirements on dimensions, i.e. the matrix is
323       remapped in memory if necessary. This implies that in general
324       views cannot be assigned using this operator. Views will be
325       mutated into normal matrices. The only exception to this
326       behaviour on views is when self-assignemnt is done, since
327       self-assignment is ignored.
328
329       \return A const reference to the resulting matrix.
330
331       \see void set(const matrix&).
332
333       \throw A GSL_error is indirectly thrown if memory allocation
334       fails, or if dimensions mis-match.
335    */
336    const matrix& operator=(const matrix& other);
337
338    /**
339       \brief Add and assign operator.
340
341       Elementwise addition of the elements of matrix \a b to the left
342       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
343       \f$.
344
345       \return A const reference to the resulting matrix.
346
347       \throw GSL_error if dimensions mis-match.
348    */
349    const matrix& operator+=(const matrix& b);
350
351    /**
352       \brief Add and assign operator
353
354       Add the scalar value \a d to the left hand side matrix, \f$
355       a_{ij} = a_{ij} + d \; \forall i,j \f$.
356    */
357    const matrix& operator+=(const double d);
358
359    /**
360       \brief Subtract and assign operator.
361
362       Elementwise subtraction of the elements of matrix \a b to the
363       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
364       i,j \f$.
365
366       \return A const reference to the resulting matrix.
367
368       \throw GSL_error if dimensions mis-match.
369    */
370    const matrix& operator-=(const matrix&);
371
372    /**
373       \brief Subtract and assign operator
374
375       Subtract the scalar value \a d to the left hand side matrix,
376       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
377    */
378    const matrix& operator-=(const double d);
379
380    /**
381       \brief Multiply and assigment operator.
382
383       \return Const reference to the resulting matrix.
384
385       \throw GSL_error if memory allocation fails.
386    */
387    const matrix& operator*=(const matrix&);
388
389    /**
390       \brief Multiply and assignment operator
391
392       Multiply the elements of the left hand side matrix with a
393       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
394
395       \throw GSL_error if memory allocation fails.
396    */
397    const matrix& operator*=(double d);
398
399  private:
400
401    /**
402       \brief Create a new copy of the internal GSL matrix.
403
404       Necessary memory for the new GSL matrix is allocated and the
405       caller is responsible for freeing the allocated memory.
406
407       \return A pointer to a copy of the internal GSL matrix.
408
409       \throw GSL_error if memory cannot be allocated for the new
410       copy, or if dimensions mis-match.
411    */
412    gsl_matrix* create_gsl_matrix_copy(void) const;
413
414    // blas_result_ is used to temporarily store result in BLAS calls
415    // and when not NULL it should always have the same size as
416    // m_. Memory is not allocated for blas_result_ until it is
417    // needed.
418    gsl_matrix* blas_result_;
419    gsl_matrix* m_;
420    gsl_matrix_view* view_;
421  };
422
423  /**
424     \brief Check if all elements of the matrix are zero.
425
426     \return True if all elements in the matrix is zero, false
427     othwerwise.
428  */
429  bool isnull(const matrix&);
430
431  /**
432     \brief Get the maximum value of the matrix.
433
434     \return The maximum value of the matrix.
435  */
436  double max(const matrix&);
437
438  /**
439     \brief Get the minimum value of the matrix.
440
441     \return The minimum value of the matrix.
442  */
443  double min(const matrix&);
444
445  /**
446     \brief Locate the maximum and minumum element in the matrix.
447
448     \return The indecies to the element with the minimum and maximum
449     values of the matrix, respectively.
450
451     \note Lower index has precedence (searching in row-major order).
452  */
453  void minmax_index(const matrix&,
454                    std::pair<size_t,size_t>& min,
455                    std::pair<size_t,size_t>& max);
456
457  /**
458     \brief Create a matrix \a flag indicating NaN's in another matrix
459     \a templat.
460
461     The \a flag matrix is changed to contain 1's and 0's only. A 1
462     means that the corresponding element in the \a templat matrix is
463     valid and a zero means that the corresponding element is a NaN.
464
465     \note Space for matrix \a flag is reallocated to fit the size of
466     matrix \a templat if sizes mismatch.
467
468     \return True if the \a templat matrix contains at least one NaN.
469  */
470  bool nan(const matrix& templat, matrix& flag);
471
472  /**
473     \brief Exchange all elements between the matrices by copying.
474
475     The two matrices must have the same size.
476
477     \throw GSL_error if either index is out of bounds.
478  */
479  void swap(matrix&, matrix&);
480
481  /**
482     \brief The output operator for the matrix class.
483  */
484  std::ostream& operator<< (std::ostream& s, const matrix&);
485
486  /**
487     \brief vector matrix multiplication
488   */
489  vector operator*(const matrix&, const vector&);
490
491  /**
492     \brief matrix vector multiplication
493   */
494  vector operator*(const vector&, const matrix&);
495
496}}} // of namespace utility, yat, and theplu
497
498#endif
Note: See TracBrowser for help on using the repository browser.