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

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

Fixes #138. Removed add, sub, scale and added appropriate operators.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.7 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 773 2007-03-01 00:51:25Z 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    /// 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    bool equal(const matrix&, const double precision=0) const;
167
168    ///
169    /// @return A const pointer to the internal GSL matrix.
170    ///
171    const gsl_matrix* gsl_matrix_p(void) const;
172
173    ///
174    /// @return A pointer to the internal GSL matrix.
175    ///
176    gsl_matrix* gsl_matrix_p(void);
177
178    ///
179    /// @return True if all elements in the matrix is zero, false
180    /// othwerwise;
181    ///
182    bool isnull(void) const;
183
184    ///
185    /// @brief Check if the matrix object is a view (sub-matrix) to
186    /// another matrix.
187    ///
188    /// @return True if the object is a view, false othwerwise.
189    ///
190    bool isview(void) const;
191
192    ///
193    /// @return The maximum value of the matrix.
194    ///
195    double max(void) const;
196
197    ///
198    /// @return The minimum value of the matrix.
199    ///
200    double min(void) const;
201
202    ///
203    /// @return The indecies to the element with the minimum and
204    /// maximum values of the matrix, respectively. The lowest index
205    /// has precedence (searching in row-major order).
206    ///
207    void minmax_index(std::pair<size_t,size_t>& min,
208                      std::pair<size_t,size_t>& max) const;
209
210    ///
211    /// The argument is changed into a matrix containing 1's and 0's
212    /// only. 1 means element in matrix is valid and a zero means
213    /// element is a NaN.
214    ///
215    /// @return true if matrix contains at least one NaN.
216    ///
217    bool nan(matrix&) const;
218
219    /**
220       Multiply the elements of matrix \a b with the elements of the
221       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
222       \f$. The result is stored into the calling matrix.
223
224       \throw GSL_error if dimensions mis-match.
225    */
226    void mul_elements(const matrix& b);
227
228    ///
229    /// @return The number of rows in the matrix.
230    ///
231    size_t rows(void) const;
232
233    /**
234       \brief Set element values to values in \a mat.
235
236       This function is needed for assignment of viewed elements.
237
238       \see const matrix& operator=(const matrix&)
239
240       \throw GSL_error if dimensions mis-match.
241    */
242    void set(const matrix& mat);
243
244    ///
245    /// Set all elements to \a value.
246    ///
247    void set_all(const double value);
248
249    /**
250       \brief Set \a column values to values in \a vec.
251
252       \note No check on size is done.
253
254       \throw GSL_error if index is out of range or mis-match in
255       sizes.
256    */
257    void set_column(const size_t column, const vector& vec);
258
259    /**
260       \brief Set \a row values to values in \a vec.
261
262       \note No check on size is done.
263
264       \throw GSL_error if index is out of range or mis-match in
265       sizes.
266    */
267    void set_row(const size_t row, const vector& vec);
268
269    /**
270       \brief Exchange the elements of the this and \a other by copying.
271
272       The two matrices must have the same size.
273
274       \throw GSL_error if either index is out of bounds.
275    */
276    void swap(matrix& other);
277
278    /**
279       \brief Swap columns \a i and \a j.
280
281       \throw GSL_error if either index is out of bounds.
282    */
283    void swap_columns(const size_t i,const size_t j);
284
285    /**
286       \brief Swap row \a i and column \a j.
287
288       The matrix must be square.
289
290       \throw GSL_error if either index is out of bounds, or if matrix
291       is not square.
292    */
293    void swap_rowcol(const size_t i,const size_t j);
294
295    /**
296       \brief Swap rows \a i and \a j.
297
298       \throw GSL_error if either index is out of bounds.
299    */
300    void swap_rows(const size_t i, const size_t j);
301
302    /**
303       \brief Transpose the matrix.
304
305       \throw GSL_error if memory allocation fails for the new
306       transposed matrix.
307    */
308    void transpose(void);
309
310    /**
311       \brief Element access operator.
312
313       \return Reference to the element position (\a row, \a column).
314
315       \throw If GSL range checks are enabled in the underlying GSL
316       library a GSL_error exception is thrown if either index is out
317       of range.
318    */
319    double& operator()(size_t row,size_t column);
320
321    /**
322       \brief Element access operator.
323
324       \return Const reference to the element position (\a row, \a
325       column).
326
327       \throw If GSL range checks are enabled in the underlying GSL
328       library a GSL_error exception is thrown if either index is out
329       of range.
330    */
331    const double& operator()(size_t row,size_t column) const;
332
333    ///
334    /// @brief Comparison operator.
335    ///
336    /// @return True if all elements are equal otherwise False.
337    ///
338    /// @see equal
339    ///
340    bool operator==(const matrix& other) const;
341
342    ///
343    /// @brief Comparison operator.
344    ///
345    /// @return False if all elements are equal otherwise True.
346    ///
347    /// @see equal
348    ///
349    bool operator!=(const matrix& other) const;
350
351    /**
352       \brief The assignment operator.
353
354       There is no requirements on dimensions, i.e. the matrix is
355       remapped in memory if necessary. This implies that in general
356       views cannot be assigned using this operator. Views will be
357       mutated into normal matrices. The only exception to this
358       behaviour on views is when self-assignemnt is done, since
359       self-assignment is ignored.
360
361       \return A const reference to the resulting matrix.
362
363       \see void set(const matrix&).
364
365       \throw A GSL_error is indirectly thrown if memory allocation
366       fails, or if dimensions mis-match.
367    */
368    const matrix& operator=(const matrix& other);
369
370    /**
371       \brief Add and assign operator.
372
373       Elementwise addition of the elements of matrix \a b to the left
374       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
375       \f$.
376
377       \return A const reference to the resulting matrix.
378
379       \throw GSL_error if dimensions mis-match.
380    */
381    const matrix& operator+=(const matrix& b);
382
383    /**
384       \brief Add and assign operator
385
386       Add the scalar value \a d to the left hand side matrix, \f$
387       a_{ij} = a_{ij} + d \; \forall i,j \f$.
388    */
389    const matrix& operator+=(const double d);
390
391    /**
392       \brief Subtract and assign operator.
393
394       Elementwise subtraction of the elements of matrix \a b to the
395       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
396       i,j \f$.
397
398       \return A const reference to the resulting matrix.
399
400       \throw GSL_error if dimensions mis-match.
401    */
402    const matrix& operator-=(const matrix&);
403
404    /**
405       \brief Subtract and assign operator
406
407       Subtract the scalar value \a d to the left hand side matrix,
408       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
409    */
410    const matrix& operator-=(const double d);
411
412    /**
413       \brief Multiply and assigment operator.
414
415       \return Const reference to the resulting matrix.
416
417       \throw GSL_error if memory allocation fails.
418    */
419    const matrix& operator*=(const matrix&);
420
421    /**
422       \brief Multiply and assignment operator
423
424       Multiply the elements of the left hand side matrix with a
425       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
426
427       \throw GSL_error if memory allocation fails.
428    */
429    const matrix& operator*=(double d);
430
431  private:
432
433    /**
434       \brief Create a new copy of the internal GSL matrix.
435
436       Necessary memory for the new GSL matrix is allocated and the
437       caller is responsible for freeing the allocated memory.
438
439       \return A pointer to a copy of the internal GSL matrix.
440
441       \throw GSL_error if memory cannot be allocated for the new
442       copy, or if dimensions mis-match.
443    */
444    gsl_matrix* create_gsl_matrix_copy(void) const;
445
446    // blas_result_ is used to temporarily store result in BLAS calls
447    // and when not NULL it should always have the same size as
448    // m_. Memory is not allocated for blas_result_ until it is
449    // needed.
450    gsl_matrix* blas_result_;
451    gsl_matrix* m_;
452    gsl_matrix_view* view_;
453  };
454
455  ///
456  /// The output operator for the matrix class.
457  ///
458  std::ostream& operator<< (std::ostream& s, const matrix&);
459
460  /**
461     @brief vector matrix multiplication
462   */
463  vector operator*(const matrix&, const vector&);
464
465  /**
466     @brief matrix vector multiplication
467   */
468  vector operator*(const vector&, const matrix&);
469
470}}} // of namespace utility, yat, and theplu
471
472#endif
Note: See TracBrowser for help on using the repository browser.