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

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

Modified istream constructors to support both the old format and CSV files. Also support has been added for NaN:s

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