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

Last change on this file since 1009 was 1009, checked in by Peter, 14 years ago

merging branch peter-dev into trunk delta 1008:994

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 1009 2008-02-01 04:15:44Z peter $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
10  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
11
12  This file is part of the yat library, http://trac.thep.lu.se/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 "vectorView.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    /// Set all elements to \a value.
146    ///
147    void all(const double value);
148
149    /**
150       \brief Make a copy of \a other.
151
152       This function will make a deep copy of \a other. Memory is
153       resized and view state is changed to same as in other.
154    */
155    const matrix& clone(const matrix& other);
156
157    /**
158     */
159    vectorView column_vec(size_t);
160
161    /**
162     */
163    const vectorView column_vec(size_t) const;
164
165    ///
166    /// @return The number of columns in the matrix.
167    ///
168    size_t columns(void) const;
169
170    /**
171       Elementwise division of the elemnts of the calling matrix by
172       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
173       \forall i,j \f$. The result is stored into the calling matrix.
174
175       \throw GSL_error if dimensions mis-match.
176    */
177    void div(const matrix& b);
178
179    /**
180       \brief Check whether matrices are equal within a user defined
181       precision, set by \a precision.
182
183       \return True if each element deviates less or equal than \a
184       d. If any matrix contain a NaN, false is always returned.
185
186       \see operator== and operator!=
187    */
188    bool equal(const matrix&, const double precision=0) const;
189
190    ///
191    /// @return A const pointer to the internal GSL matrix.
192    ///
193    const gsl_matrix* gsl_matrix_p(void) const;
194
195    ///
196    /// @return A pointer to the internal GSL matrix.
197    ///
198    gsl_matrix* gsl_matrix_p(void);
199
200    ///
201    /// @brief Check if the matrix object is a view (sub-matrix) to
202    /// another matrix.
203    ///
204    /// @return True if the object is a view, false othwerwise.
205    ///
206    bool isview(void) const;
207
208    /**
209       Multiply the elements of matrix \a b with the elements of the
210       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
211       \f$. The result is stored into the calling matrix.
212
213       \throw GSL_error if dimensions mis-match.
214    */
215    void mul(const matrix& b);
216
217    /**
218       \brief Resize matrix
219       
220       All elements are set to @a init_value.
221
222       \note underlying GSL matrix is destroyed and a view into this
223       matrix becomes invalid.
224    */
225    void resize(size_t, size_t, double init_value=0);
226
227    ///
228    /// @return The number of rows in the matrix.
229    ///
230    size_t rows(void) const;
231
232    /**
233     */
234    vectorView row_vec(size_t);
235
236    /**
237     */
238    const vectorView row_vec(size_t) const;
239
240    /**
241       \brief Swap columns \a i and \a j.
242
243       \throw GSL_error if either index is out of bounds.
244    */
245    void swap_columns(const size_t i,const size_t j);
246
247    /**
248       \brief Swap row \a i and column \a j.
249
250       The matrix must be square.
251
252       \throw GSL_error if either index is out of bounds, or if matrix
253       is not square.
254    */
255    void swap_rowcol(const size_t i,const size_t j);
256
257    /**
258       \brief Swap rows \a i and \a j.
259
260       \throw GSL_error if either index is out of bounds.
261    */
262    void swap_rows(const size_t i, const size_t j);
263
264    /**
265       \brief Transpose the matrix.
266
267       \throw GSL_error if memory allocation fails for the new
268       transposed matrix.
269    */
270    void transpose(void);
271
272    /**
273       \brief Element access operator.
274
275       \return Reference to the element position (\a row, \a column).
276
277       \throw If GSL range checks are enabled in the underlying GSL
278       library a GSL_error exception is thrown if either index is out
279       of range.
280    */
281    double& operator()(size_t row,size_t column);
282
283    /**
284       \brief Element access operator.
285
286       \return Const reference to the element position (\a row, \a
287       column).
288
289       \throw If GSL range checks are enabled in the underlying GSL
290       library a GSL_error exception is thrown if either index is out
291       of range.
292    */
293    const double& operator()(size_t row,size_t column) const;
294
295    /**
296       \brief Comparison operator. Takes squared time.
297
298       Checks are performed with exact matching, i.e., rounding off
299       effects may destroy comparison. Use the equal function for
300       comparing elements within a user defined precision.
301
302       \return True if all elements are equal otherwise false.
303
304       \see equal
305    */
306    bool operator==(const matrix& other) const;
307
308    /**
309       \brief Comparison operator. Takes squared time.
310
311       Checks are performed with exact matching, i.e., rounding off
312       effects may destroy comparison. Use the equal function for
313       comparing elements within a user defined precision.
314
315       \return False if all elements are equal otherwise true.
316
317       \see equal
318    */
319    bool operator!=(const matrix& other) const;
320
321    /**
322       \brief The assignment operator.
323
324       Dimensions of the matrices must match. If the LHS matrix is a
325       view, the underlying data will be changed.
326
327       \return A const reference to the resulting matrix.
328
329       \see void set(const matrix&).
330
331       \throw GSL_error if dimensions mis-match.
332    */
333    const matrix& operator=(const matrix& other);
334
335    /**
336       \brief Add and assign operator.
337
338       Elementwise addition of the elements of matrix \a b to the left
339       hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
340       \f$.
341
342       \return A const reference to the resulting matrix.
343
344       \throw GSL_error if dimensions mis-match.
345    */
346    const matrix& operator+=(const matrix& b);
347
348    /**
349       \brief Add and assign operator
350
351       Add the scalar value \a d to the left hand side matrix, \f$
352       a_{ij} = a_{ij} + d \; \forall i,j \f$.
353    */
354    const matrix& operator+=(const double d);
355
356    /**
357       \brief Subtract and assign operator.
358
359       Elementwise subtraction of the elements of matrix \a b to the
360       left hand side matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
361       i,j \f$.
362
363       \return A const reference to the resulting matrix.
364
365       \throw GSL_error if dimensions mis-match.
366    */
367    const matrix& operator-=(const matrix&);
368
369    /**
370       \brief Subtract and assign operator
371
372       Subtract the scalar value \a d to the left hand side matrix,
373       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
374    */
375    const matrix& operator-=(const double d);
376
377    /**
378       \brief Multiply and assigment operator.
379
380       \return Const reference to the resulting matrix.
381
382       \throw GSL_error if memory allocation fails.
383    */
384    const matrix& operator*=(const matrix&);
385
386    /**
387       \brief Multiply and assignment operator
388
389       Multiply the elements of the left hand side matrix with a
390       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
391
392       \throw GSL_error if memory allocation fails.
393    */
394    const matrix& operator*=(double d);
395
396  private:
397
398    /**
399       \brief Create a new copy of the internal GSL matrix.
400
401       Necessary memory for the new GSL matrix is allocated and the
402       caller is responsible for freeing the allocated memory.
403
404       \return A pointer to a copy of the internal GSL matrix.
405
406       \throw GSL_error if memory cannot be allocated for the new
407       copy, or if dimensions mis-match.
408    */
409    gsl_matrix* create_gsl_matrix_copy(void) const;
410
411    /**
412       \brief Clear all dynamically allocated memory.
413
414       Internal utility function.
415    */
416    void delete_allocated_memory(void);
417
418    // blas_result_ is used to temporarily store result in BLAS calls
419    // and when not NULL it should always have the same size as
420    // m_. Memory is not allocated for blas_result_ until it is
421    // needed.
422    gsl_matrix* blas_result_;
423    gsl_matrix* m_;
424    gsl_matrix_view* view_;
425    const gsl_matrix_const_view* view_const_;
426    // proxy_m_ is used to access the proper underlying gsl_matrix in
427    // all const member functions. It is not used by design for
428    // non-const vector functions and operators. This is to make sure
429    // that runtime errors occur if a const matrix is used in an
430    // inappropriate manner such as on left hand side in assignment
431    // (remember, m_ is null for const matrix views).
432    const gsl_matrix* proxy_m_;
433  };
434
435  /**
436     \brief Check if all elements of the matrix are zero.
437
438     \return True if all elements in the matrix is zero, false
439     othwerwise.
440  */
441  bool isnull(const matrix&);
442
443  /**
444     \brief Get the maximum value of the matrix.
445
446     \return The maximum value of the matrix.
447  */
448  double max(const matrix&);
449
450  /**
451     \brief Get the minimum value of the matrix.
452
453     \return The minimum value of the matrix.
454  */
455  double min(const matrix&);
456
457  /**
458     \brief Locate the maximum and minumum element in the matrix.
459
460     \return The indecies to the element with the minimum and maximum
461     values of the matrix, respectively.
462
463     \note Lower index has precedence (searching in row-major order).
464  */
465  void minmax_index(const matrix&,
466                    std::pair<size_t,size_t>& min,
467                    std::pair<size_t,size_t>& max);
468
469  /**
470     \brief Create a matrix \a flag indicating NaN's in another matrix
471     \a templat.
472
473     The \a flag matrix is changed to contain 1's and 0's only. A 1
474     means that the corresponding element in the \a templat matrix is
475     valid and a zero means that the corresponding element is a NaN.
476
477     \note Space for matrix \a flag is reallocated to fit the size of
478     matrix \a templat if sizes mismatch.
479
480     \return True if the \a templat matrix contains at least one NaN.
481  */
482  bool nan(const matrix& templat, matrix& flag);
483
484  /**
485     \brief Exchange all elements between the matrices by copying.
486
487     The two matrices must have the same size.
488
489     \throw GSL_error if either index is out of bounds.
490  */
491  void swap(matrix&, matrix&);
492
493  /**
494     \brief The output operator for the matrix class.
495  */
496  std::ostream& operator<< (std::ostream& s, const matrix&);
497
498  /**
499     \brief vector matrix multiplication
500   */
501  vector operator*(const matrix&, const vector&);
502
503  /**
504     \brief matrix vector multiplication
505   */
506  vector operator*(const vector&, const matrix&);
507
508}}} // of namespace utility, yat, and theplu
509
510#endif
Note: See TracBrowser for help on using the repository browser.