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

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

Fixes #206.

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