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

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

Fixes #299. Memory leak in matrix was found and removed.

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