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

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

Merged better_matrix_class changes r402:419 into the trunk.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
Line 
1// $Id: matrix.h 420 2005-12-02 00:15:50Z jari $
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 with white space
100    /// characters except the new line (\\n) character. Rows, and
101    /// columns, are read consecutively and only rectangular matrices
102    /// are supported. Empty lines are ignored. End of input is at end
103    /// of file marker.
104    ///
105    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
106
107    ///
108    /// The destructor.
109    ///
110    ~matrix(void);
111
112    ///
113    /// Elementwise addition of the elements of matrix \a b to the
114    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
115    /// \forall i,j\f$. The result is stored into the calling matrix.
116    ///
117    /// @return Whatever GSL returns.
118    ///
119    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
120
121    ///
122    /// Add the scalar value \a d to the elements of the calling
123    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
124    /// is stored into the calling matrix.
125    ///
126    /// @return Whatever GSL returns.
127    ///
128    inline int
129    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
130
131    ///
132    /// @return The number of columns in the matrix.
133    ///
134    inline size_t columns(void) const { return m_->size2; }
135
136    ///
137    /// Elementwise division of the elemnts of the calling matrix by
138    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
139    /// \forall i,j\f$. The result is stored into the calling matrix.
140    ///
141    /// @return Whatever GSL returns.
142    ///
143    inline int
144    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
145
146    ///
147    /// Check whether matrices are equal within a user defined
148    /// precision, set by \a precision.
149    ///
150    /// @return True if each element deviates less or equal than \a
151    /// d. If any matrix contain a NaN, false is always returned.
152    ///
153    bool equal(const matrix&, const double precision=0) const;
154
155    ///
156    /// @return A const pointer to the internal GSL matrix.
157    ///
158    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }
159
160    ///
161    /// @return A pointer to the internal GSL matrix.
162    ///
163    // Jari, is this needed?
164    inline gsl_matrix* gsl_matrix_p(void) { return m_; };
165
166    ///
167    /// @return True if all elements in the matrix is zero, false
168    /// othwerwise;
169    ///
170    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
171
172    ///
173    /// Check if the matrix object is a view (sub-matrix) to another
174    /// matrix.
175    ///
176    /// @return True if the object is a view, false othwerwise.
177    ///
178    inline bool isview(void) const { return view_; }
179
180    ///
181    /// @return The maximum value of the matrix.
182    ///
183    inline double max(void) const { return gsl_matrix_max(m_); }
184
185    ///
186    /// @return The index to the element with the maximum value of the
187    /// matrix. The lowest index has precedence (searching in
188    /// row-major order).
189    ///
190    inline std::pair<size_t,size_t> max_index(void) const;
191
192    ///
193    /// @return The minimum value of the matrix.
194    ///
195    inline double min(void) const { return gsl_matrix_min(m_); }
196
197    ///
198    /// @return The index to the element with the minimum value of the
199    /// matrix. The lowest index has precedence (searching in
200    /// row-major order).
201    ///
202    inline std::pair<size_t,size_t> min_index(void) const;
203
204    ///
205    /// @return The indecies to the element with the minimum and
206    /// maximum values of the matrix, respectively. The lowest index
207    /// has precedence (searching in row-major order).
208    ///
209    inline void minmax_index(std::pair<size_t,size_t>& min,
210                             std::pair<size_t,size_t>& max) const
211      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
212                                &max.first,&max.second); }
213
214    ///
215    /// @return The minimum and maximum values of the matrix, as the
216    /// \a first and \a second member of the returned \a pair,
217    /// respectively.
218    ///
219    std::pair<double,double> minmax(void) const;
220
221    ///
222    /// Multiply the elements of matrix \a b with the elements of the
223    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
224    /// i,j\f$. The result is stored into the calling matrix.
225    ///
226    /// @return Whatever GSL returns.
227    ///
228    inline int
229    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
230
231    ///
232    /// @return The number of rows in the matrix.
233    ///
234    inline size_t rows(void) const { return m_->size1; }
235
236    ///
237    /// Multiply the elements of the calling matrix with a scalar \a
238    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
239    /// stored into the calling matrix.
240    ///
241    /// @return Whatever GSL returns.
242    ///
243    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
244
245    ///
246    /// Set element values to values in \a mat. This function is
247    /// needed for assignment of viewed elements.
248    ///
249    /// @return Whatever GSL returns.
250    ///
251    /// @note No check on size is done.
252    ///
253    /// @see const matrix& operator=(const matrix&)
254    ///
255    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
256
257    ///
258    /// Set all elements to \a value.
259    ///
260    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
261
262    ///
263    /// Set \a column values to values in \a vec.
264    ///
265    /// @return Whatever GSL returns.
266    ///
267    /// @note No check on size is done.
268    ///
269    inline int set_column(const size_t column, const vector& vec)
270      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
271
272    ///
273    /// Set \a row values to values in \a vec.
274    ///
275    /// @return Whatever GSL returns.
276    ///
277    /// @note No check on size is done.
278    ///
279    inline int set_row(const size_t row, const vector& vec)
280      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
281
282    ///
283    /// Subtract the elements of matrix \a b from the elements of the
284    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
285    /// i,j\f$. The result is stored into the calling matrix.
286    ///
287    /// @return Whatever GSL returns.
288    ///
289    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
290
291    ///
292    /// Exchange the elements of the this and \a other by copying. The
293    /// two matrices must have the same size.
294    ///
295    /// @return Whatever GSL returns.
296    ///
297    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
298
299    ///
300    /// Swap columns \a i and \a j.
301    ///
302    inline int swap_columns(const size_t i,const size_t j)
303      { return gsl_matrix_swap_columns(m_,i,j); }
304
305    ///
306    /// Swap row \a i and column \a j.
307    ///
308    inline int swap_rowcol(const size_t i,const size_t j)
309      { return gsl_matrix_swap_rowcol(m_,i,j); }
310
311    ///
312    /// Swap rows \a i and \a j.
313    ///
314    inline int swap_rows(const size_t i, const size_t j)
315      { return gsl_matrix_swap_rows(m_,i,j); }
316
317    ///
318    /// Transpose the matrix.
319    ///
320    void transpose(void);
321
322    ///
323    /// @return Reference to the element position (\a row, \a column).
324    ///
325    inline double& operator()(size_t row,size_t column)
326      { return (*gsl_matrix_ptr(m_,row,column)); }
327
328    ///
329    /// @return Const reference to the element position (\a row, \a
330    /// column).
331    ///
332    inline const double& operator()(size_t row,size_t column) const
333      { return (*gsl_matrix_const_ptr(m_,row,column)); }
334
335    ///
336    /// Matrix-vector multiplication.
337    ///
338    /// @return The resulting vector.
339    ///
340    // Jari, where should this go?
341    const vector operator*(const vector&) const;
342
343    ///
344    /// Comparison operator.
345    ///
346    /// @return True if all elements are equal otherwise False.
347    ///
348    /// @see equal
349    ///
350    inline bool operator==(const matrix& other) const { return equal(other); }
351
352    ///
353    /// Comparison operator.
354    ///
355    /// @return False if all elements are equal otherwise True.
356    ///
357    /// @see equal
358    ///
359    inline bool operator!=(const matrix& other) const { return !equal(other); }
360
361    ///
362    /// The assignment operator. There is no requirements on
363    /// dimensions, i.e. the matrix is remapped in memory if
364    /// necessary. This implies that in general views cannot be
365    /// assigned using this operator. Views will be mutated into
366    /// normal matrices. The only exception to this behaviour on views
367    /// is when self-assignemnt is done, since self-assignment is
368    /// ignored.
369    ///
370    /// @return A const reference to the resulting matrix.
371    ///
372    /// @see int set(const matrix&)
373    ///
374    const matrix& operator=(const matrix& other);
375
376    ///
377    /// Add and assign operator.
378    ///
379    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
380
381    ///
382    /// Add and assign operator.
383    ///
384    inline const matrix&
385    operator+=(const double d) { add_constant(d); return *this; }
386
387    ///
388    /// Subtract and assign operator.
389    ///
390    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
391
392    ///
393    /// Multiply and assigment operator.
394    ///
395    /// @return Const reference to the resulting matrix.
396    ///
397    const matrix& operator*=(const matrix&);
398
399    ///
400    /// Multiply and assign operator.
401    ///
402    inline const matrix& operator*=(const double d) { scale(d); return *this; }
403
404
405  private:
406
407    ///
408    /// Create a new copy of the internal GSL matrix.
409    ///
410    /// Necessary memory for the new GSL matrix is allocated and the
411    /// caller is responsible for freeing the allocated memory.
412    ///
413    /// @return A pointer to a copy of the internal GSL matrix.
414    ///
415    gsl_matrix* create_gsl_matrix_copy(void) const;
416
417    gsl_matrix* m_;
418    gsl_matrix_view* view_;
419  };
420
421  ///
422  /// The output operator for the matrix class.
423  ///
424  std::ostream& operator<< (std::ostream& s, const matrix&);
425
426
427}} // of namespace gslapi and namespace theplu
428
429#endif
Note: See TracBrowser for help on using the repository browser.