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

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

Fixes #194.

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