source: trunk/lib/gslapi/matrix.h @ 536

Last change on this file since 536 was 536, checked in by Peter, 16 years ago

fixed problem with ConsensusInputranker? - Now all tests do pass

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