source: trunk/yat/utility/matrix.h @ 735

Last change on this file since 735 was 735, checked in by Peter, 15 years ago

forgot to check in matrix.h

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