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

Last change on this file since 808 was 808, checked in by Peter, 15 years ago

previous argument was invalid, but here is an implementation. Fixes #205

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