source: branches/better_matrix_class/lib/gslapi/matrix.h @ 415

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

Removed some gslapi copy intensive operators.
Made gslapi assignment operators ignore views, and change them
into normal vectors if assigned.
Cleaning up of code aswell.

  • 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 415 2005-12-01 15:52:36Z 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.