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

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

Cleaning.

  • 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 684 2006-10-13 17:35:49Z 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
12  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
13
14  The yat library is free software; you can redistribute it and/or
15  modify it under the terms of the GNU General Public License as
16  published by the Free Software Foundation; either version 2 of the
17  License, or (at your option) any later version.
18
19  The yat library is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22  General Public License for more details.
23
24  You should have received a copy of the GNU General Public License
25  along with this program; if not, write to the Free Software
26  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
27  02111-1307, USA.
28*/
29
30#include "vector.h"
31#include "Exception.h"
32
33#include <gsl/gsl_matrix.h>
34#include <iostream>
35
36namespace theplu {
37namespace yat {
38namespace utility {
39
40  class vector;
41
42  ///
43  /// This is the yat interface to GSL matrix. 'double' is the
44  /// only type supported, maybe we should add a 'complex' type as
45  /// well in the future.
46  ///
47  /// \par[File streams] Reading and writing vectors to file streams
48  /// are of course supported. These are implemented without using GSL
49  /// functionality, and thus binary read and write to streams are not
50  /// supported.
51  ///
52  /// \par[Matrix views] GSL matrix views are supported and these are
53  /// disguised as ordinary utility::matrix'es. A support function is
54  /// added, utility::matrix::isview(), that can be used to check if a
55  /// matrix object is a view. Note that view matricess do not own the
56  /// undelying data, and a view is not valid if the matrix owning the
57  /// data is deallocated.
58  ///
59  /// @note All GSL matrix related functions are not implement but
60  /// most functionality defined for GSL matrices can be achieved with
61  /// this interface class. Most notable GSL functionality not
62  /// supported are no binary file support and views on arrays,
63  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
64  /// superdiagonals. If there is an interest from the user community,
65  /// the omitted functionality could be included.
66  ///
67  /// @todo Maybe it would be smart to create temporary objects need
68  /// for BLAS in constructors. As it is now temporary objects are
69  /// called before BLAS functionality iss used, cf. const matrix&
70  /// matrix::operator*=(const matrix& other)
71  ///
72  class matrix
73  {
74  public:
75
76    /**
77       @brief The default constructor.
78   
79       This contructor does not initialize underlying (essential)
80       structures.
81    */
82    inline matrix(void) : m_(NULL), view_(NULL) {}
83
84    ///
85    /// Constructor. Allocates memory space for \a r times \a c
86    /// elements, and sets all elements to \a init_value.
87    ///
88    inline matrix(const size_t& r, const size_t& c, double init_value=0)
89      : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); }
90
91    ///
92    /// The copy constructor.
93    ///
94    /// @note If the object to be copied is a matrix view, the values
95    /// of the view will be copied, i.e. the view is not copied.
96    ///
97    inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();}
98
99    ///
100    /// 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
104    /// column 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
115    /// use is undefined.
116    ///
117    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
118           size_t n_column);
119
120    ///
121    /// The istream constructor.
122    ///
123    /// Matrix rows are sepearated with the new line character, and
124    /// column elements in a row must be separated either with white space
125    /// characters except the new line (\\n) character (default), or
126    /// by the delimiter \a sep.
127    ///. Rows, and
128    /// columns, are read consecutively and only rectangular matrices
129    /// are supported. Empty lines are ignored. End of input is at end
130    /// of file marker.
131    ///
132    explicit matrix(std::istream &, char sep='\0') 
133      throw(utility::IO_error, std::exception);
134
135    ///
136    /// The destructor.
137    ///
138    ~matrix(void);
139
140    ///
141    /// Elementwise addition of the elements of matrix \a b to the
142    /// elements of the calling matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \;
143    /// \forall i,j \f$. The result is stored into the calling matrix.
144    ///
145    /// @return Whatever GSL returns.
146    ///
147    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
148
149    ///
150    /// Add the scalar value \a d to the elements of the calling
151    /// matrix, \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$. The result
152    /// is stored into the calling matrix.
153    ///
154    /// @return Whatever GSL returns.
155    ///
156    inline int
157    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
158
159    ///
160    /// @return The number of columns in the matrix.
161    ///
162    inline size_t columns(void) const { return m_->size2; }
163
164    ///
165    /// Elementwise division of the elemnts of the calling matrix by
166    /// the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
167    /// \forall i,j \f$. The result is stored into the calling matrix.
168    ///
169    /// @return Whatever GSL returns.
170    ///
171    inline int
172    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
173
174    ///
175    /// Check whether matrices are equal within a user defined
176    /// precision, set by \a precision.
177    ///
178    /// @return True if each element deviates less or equal than \a
179    /// d. If any matrix contain a NaN, false is always returned.
180    ///
181    bool equal(const matrix&, const double precision=0) const;
182
183    ///
184    /// @return A const pointer to the internal GSL matrix.
185    ///
186    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }
187
188    ///
189    /// @return A pointer to the internal GSL matrix.
190    ///
191    // Jari, is this needed?
192    inline gsl_matrix* gsl_matrix_p(void) { return m_; };
193
194    ///
195    /// @return True if all elements in the matrix is zero, false
196    /// othwerwise;
197    ///
198    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
199
200    ///
201    /// Check if the matrix object is a view (sub-matrix) to another
202    /// matrix.
203    ///
204    /// @return True if the object is a view, false othwerwise.
205    ///
206    inline bool isview(void) const { return view_; }
207
208    ///
209    /// @return The maximum value of the matrix.
210    ///
211    inline double max(void) const { return gsl_matrix_max(m_); }
212
213    ///
214    /// @return The index to the element with the maximum value of the
215    /// matrix. The lowest index has precedence (searching in
216    /// row-major order).
217    ///
218    inline std::pair<size_t,size_t> max_index(void) const;
219
220    ///
221    /// @return The minimum value of the matrix.
222    ///
223    inline double min(void) const { return gsl_matrix_min(m_); }
224
225    ///
226    /// @return The index to the element with the minimum value of the
227    /// matrix. The lowest index has precedence (searching in
228    /// row-major order).
229    ///
230    inline std::pair<size_t,size_t> min_index(void) const;
231
232    ///
233    /// @return The indecies to the element with the minimum and
234    /// maximum values of the matrix, respectively. The lowest index
235    /// has precedence (searching in row-major order).
236    ///
237    inline void minmax_index(std::pair<size_t,size_t>& min,
238                             std::pair<size_t,size_t>& max) const
239      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
240                                &max.first,&max.second); }
241
242    ///
243    /// @return The minimum and maximum values of the matrix, as the
244    /// \a first and \a second member of the returned \a pair,
245    /// respectively.
246    ///
247    std::pair<double,double> minmax(void) const;
248
249    ///
250    /// The argument is changed into a matrix containing 1's and 0's
251    /// only. 1 means element in matrix is valid and a zero means
252    /// element is a NaN.
253    ///
254    /// @return true if matrix contains at least one NaN.
255    ///
256    bool nan(matrix&) const;
257
258    ///
259    /// Multiply the elements of matrix \a b with the elements of the
260    /// calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall
261    /// i,j \f$. The result is stored into the calling matrix.
262    ///
263    /// @return Whatever GSL returns.
264    ///
265    inline int
266    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
267
268    ///
269    /// @return The number of rows in the matrix.
270    ///
271    inline size_t rows(void) const { return m_->size1; }
272
273    ///
274    /// Multiply the elements of the calling matrix with a scalar \a
275    /// d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$. The result is
276    /// stored into the calling matrix.
277    ///
278    /// @return Whatever GSL returns.
279    ///
280    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
281
282    ///
283    /// Set element values to values in \a mat. This function is
284    /// needed for assignment of viewed elements.
285    ///
286    /// @return Whatever GSL returns.
287    ///
288    /// @note No check on size is done.
289    ///
290    /// @see const matrix& operator=(const matrix&)
291    ///
292    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
293
294    ///
295    /// Set all elements to \a value.
296    ///
297    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
298
299    ///
300    /// Set \a column values to values in \a vec.
301    ///
302    /// @return Whatever GSL returns.
303    ///
304    /// @note No check on size is done.
305    ///
306    inline int set_column(const size_t column, const vector& vec)
307      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
308
309    ///
310    /// Set \a row values to values in \a vec.
311    ///
312    /// @return Whatever GSL returns.
313    ///
314    /// @note No check on size is done.
315    ///
316    inline int set_row(const size_t row, const vector& vec)
317      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
318
319    ///
320    /// Subtract the elements of matrix \a b from the elements of the
321    /// calling matrix ,\f$ a_{ij} = a_{ij} - b_{ij} \; \forall
322    /// i,j \f$. The result is stored into the calling matrix.
323    ///
324    /// @return Whatever GSL returns.
325    ///
326    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
327
328    ///
329    /// Exchange the elements of the this and \a other by copying. The
330    /// two matrices must have the same size.
331    ///
332    /// @return Whatever GSL returns.
333    ///
334    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
335
336    ///
337    /// Swap columns \a i and \a j.
338    ///
339    inline int swap_columns(const size_t i,const size_t j)
340      { return gsl_matrix_swap_columns(m_,i,j); }
341
342    ///
343    /// Swap row \a i and column \a j.
344    ///
345    inline int swap_rowcol(const size_t i,const size_t j)
346      { return gsl_matrix_swap_rowcol(m_,i,j); }
347
348    ///
349    /// Swap rows \a i and \a j.
350    ///
351    inline int swap_rows(const size_t i, const size_t j)
352      { return gsl_matrix_swap_rows(m_,i,j); }
353
354    ///
355    /// Transpose the matrix.
356    ///
357    void transpose(void);
358
359    ///
360    /// @return Reference to the element position (\a row, \a column).
361    ///
362    inline double& operator()(size_t row,size_t column)
363    { return (*gsl_matrix_ptr(m_,row,column)); }
364
365    ///
366    /// @return Const reference to the element position (\a row, \a
367    /// column).
368    ///
369    inline const double& operator()(size_t row,size_t column) const
370    { return (*gsl_matrix_const_ptr(m_,row,column)); }
371
372    ///
373    /// Matrix-vector multiplication.
374    ///
375    /// @return The resulting vector.
376    ///
377    // Jari, where should this go?
378    const vector operator*(const vector&) const;
379
380    ///
381    /// Comparison operator.
382    ///
383    /// @return True if all elements are equal otherwise False.
384    ///
385    /// @see equal
386    ///
387    inline bool operator==(const matrix& other) const { return equal(other); }
388
389    ///
390    /// Comparison operator.
391    ///
392    /// @return False if all elements are equal otherwise True.
393    ///
394    /// @see equal
395    ///
396    inline bool operator!=(const matrix& other) const { return !equal(other); }
397
398    ///
399    /// The assignment operator. There is no requirements on
400    /// dimensions, i.e. the matrix is remapped in memory if
401    /// necessary. This implies that in general views cannot be
402    /// assigned using this operator. Views will be mutated into
403    /// normal matrices. The only exception to this behaviour on views
404    /// is when self-assignemnt is done, since self-assignment is
405    /// ignored.
406    ///
407    /// @return A const reference to the resulting matrix.
408    ///
409    /// @see int set(const matrix&)
410    ///
411    const matrix& operator=(const matrix& other);
412
413    ///
414    /// Add and assign operator.
415    ///
416    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
417
418    ///
419    /// Add and assign operator.
420    ///
421    inline const matrix&
422    operator+=(const double d) { add_constant(d); return *this; }
423
424    ///
425    /// Subtract and assign operator.
426    ///
427    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
428
429    ///
430    /// Multiply and assigment operator.
431    ///
432    /// @return Const reference to the resulting matrix.
433    ///
434    const matrix& operator*=(const matrix&);
435
436    ///
437    /// Multiply and assign operator.
438    ///
439    inline const matrix& operator*=(const double d) { scale(d); return *this; }
440
441
442  private:
443
444    ///
445    /// Create a new copy of the internal GSL matrix.
446    ///
447    /// Necessary memory for the new GSL matrix is allocated and the
448    /// caller is responsible for freeing the allocated memory.
449    ///
450    /// @return A pointer to a copy of the internal GSL matrix.
451    ///
452    gsl_matrix* create_gsl_matrix_copy(void) const;
453
454    gsl_matrix* m_;
455    gsl_matrix_view* view_;
456  };
457
458  ///
459  /// The output operator for the matrix class.
460  ///
461  std::ostream& operator<< (std::ostream& s, const matrix&);
462
463
464}}} // of namespace utility, yat and theplu
465
466#endif
Note: See TracBrowser for help on using the repository browser.