source: trunk/c++_tools/gslapi/matrix.h @ 611

Last change on this file since 611 was 611, checked in by Markus Ringnér, 15 years ago

Fixes #100

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.0 KB
Line 
1// $Id: matrix.h 611 2006-08-30 12:39:20Z markus $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
7  Copyright (C) 2006 Jari Häkkinen, Peter Johansson
8
9  This file is part of the thep c++ tools library,
10                                http://lev.thep.lu.se/trac/c++_tools
11
12  The c++ tools library is free software; you can redistribute it
13  and/or modify it under the terms of the GNU General Public License
14  as published by the Free Software Foundation; either version 2 of
15  the License, or (at your option) any later version.
16
17  The c++ tools library is distributed in the hope that it will be
18  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
19  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20  General Public License for more details.
21
22  You should have received a copy of the GNU General Public License
23  along with this program; if not, write to the Free Software
24  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
25  02111-1307, USA.
26*/
27
28#ifndef _theplu_gslapi_matrix_
29#define _theplu_gslapi_matrix_
30
31#include <c++_tools/gslapi/vector.h>
32#include <c++_tools/utility/Exception.h>
33
34#include <gsl/gsl_matrix.h>
35#include <iostream>
36
37namespace theplu {
38namespace gslapi {
39
40  class vector;
41
42  ///
43  /// This is the C++ tools 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 gslapi::matrix'es. A support function is
54  /// added, gslapi::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  /// gslapi::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    /// 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 &) throw (utility::IO_error,std::exception);
133    explicit matrix(std::istream &, char sep='\0') 
134      throw (utility::IO_error,std::exception);
135
136
137    ///
138    /// The destructor.
139    ///
140    ~matrix(void);
141
142    ///
143    /// Elementwise addition of the elements of matrix \a b to the
144    /// elements of the calling matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \;
145    /// \forall i,j \f$. The result is stored into the calling matrix.
146    ///
147    /// @return Whatever GSL returns.
148    ///
149    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
150
151    ///
152    /// Add the scalar value \a d to the elements of the calling
153    /// matrix, \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$. The result
154    /// is stored into the calling matrix.
155    ///
156    /// @return Whatever GSL returns.
157    ///
158    inline int
159    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
160
161    ///
162    /// @return The number of columns in the matrix.
163    ///
164    inline size_t columns(void) const { return m_->size2; }
165
166    ///
167    /// Elementwise division of the elemnts of the calling matrix by
168    /// the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
169    /// \forall i,j \f$. The result is stored into the calling matrix.
170    ///
171    /// @return Whatever GSL returns.
172    ///
173    inline int
174    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
175
176    ///
177    /// Check whether matrices are equal within a user defined
178    /// precision, set by \a precision.
179    ///
180    /// @return True if each element deviates less or equal than \a
181    /// d. If any matrix contain a NaN, false is always returned.
182    ///
183    bool equal(const matrix&, const double precision=0) const;
184
185    ///
186    /// @return A const pointer to the internal GSL matrix.
187    ///
188    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }
189
190    ///
191    /// @return A pointer to the internal GSL matrix.
192    ///
193    // Jari, is this needed?
194    inline gsl_matrix* gsl_matrix_p(void) { return m_; };
195
196    ///
197    /// @return True if all elements in the matrix is zero, false
198    /// othwerwise;
199    ///
200    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
201
202    ///
203    /// Check if the matrix object is a view (sub-matrix) to another
204    /// matrix.
205    ///
206    /// @return True if the object is a view, false othwerwise.
207    ///
208    inline bool isview(void) const { return view_; }
209
210    ///
211    /// @return The maximum value of the matrix.
212    ///
213    inline double max(void) const { return gsl_matrix_max(m_); }
214
215    ///
216    /// @return The index to the element with the maximum value of the
217    /// matrix. The lowest index has precedence (searching in
218    /// row-major order).
219    ///
220    inline std::pair<size_t,size_t> max_index(void) const;
221
222    ///
223    /// @return The minimum value of the matrix.
224    ///
225    inline double min(void) const { return gsl_matrix_min(m_); }
226
227    ///
228    /// @return The index to the element with the minimum value of the
229    /// matrix. The lowest index has precedence (searching in
230    /// row-major order).
231    ///
232    inline std::pair<size_t,size_t> min_index(void) const;
233
234    ///
235    /// @return The indecies to the element with the minimum and
236    /// maximum values of the matrix, respectively. The lowest index
237    /// has precedence (searching in row-major order).
238    ///
239    inline void minmax_index(std::pair<size_t,size_t>& min,
240                             std::pair<size_t,size_t>& max) const
241      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
242                                &max.first,&max.second); }
243
244    ///
245    /// @return The minimum and maximum values of the matrix, as the
246    /// \a first and \a second member of the returned \a pair,
247    /// respectively.
248    ///
249    std::pair<double,double> minmax(void) const;
250
251    ///
252    /// The argument is changed into a matrix containing 1's and 0's
253    /// only. 1 means element in matrix is valid and a zero means
254    /// element is a NaN.
255    ///
256    /// @return true if matrix contains at least one NaN.
257    ///
258    bool nan(matrix&) const;
259
260    ///
261    /// Multiply the elements of matrix \a b with the elements of the
262    /// calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall
263    /// i,j \f$. The result is stored into the calling matrix.
264    ///
265    /// @return Whatever GSL returns.
266    ///
267    inline int
268    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
269
270    ///
271    /// @return The number of rows in the matrix.
272    ///
273    inline size_t rows(void) const { return m_->size1; }
274
275    ///
276    /// Multiply the elements of the calling matrix with a scalar \a
277    /// d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$. The result is
278    /// stored into the calling matrix.
279    ///
280    /// @return Whatever GSL returns.
281    ///
282    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
283
284    ///
285    /// Set element values to values in \a mat. This function is
286    /// needed for assignment of viewed elements.
287    ///
288    /// @return Whatever GSL returns.
289    ///
290    /// @note No check on size is done.
291    ///
292    /// @see const matrix& operator=(const matrix&)
293    ///
294    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
295
296    ///
297    /// Set all elements to \a value.
298    ///
299    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
300
301    ///
302    /// Set \a column values to values in \a vec.
303    ///
304    /// @return Whatever GSL returns.
305    ///
306    /// @note No check on size is done.
307    ///
308    inline int set_column(const size_t column, const vector& vec)
309      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
310
311    ///
312    /// Set \a row values to values in \a vec.
313    ///
314    /// @return Whatever GSL returns.
315    ///
316    /// @note No check on size is done.
317    ///
318    inline int set_row(const size_t row, const vector& vec)
319      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
320
321    ///
322    /// Subtract the elements of matrix \a b from the elements of the
323    /// calling matrix ,\f$ a_{ij} = a_{ij} - b_{ij} \; \forall
324    /// i,j \f$. The result is stored into the calling matrix.
325    ///
326    /// @return Whatever GSL returns.
327    ///
328    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
329
330    ///
331    /// Exchange the elements of the this and \a other by copying. The
332    /// two matrices must have the same size.
333    ///
334    /// @return Whatever GSL returns.
335    ///
336    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
337
338    ///
339    /// Swap columns \a i and \a j.
340    ///
341    inline int swap_columns(const size_t i,const size_t j)
342      { return gsl_matrix_swap_columns(m_,i,j); }
343
344    ///
345    /// Swap row \a i and column \a j.
346    ///
347    inline int swap_rowcol(const size_t i,const size_t j)
348      { return gsl_matrix_swap_rowcol(m_,i,j); }
349
350    ///
351    /// Swap rows \a i and \a j.
352    ///
353    inline int swap_rows(const size_t i, const size_t j)
354      { return gsl_matrix_swap_rows(m_,i,j); }
355
356    ///
357    /// Transpose the matrix.
358    ///
359    void transpose(void);
360
361    ///
362    /// @return Reference to the element position (\a row, \a column).
363    ///
364    inline double& operator()(size_t row,size_t column)
365    { return (*gsl_matrix_ptr(m_,row,column)); }
366
367    ///
368    /// @return Const reference to the element position (\a row, \a
369    /// column).
370    ///
371    inline const double& operator()(size_t row,size_t column) const
372    { return (*gsl_matrix_const_ptr(m_,row,column)); }
373
374    ///
375    /// Matrix-vector multiplication.
376    ///
377    /// @return The resulting vector.
378    ///
379    // Jari, where should this go?
380    const vector operator*(const vector&) const;
381
382    ///
383    /// Comparison operator.
384    ///
385    /// @return True if all elements are equal otherwise False.
386    ///
387    /// @see equal
388    ///
389    inline bool operator==(const matrix& other) const { return equal(other); }
390
391    ///
392    /// Comparison operator.
393    ///
394    /// @return False if all elements are equal otherwise True.
395    ///
396    /// @see equal
397    ///
398    inline bool operator!=(const matrix& other) const { return !equal(other); }
399
400    ///
401    /// The assignment operator. There is no requirements on
402    /// dimensions, i.e. the matrix is remapped in memory if
403    /// necessary. This implies that in general views cannot be
404    /// assigned using this operator. Views will be mutated into
405    /// normal matrices. The only exception to this behaviour on views
406    /// is when self-assignemnt is done, since self-assignment is
407    /// ignored.
408    ///
409    /// @return A const reference to the resulting matrix.
410    ///
411    /// @see int set(const matrix&)
412    ///
413    const matrix& operator=(const matrix& other);
414
415    ///
416    /// Add and assign operator.
417    ///
418    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
419
420    ///
421    /// Add and assign operator.
422    ///
423    inline const matrix&
424    operator+=(const double d) { add_constant(d); return *this; }
425
426    ///
427    /// Subtract and assign operator.
428    ///
429    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
430
431    ///
432    /// Multiply and assigment operator.
433    ///
434    /// @return Const reference to the resulting matrix.
435    ///
436    const matrix& operator*=(const matrix&);
437
438    ///
439    /// Multiply and assign operator.
440    ///
441    inline const matrix& operator*=(const double d) { scale(d); return *this; }
442
443
444  private:
445
446    ///
447    /// Create a new copy of the internal GSL matrix.
448    ///
449    /// Necessary memory for the new GSL matrix is allocated and the
450    /// caller is responsible for freeing the allocated memory.
451    ///
452    /// @return A pointer to a copy of the internal GSL matrix.
453    ///
454    gsl_matrix* create_gsl_matrix_copy(void) const;
455
456    gsl_matrix* m_;
457    gsl_matrix_view* view_;
458  };
459
460  ///
461  /// The output operator for the matrix class.
462  ///
463  std::ostream& operator<< (std::ostream& s, const matrix&);
464
465
466}} // of namespace gslapi and namespace theplu
467
468#endif
Note: See TracBrowser for help on using the repository browser.