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

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

dox

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.2 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 811 2007-03-16 01:00:31Z 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  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
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.
67  ///
68  class matrix
69  {
70  public:
71
72    /**
73       @brief The default constructor.
74   
75       This contructor does not initialize underlying (essential)
76       structures.
77    */
78    matrix(void);
79
80    /**
81       \brief Constructor allocating memory space for \a r times \a c
82       elements, and sets all elements to \a init_value.
83
84       \throw GSL_error if memory allocation fails.
85    */
86    matrix(const size_t& r, const size_t& c, double init_value=0);
87
88    /**
89       \brief The copy constructor.
90
91       \note If the object to be copied is a matrix view, the values
92       of the view will be copied, i.e. the view is not copied.
93
94       \throw A GSL_error is indirectly thrown if memory allocation
95       fails.
96    */
97    matrix(const matrix&);
98
99    /**
100       \brief 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 column
104       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 use
115       is undefined.
116
117       \throw GSL_error if a view cannot be set up.
118    */
119    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
120           size_t n_column);
121
122    /**
123       \brief The istream constructor.
124
125       Matrix rows are sepearated with the new line character, and
126       column elements in a row must be separated either with white
127       space characters except the new line (\\n) character (default),
128       or by the delimiter \a sep.  Rows, and columns, are read
129       consecutively and only rectangular matrices are
130       supported. Empty lines are ignored. End of input is at end of
131       file marker.
132
133       \throw GSL_error if memory allocation fails, IO_error if
134       unexpected input is found in the input stream.
135    */
136    explicit matrix(std::istream &, char sep='\0') 
137      throw(utility::IO_error, std::exception);
138
139    /**
140       \brief The destructor.
141    */
142    ~matrix(void);
143
144    /**
145       \brief Make a copy of \a other.
146
147       This function will make a deep copy of \a other. Memory is
148       resized and view state is changed to same as in other.
149    */
150    const matrix& clone(const matrix& other);
151
152    ///
153    /// @return The number of columns in the matrix.
154    ///
155    size_t columns(void) const;
156
157    /**
158       Elementwise division of the elemnts of the calling matrix by
159       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
160       \forall i,j \f$. The result is stored into the calling matrix.
161
162       \throw GSL_error if dimensions mis-match.
163    */
164    void div(const matrix& b);
165
166    /**
167       \brief Check whether matrices are equal within a user defined
168       precision, set by \a precision.
169
170       \return True if each element deviates less or equal than \a
171       d. If any matrix contain a NaN, false is always returned.
172
173       \see operator== and operator!=
174    */
175    bool equal(const matrix&, const double precision=0) const;
176
177    ///
178    /// @return A const pointer to the internal GSL matrix.
179    ///
180    const gsl_matrix* gsl_matrix_p(void) const;
181
182    ///
183    /// @return A pointer to the internal GSL matrix.
184    ///
185    gsl_matrix* gsl_matrix_p(void);
186
187    ///
188    /// @brief Check if the matrix object is a view (sub-matrix) to
189    /// another matrix.
190    ///
191    /// @return True if the object is a view, false othwerwise.
192    ///
193    bool isview(void) const;
194
195    /**
196       Multiply the elements of matrix \a b with the elements of the
197       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
198       \f$. The result is stored into the calling matrix.
199
200       \throw GSL_error if dimensions mis-match.
201    */
202    void mul(const matrix& b);
203
204    /**
205       \brief Resize matrix
206       
207       All elements are set to @a init_value.
208
209       \note underlying GSL matrix is destroyed and a view into this
210       matrix becomes invalid.
211    */
212    void resize(size_t, size_t, double init_value=0);
213
214    ///
215    /// @return The number of rows in the matrix.
216    ///
217    size_t rows(void) const;
218
219    /**
220       \brief Set element values to values in \a mat.
221
222       This function is needed for assignment of viewed elements.
223
224       \see const matrix& operator=(const matrix&)
225
226       \throw GSL_error if dimensions mis-match.
227    */
228    void set(const matrix& mat);
229
230    ///
231    /// Set all elements to \a value.
232    ///
233    void set_all(const double value);
234
235    /**
236       \brief Set \a column values to values in \a vec.
237
238       \note No check on size is done.
239
240       \throw GSL_error if index is out of range or mis-match in
241       sizes.
242    */
243    void set_column(const size_t column, const vector& vec);
244
245    /**
246       \brief Set \a row values to values in \a vec.
247
248       \note No check on size is done.
249
250       \throw GSL_error if index is out of range or mis-match in
251       sizes.
252    */
253    void set_row(const size_t row, const vector& vec);
254
255    /**
256       \brief Swap columns \a i and \a j.
257
258       \throw GSL_error if either index is out of bounds.
259    */
260    void swap_columns(const size_t i,const size_t j);
261
262    /**
263       \brief Swap row \a i and column \a j.
264
265       The matrix must be square.
266
267       \throw GSL_error if either index is out of bounds, or if matrix
268       is not square.
269    */
270    void swap_rowcol(const size_t i,const size_t j);
271
272    /**
273       \brief Swap rows \a i and \a j.
274
275       \throw GSL_error if either index is out of bounds.
276    */
277    void swap_rows(const size_t i, const size_t j);
278
279    /**
280       \brief Transpose the matrix.
281
282       \throw GSL_error if memory allocation fails for the new
283       transposed matrix.
284    */
285    void transpose(void);
286
287    /**
288       \brief Element access operator.
289
290       \return Reference to the element position (\a row, \a column).
291
292       \throw If GSL range checks are enabled in the underlying GSL
293       library a GSL_error exception is thrown if either index is out
294       of range.
295    */
296    double& operator()(size_t row,size_t column);
297
298    /**
299       \brief Element access operator.
300
301       \return Const reference to the element position (\a row, \a
302       column).
303
304       \throw If GSL range checks are enabled in the underlying GSL
305       library a GSL_error exception is thrown if either index is out
306       of range.
307    */
308    const double& operator()(size_t row,size_t column) const;
309
310    /**
311       \brief Comparison operator. Takes squared time.
312
313       Checks are performed with exact matching, i.e., rounding off
314       effects may destroy comparison. Use the equal function for
315       comparing elements within a user defined precision.
316
317       \return True if all elements are equal otherwise false.
318
319       \see equal
320    */
321    bool operator==(const matrix& other) const;
322
323    /**
324       \brief Comparison operator. Takes squared time.
325
326       Checks are performed with exact matching, i.e., rounding off
327       effects may destroy comparison. Use the equal function for
328       comparing elements within a user defined precision.
329
330       \return False if all elements are equal otherwise true.
331
332       \see equal
333    */
334    bool operator!=(const matrix& other) const;
335
336    /**
337       \brief The assignment operator.
338
339       Dimensions of the matrices must match. If the LHS matrix is a
340       view, the underlying data will be changed.
341
342       \return A const reference to the resulting matrix.
343
344       \see void set(const matrix&).
345
346       \throw GSL_error if dimensions mis-match.
347    */
348    const matrix& operator=(const matrix& other);
349
350    /**
351       \brief Add and assign operator.
352
353       Elementwise addition of the elements of matrix \a b to the left
354       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
355       \f$.
356
357       \return A const reference to the resulting matrix.
358
359       \throw GSL_error if dimensions mis-match.
360    */
361    const matrix& operator+=(const matrix& b);
362
363    /**
364       \brief Add and assign operator
365
366       Add the scalar value \a d to the left hand side matrix, \f$
367       a_{ij} = a_{ij} + d \; \forall i,j \f$.
368    */
369    const matrix& operator+=(const double d);
370
371    /**
372       \brief Subtract and assign operator.
373
374       Elementwise subtraction of the elements of matrix \a b to the
375       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
376       i,j \f$.
377
378       \return A const reference to the resulting matrix.
379
380       \throw GSL_error if dimensions mis-match.
381    */
382    const matrix& operator-=(const matrix&);
383
384    /**
385       \brief Subtract and assign operator
386
387       Subtract the scalar value \a d to the left hand side matrix,
388       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
389    */
390    const matrix& operator-=(const double d);
391
392    /**
393       \brief Multiply and assigment operator.
394
395       \return Const reference to the resulting matrix.
396
397       \throw GSL_error if memory allocation fails.
398    */
399    const matrix& operator*=(const matrix&);
400
401    /**
402       \brief Multiply and assignment operator
403
404       Multiply the elements of the left hand side matrix with a
405       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
406
407       \throw GSL_error if memory allocation fails.
408    */
409    const matrix& operator*=(double d);
410
411  private:
412
413    /**
414       \brief Create a new copy of the internal GSL matrix.
415
416       Necessary memory for the new GSL matrix is allocated and the
417       caller is responsible for freeing the allocated memory.
418
419       \return A pointer to a copy of the internal GSL matrix.
420
421       \throw GSL_error if memory cannot be allocated for the new
422       copy, or if dimensions mis-match.
423    */
424    gsl_matrix* create_gsl_matrix_copy(void) const;
425
426    /**
427       \brief Clear all dynamically allocated memory.
428
429       Internal utility function.
430    */
431    void delete_allocated_memory(void);
432
433    // blas_result_ is used to temporarily store result in BLAS calls
434    // and when not NULL it should always have the same size as
435    // m_. Memory is not allocated for blas_result_ until it is
436    // needed.
437    gsl_matrix* blas_result_;
438    gsl_matrix* m_;
439    gsl_matrix_view* view_;
440    const gsl_matrix_const_view* view_const_;
441    // proxy_m_ is used to access the proper underlying gsl_matrix in
442    // all const member functions. It is not used by design for
443    // non-const vector functions and operators. This is to make sure
444    // that runtime errors occur if a const matrix is used in an
445    // inappropriate manner such as on left hand side in assignment
446    // (remember, m_ is null for const matrix views).
447    const gsl_matrix* proxy_m_;
448  };
449
450  /**
451     \brief Check if all elements of the matrix are zero.
452
453     \return True if all elements in the matrix is zero, false
454     othwerwise.
455  */
456  bool isnull(const matrix&);
457
458  /**
459     \brief Get the maximum value of the matrix.
460
461     \return The maximum value of the matrix.
462  */
463  double max(const matrix&);
464
465  /**
466     \brief Get the minimum value of the matrix.
467
468     \return The minimum value of the matrix.
469  */
470  double min(const matrix&);
471
472  /**
473     \brief Locate the maximum and minumum element in the matrix.
474
475     \return The indecies to the element with the minimum and maximum
476     values of the matrix, respectively.
477
478     \note Lower index has precedence (searching in row-major order).
479  */
480  void minmax_index(const matrix&,
481                    std::pair<size_t,size_t>& min,
482                    std::pair<size_t,size_t>& max);
483
484  /**
485     \brief Create a matrix \a flag indicating NaN's in another matrix
486     \a templat.
487
488     The \a flag matrix is changed to contain 1's and 0's only. A 1
489     means that the corresponding element in the \a templat matrix is
490     valid and a zero means that the corresponding element is a NaN.
491
492     \note Space for matrix \a flag is reallocated to fit the size of
493     matrix \a templat if sizes mismatch.
494
495     \return True if the \a templat matrix contains at least one NaN.
496  */
497  bool nan(const matrix& templat, matrix& flag);
498
499  /**
500     \brief Exchange all elements between the matrices by copying.
501
502     The two matrices must have the same size.
503
504     \throw GSL_error if either index is out of bounds.
505  */
506  void swap(matrix&, matrix&);
507
508  /**
509     \brief The output operator for the matrix class.
510  */
511  std::ostream& operator<< (std::ostream& s, const matrix&);
512
513  /**
514     \brief vector matrix multiplication
515   */
516  vector operator*(const matrix&, const vector&);
517
518  /**
519     \brief matrix vector multiplication
520   */
521  vector operator*(const vector&, const matrix&);
522
523}}} // of namespace utility, yat, and theplu
524
525#endif
Note: See TracBrowser for help on using the repository browser.