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

Last change on this file since 610 was 597, checked in by Markus Ringnér, 16 years ago

Fixed comments so they pass without some of the complaits from doxygen. Have not looked at the actual contents of comments

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.9 KB
Line 
1// $Id: matrix.h 597 2006-08-28 13:03:54Z 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    /// @return matrix containing 1 and 0 only. 1 means element in
253    /// matrix is valid and a zero means element is a NaN.
254    ///
255    matrix nan(void) const;
256
257    ///
258    /// Multiply the elements of matrix \a b with the elements of the
259    /// calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall
260    /// i,j \f$. The result is stored into the calling matrix.
261    ///
262    /// @return Whatever GSL returns.
263    ///
264    inline int
265    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
266
267    ///
268    /// @return The number of rows in the matrix.
269    ///
270    inline size_t rows(void) const { return m_->size1; }
271
272    ///
273    /// Multiply the elements of the calling matrix with a scalar \a
274    /// d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$. The result is
275    /// stored into the calling matrix.
276    ///
277    /// @return Whatever GSL returns.
278    ///
279    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
280
281    ///
282    /// Set element values to values in \a mat. This function is
283    /// needed for assignment of viewed elements.
284    ///
285    /// @return Whatever GSL returns.
286    ///
287    /// @note No check on size is done.
288    ///
289    /// @see const matrix& operator=(const matrix&)
290    ///
291    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
292
293    ///
294    /// Set all elements to \a value.
295    ///
296    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
297
298    ///
299    /// Set \a column values to values in \a vec.
300    ///
301    /// @return Whatever GSL returns.
302    ///
303    /// @note No check on size is done.
304    ///
305    inline int set_column(const size_t column, const vector& vec)
306      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
307
308    ///
309    /// Set \a row values to values in \a vec.
310    ///
311    /// @return Whatever GSL returns.
312    ///
313    /// @note No check on size is done.
314    ///
315    inline int set_row(const size_t row, const vector& vec)
316      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
317
318    ///
319    /// Subtract the elements of matrix \a b from the elements of the
320    /// calling matrix ,\f$ a_{ij} = a_{ij} - b_{ij} \; \forall
321    /// i,j \f$. The result is stored into the calling matrix.
322    ///
323    /// @return Whatever GSL returns.
324    ///
325    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
326
327    ///
328    /// Exchange the elements of the this and \a other by copying. The
329    /// two matrices must have the same size.
330    ///
331    /// @return Whatever GSL returns.
332    ///
333    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
334
335    ///
336    /// Swap columns \a i and \a j.
337    ///
338    inline int swap_columns(const size_t i,const size_t j)
339      { return gsl_matrix_swap_columns(m_,i,j); }
340
341    ///
342    /// Swap row \a i and column \a j.
343    ///
344    inline int swap_rowcol(const size_t i,const size_t j)
345      { return gsl_matrix_swap_rowcol(m_,i,j); }
346
347    ///
348    /// Swap rows \a i and \a j.
349    ///
350    inline int swap_rows(const size_t i, const size_t j)
351      { return gsl_matrix_swap_rows(m_,i,j); }
352
353    ///
354    /// Transpose the matrix.
355    ///
356    void transpose(void);
357
358    ///
359    /// @return Reference to the element position (\a row, \a column).
360    ///
361    inline double& operator()(size_t row,size_t column)
362    { return (*gsl_matrix_ptr(m_,row,column)); }
363
364    ///
365    /// @return Const reference to the element position (\a row, \a
366    /// column).
367    ///
368    inline const double& operator()(size_t row,size_t column) const
369    { return (*gsl_matrix_const_ptr(m_,row,column)); }
370
371    ///
372    /// Matrix-vector multiplication.
373    ///
374    /// @return The resulting vector.
375    ///
376    // Jari, where should this go?
377    const vector operator*(const vector&) const;
378
379    ///
380    /// Comparison operator.
381    ///
382    /// @return True if all elements are equal otherwise False.
383    ///
384    /// @see equal
385    ///
386    inline bool operator==(const matrix& other) const { return equal(other); }
387
388    ///
389    /// Comparison operator.
390    ///
391    /// @return False if all elements are equal otherwise True.
392    ///
393    /// @see equal
394    ///
395    inline bool operator!=(const matrix& other) const { return !equal(other); }
396
397    ///
398    /// The assignment operator. There is no requirements on
399    /// dimensions, i.e. the matrix is remapped in memory if
400    /// necessary. This implies that in general views cannot be
401    /// assigned using this operator. Views will be mutated into
402    /// normal matrices. The only exception to this behaviour on views
403    /// is when self-assignemnt is done, since self-assignment is
404    /// ignored.
405    ///
406    /// @return A const reference to the resulting matrix.
407    ///
408    /// @see int set(const matrix&)
409    ///
410    const matrix& operator=(const matrix& other);
411
412    ///
413    /// Add and assign operator.
414    ///
415    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
416
417    ///
418    /// Add and assign operator.
419    ///
420    inline const matrix&
421    operator+=(const double d) { add_constant(d); return *this; }
422
423    ///
424    /// Subtract and assign operator.
425    ///
426    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
427
428    ///
429    /// Multiply and assigment operator.
430    ///
431    /// @return Const reference to the resulting matrix.
432    ///
433    const matrix& operator*=(const matrix&);
434
435    ///
436    /// Multiply and assign operator.
437    ///
438    inline const matrix& operator*=(const double d) { scale(d); return *this; }
439
440
441  private:
442
443    ///
444    /// Create a new copy of the internal GSL matrix.
445    ///
446    /// Necessary memory for the new GSL matrix is allocated and the
447    /// caller is responsible for freeing the allocated memory.
448    ///
449    /// @return A pointer to a copy of the internal GSL matrix.
450    ///
451    gsl_matrix* create_gsl_matrix_copy(void) const;
452
453    gsl_matrix* m_;
454    gsl_matrix_view* view_;
455  };
456
457  ///
458  /// The output operator for the matrix class.
459  ///
460  std::ostream& operator<< (std::ostream& s, const matrix&);
461
462
463}} // of namespace gslapi and namespace theplu
464
465#endif
Note: See TracBrowser for help on using the repository browser.