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

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

trac moved to new location.

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