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

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

More view fixing.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 KB
Line 
1// $Id: matrix.h 409 2005-11-28 00:15:41Z 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  class matrix
43  {
44  public:
45
46    ///
47    /// The default constructor.
48    ///
49    /// This contructor does not initialize underlying (essential)
50    /// structures.
51    ///
52    matrix(void);
53
54    ///
55    /// Constructor. Allocates memory space for \a r times \a c
56    /// elements, and sets all elements to \a init_value.
57    ///
58    matrix(const size_t& r, const size_t& c, const double init_value=0);
59
60    ///
61    /// The copy constructor.
62    ///
63    /// @note If the object to be copied is a matrix view, the values
64    /// of the view will be copied, i.e. the view is not copied.
65    ///
66    matrix(const matrix&);
67
68    ///
69    /// The matrix view constructor.
70    ///
71    /// Create a view of matrix \a m, with starting row \a offset_row,
72    /// starting column \a offset_column, row size \a n_row, and
73    /// column size \a n_column.
74    ///
75    /// A matrix view can be used as any matrix with the difference
76    /// that changes made to the view will also change the object that
77    /// is viewed. Also, using the copy constructor will create a new
78    /// matrix object that is a copy of whatever is viewed. If a copy
79    /// of the view is needed then you should use this constructor to
80    /// obtain a copy.
81    ///
82    /// @note If the object viewed by the view goes out of scope or is
83    /// deleted, the view becomes invalid and the result of further
84    /// use is undefined.
85    ///
86    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
87           size_t n_column);
88
89    ///
90    /// The istream constructor.
91    ///
92    /// Matrix rows are sepearated with the new line character, and
93    /// column elements in a row must be separated with white space
94    /// characters except the new line (\\n) character. Rows, and
95    /// columns, are read consecutively and only rectangular matrices
96    /// are supported. Empty lines are ignored. End of input is at end
97    /// of file marker.
98    ///
99    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
100
101    ///
102    /// The destructor.
103    ///
104    ~matrix(void);
105
106    ///
107    /// Elementwise addition of the elements of matrix \a b to the
108    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
109    /// \forall i,j\f$. The result is stored into the calling matrix.
110    ///
111    /// @return Whatever GSL returns.
112    ///
113    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
114
115    ///
116    /// Add the scalar value \a d to the elements of the calling
117    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
118    /// is stored into the calling matrix.
119    ///
120    /// @return Whatever GSL returns.
121    ///
122    inline int
123    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
124
125    ///
126    /// @return The number of columns in the matrix.
127    ///
128    inline size_t columns(void) const { return m_->size2; }
129
130    ///
131    /// Elementwise division of the elemnts of the calling matrix by
132    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
133    /// \forall i,j\f$. The result is stored into the calling matrix.
134    ///
135    /// @return Whatever GSL returns.
136    ///
137    inline int
138    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
139
140    ///
141    /// Check whether matrices are equal within a user defined
142    /// precision, set by \a precision.
143    ///
144    /// @return True if each element deviates less or equal than \a
145    /// d. If any matrix contain a NaN, false is always returned.
146    ///
147    bool equal(const matrix&, const double precision=0) const;
148
149    ///
150    /// @return A const pointer to the internal GSL matrix.
151    ///
152    // Jari, is this needed?
153    inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }
154
155    ///
156    /// @return A pointer to the internal GSL matrix.
157    ///
158    // Jari, is this needed?
159    inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };
160
161    ///
162    /// @return True if all elements in the matrix is zero, false
163    /// othwerwise;
164    ///
165    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
166
167    ///
168    /// Check if the matrix object is a view (sub-matrix) to another
169    /// matrix.
170    ///
171    /// @return True if the object is a view, false othwerwise.
172    ///
173    inline bool isview(void) const { return view_; }
174
175    ///
176    /// @return The maximum value of the matrix.
177    ///
178    inline double max(void) const { return gsl_matrix_max(m_); }
179
180    ///
181    /// @return The index to the element with the maximum value of the
182    /// matrix. The lowest index has precedence (searching in
183    /// row-major order).
184    ///
185    inline std::pair<size_t,size_t> max_index(void) const;
186
187    ///
188    /// @return The minimum value of the matrix.
189    ///
190    inline double min(void) const { return gsl_matrix_min(m_); }
191
192    ///
193    /// @return The index to the element with the minimum value of the
194    /// matrix. The lowest index has precedence (searching in
195    /// row-major order).
196    ///
197    inline std::pair<size_t,size_t> min_index(void) const;
198
199    ///
200    /// @return The indecies to the element with the minimum and
201    /// maximum values of the matrix, respectively. The lowest index
202    /// has precedence (searching in row-major order).
203    ///
204    inline void minmax_index(std::pair<size_t,size_t>& min,
205                             std::pair<size_t,size_t>& max) const
206      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
207                                &max.first,&max.second); }
208
209    ///
210    /// @return The minimum and maximum values of the matrix, as the
211    /// \a first and \a second member of the returned \a pair,
212    /// respectively.
213    ///
214    std::pair<double,double> minmax(void) const;
215
216    ///
217    /// Multiply the elements of matrix \a b with the elements of the
218    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
219    /// i,j\f$. The result is stored into the calling matrix.
220    ///
221    /// @return Whatever GSL returns.
222    ///
223    inline int
224    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
225
226    ///
227    /// @return The number of rows in the matrix.
228    ///
229    inline size_t rows(void) const { return m_->size1; }
230
231    ///
232    /// Multiply the elements of the calling matrix with a scalar \a
233    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
234    /// stored into the calling matrix.
235    ///
236    /// @return Whatever GSL returns.
237    ///
238    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
239
240    ///
241    /// Set all elements to \a value.
242    ///
243    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
244
245    ///
246    /// Set \a column to \a vec.
247    ///
248    /// @return Whatever GSL returns.
249    ///
250    inline int set_column(const size_t column, const vector& vec)
251      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_pointer()); }
252
253    ///
254    /// Set \a row to \a vec.
255    ///
256    /// @return Whatever GSL returns.
257    ///
258    inline int set_row(const size_t row, const vector& vec)
259      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_pointer()); }
260
261    ///
262    /// Subtract the elements of matrix \a b from the elements of the
263    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
264    /// i,j\f$. The result is stored into the calling matrix.
265    ///
266    /// @return Whatever GSL returns.
267    ///
268    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
269
270    ///
271    /// Exchange the elements of the this and \a other by copying. The
272    /// two matrices must have the same size.
273    ///
274    /// @return Whatever GSL returns.
275    ///
276    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
277
278    ///
279    /// Swap columns \a i and \a j.
280    ///
281    inline int swap_columns(const size_t i,const size_t j)
282      { return gsl_matrix_swap_columns(m_,i,j); }
283
284    ///
285    /// Swap row \a i and column \a j.
286    ///
287    inline int swap_rowcol(const size_t i,const size_t j)
288      { return gsl_matrix_swap_rowcol(m_,i,j); }
289
290    ///
291    /// Swap rows \a i and \a j.
292    ///
293    inline int swap_rows(const size_t i, const size_t j)
294      { return gsl_matrix_swap_rows(m_,i,j); }
295
296    ///
297    /// Transpose the matrix.
298    ///
299    void transpose(void);
300
301    ///
302    /// @return Reference to the element position (\a row, \a column).
303    ///
304    inline double& operator()(size_t row,size_t column)
305      { return (*gsl_matrix_ptr(m_,row,column)); }
306
307    ///
308    /// @return Const reference to the element position (\a row, \a
309    /// column).
310    ///
311    inline const double& operator()(size_t row,size_t column) const
312      { return (*gsl_matrix_const_ptr(m_,row,column)); }
313
314    ///
315    /// Matrix-vector multiplication.
316    ///
317    /// @return The resulting vector.
318    ///
319    // Jari, where should this go?
320    const vector operator*(const vector&) const;
321
322    ///
323    /// Comparison operator.
324    ///
325    /// @return True if all elements are equal otherwise False.
326    ///
327    /// @see equal
328    ///
329    inline bool operator==(const matrix& other) const { return equal(other); }
330
331    ///
332    /// Comparison operator.
333    ///
334    /// @return False if all elements are equal otherwise True.
335    ///
336    /// @see equal
337    ///
338    inline bool operator!=(const matrix& other) const { return !equal(other); }
339
340    ///
341    /// The assignment operator. There is no requirements on
342    /// dimensions.
343    ///
344    const matrix& operator=(const matrix& other);
345
346    ///
347    /// Add and assign operator.
348    ///
349    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
350
351    ///
352    /// Add and assign operator.
353    ///
354    inline const matrix&
355    operator+=(const double d) { add_constant(d); return *this; }
356
357    ///
358    /// Subtract and assign operator.
359    ///
360    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
361
362    ///
363    /// Multiply and assigment operator.
364    ///
365    /// @return Const reference to the resulting matrix.
366    ///
367    const matrix& operator*=(const matrix&);
368
369    ///
370    /// Multiply and assign operator.
371    ///
372    inline const matrix& operator*=(const double d) { scale(d); return *this; }
373
374
375  private:
376
377    ///
378    /// Create a new copy of the internal GSL matrix.
379    ///
380    /// Necessary memory for the new GSL matrix is allocated and the
381    /// caller is responsible for freeing the allocated memory.
382    ///
383    /// @return A pointer to a copy of the internal GSL matrix.
384    ///
385    gsl_matrix* gsl_matrix_copy(void) const;
386
387    gsl_matrix* m_;
388    gsl_matrix_view* view_;
389  };
390
391  ///
392  /// The output operator for the matrix class.
393  ///
394  std::ostream& operator<< (std::ostream& s, const matrix&);
395
396
397}} // of namespace gslapi and namespace theplu
398
399#endif
Note: See TracBrowser for help on using the repository browser.