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

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

Addresses #193. matrix now works as outlined in the ticket
discussion. Added support for const views. Added a clone function that
facilitates resizing of matrices. clone is needed since assignement
operator functionality is changed.

  • 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 792 2007-03-11 00:14:24Z 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.
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 if needed.
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    /// @return The number of rows in the matrix.
206    ///
207    size_t rows(void) const;
208
209    /**
210       \brief Set element values to values in \a mat.
211
212       This function is needed for assignment of viewed elements.
213
214       \see const matrix& operator=(const matrix&)
215
216       \throw GSL_error if dimensions mis-match.
217    */
218    void set(const matrix& mat);
219
220    ///
221    /// Set all elements to \a value.
222    ///
223    void set_all(const double value);
224
225    /**
226       \brief Set \a column values to values in \a vec.
227
228       \note No check on size is done.
229
230       \throw GSL_error if index is out of range or mis-match in
231       sizes.
232    */
233    void set_column(const size_t column, const vector& vec);
234
235    /**
236       \brief Set \a row values to values in \a vec.
237
238       \note No check on size is done.
239
240       \throw GSL_error if index is out of range or mis-match in
241       sizes.
242    */
243    void set_row(const size_t row, const vector& vec);
244
245    /**
246       \brief Swap columns \a i and \a j.
247
248       \throw GSL_error if either index is out of bounds.
249    */
250    void swap_columns(const size_t i,const size_t j);
251
252    /**
253       \brief Swap row \a i and column \a j.
254
255       The matrix must be square.
256
257       \throw GSL_error if either index is out of bounds, or if matrix
258       is not square.
259    */
260    void swap_rowcol(const size_t i,const size_t j);
261
262    /**
263       \brief Swap rows \a i and \a j.
264
265       \throw GSL_error if either index is out of bounds.
266    */
267    void swap_rows(const size_t i, const size_t j);
268
269    /**
270       \brief Transpose the matrix.
271
272       \throw GSL_error if memory allocation fails for the new
273       transposed matrix.
274    */
275    void transpose(void);
276
277    /**
278       \brief Element access operator.
279
280       \return Reference to the element position (\a row, \a column).
281
282       \throw If GSL range checks are enabled in the underlying GSL
283       library a GSL_error exception is thrown if either index is out
284       of range.
285    */
286    double& operator()(size_t row,size_t column);
287
288    /**
289       \brief Element access operator.
290
291       \return Const reference to the element position (\a row, \a
292       column).
293
294       \throw If GSL range checks are enabled in the underlying GSL
295       library a GSL_error exception is thrown if either index is out
296       of range.
297    */
298    const double& operator()(size_t row,size_t column) const;
299
300    /**
301       \brief Comparison operator. Takes squared time.
302
303       Checks are performed with exact matching, i.e., rounding off
304       effects may destroy comparison. Use the equal function for
305       comparing elements within a user defined precision.
306
307       \return True if all elements are equal otherwise false.
308
309       \see equal
310    */
311    bool operator==(const matrix& other) const;
312
313    /**
314       \brief Comparison operator. Takes squared time.
315
316       Checks are performed with exact matching, i.e., rounding off
317       effects may destroy comparison. Use the equal function for
318       comparing elements within a user defined precision.
319
320       \return False if all elements are equal otherwise true.
321
322       \see equal
323    */
324    bool operator!=(const matrix& other) const;
325
326    /**
327       \brief The assignment operator.
328
329       Dimensions of the matrices must match. If the LHS matrix is a
330       view, the underlying data will be changed.
331
332       \return A const reference to the resulting matrix.
333
334       \see void set(const matrix&).
335
336       \throw GSL_error if dimensions mis-match.
337    */
338    const matrix& operator=(const matrix& other);
339
340    /**
341       \brief Add and assign operator.
342
343       Elementwise addition of the elements of matrix \a b to the left
344       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
345       \f$.
346
347       \return A const reference to the resulting matrix.
348
349       \throw GSL_error if dimensions mis-match.
350    */
351    const matrix& operator+=(const matrix& b);
352
353    /**
354       \brief Add and assign operator
355
356       Add the scalar value \a d to the left hand side matrix, \f$
357       a_{ij} = a_{ij} + d \; \forall i,j \f$.
358    */
359    const matrix& operator+=(const double d);
360
361    /**
362       \brief Subtract and assign operator.
363
364       Elementwise subtraction of the elements of matrix \a b to the
365       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
366       i,j \f$.
367
368       \return A const reference to the resulting matrix.
369
370       \throw GSL_error if dimensions mis-match.
371    */
372    const matrix& operator-=(const matrix&);
373
374    /**
375       \brief Subtract and assign operator
376
377       Subtract the scalar value \a d to the left hand side matrix,
378       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
379    */
380    const matrix& operator-=(const double d);
381
382    /**
383       \brief Multiply and assigment operator.
384
385       \return Const reference to the resulting matrix.
386
387       \throw GSL_error if memory allocation fails.
388    */
389    const matrix& operator*=(const matrix&);
390
391    /**
392       \brief Multiply and assignment operator
393
394       Multiply the elements of the left hand side matrix with a
395       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
396
397       \throw GSL_error if memory allocation fails.
398    */
399    const matrix& operator*=(double d);
400
401  private:
402
403    /**
404       \brief Create a new copy of the internal GSL matrix.
405
406       Necessary memory for the new GSL matrix is allocated and the
407       caller is responsible for freeing the allocated memory.
408
409       \return A pointer to a copy of the internal GSL matrix.
410
411       \throw GSL_error if memory cannot be allocated for the new
412       copy, or if dimensions mis-match.
413    */
414    gsl_matrix* create_gsl_matrix_copy(void) const;
415
416    // blas_result_ is used to temporarily store result in BLAS calls
417    // and when not NULL it should always have the same size as
418    // m_. Memory is not allocated for blas_result_ until it is
419    // needed.
420    gsl_matrix* blas_result_;
421    gsl_matrix* m_;
422    gsl_matrix_view* view_;
423    const gsl_matrix_const_view* view_const_;
424    // proxy_m_ is used to access the proper underlying gsl_matrix in
425    // all const member functions. It is not used by design for
426    // non-const vector functions and operators. This is to make sure
427    // that runtime errors occur if a const matrix is used in an
428    // inappropriate manner such as on left hand side in assignment
429    // (remember, m_ is null for const matrix views).
430    const gsl_matrix* proxy_m_;
431  };
432
433  /**
434     \brief Check if all elements of the matrix are zero.
435
436     \return True if all elements in the matrix is zero, false
437     othwerwise.
438  */
439  bool isnull(const matrix&);
440
441  /**
442     \brief Get the maximum value of the matrix.
443
444     \return The maximum value of the matrix.
445  */
446  double max(const matrix&);
447
448  /**
449     \brief Get the minimum value of the matrix.
450
451     \return The minimum value of the matrix.
452  */
453  double min(const matrix&);
454
455  /**
456     \brief Locate the maximum and minumum element in the matrix.
457
458     \return The indecies to the element with the minimum and maximum
459     values of the matrix, respectively.
460
461     \note Lower index has precedence (searching in row-major order).
462  */
463  void minmax_index(const matrix&,
464                    std::pair<size_t,size_t>& min,
465                    std::pair<size_t,size_t>& max);
466
467  /**
468     \brief Create a matrix \a flag indicating NaN's in another matrix
469     \a templat.
470
471     The \a flag matrix is changed to contain 1's and 0's only. A 1
472     means that the corresponding element in the \a templat matrix is
473     valid and a zero means that the corresponding element is a NaN.
474
475     \note Space for matrix \a flag is reallocated to fit the size of
476     matrix \a templat if sizes mismatch.
477
478     \return True if the \a templat matrix contains at least one NaN.
479  */
480  bool nan(const matrix& templat, matrix& flag);
481
482  /**
483     \brief Exchange all elements between the matrices by copying.
484
485     The two matrices must have the same size.
486
487     \throw GSL_error if either index is out of bounds.
488  */
489  void swap(matrix&, matrix&);
490
491  /**
492     \brief The output operator for the matrix class.
493  */
494  std::ostream& operator<< (std::ostream& s, const matrix&);
495
496  /**
497     \brief vector matrix multiplication
498   */
499  vector operator*(const matrix&, const vector&);
500
501  /**
502     \brief matrix vector multiplication
503   */
504  vector operator*(const vector&, const matrix&);
505
506}}} // of namespace utility, yat, and theplu
507
508#endif
Note: See TracBrowser for help on using the repository browser.