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

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

Fixes #76. Creating/recreating BLAS support matrix only when needed.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: matrix.h 762 2007-02-20 19:18: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  /// This is the yat interface to GSL matrix. 'double' is the
46  /// only type supported, maybe we should add a 'complex' type as
47  /// well in the future.
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       Elementwise addition of the elements of matrix \a b to the
147       elements of the calling matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \;
148       \forall i,j \f$. The result is stored into the calling matrix.
149
150       \throw GSL_error if dimensions mis-match.
151    */
152    void add(const matrix& b);
153
154    /**
155       Add the scalar value \a d to the elements of the calling
156       matrix, \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$. The result
157       is stored into the calling matrix.
158    */
159    void add(double d);
160
161    ///
162    /// @return The number of columns in the matrix.
163    ///
164    size_t columns(void) const;
165
166    /**
167       Elementwise division of the elemnts of the calling matrix by
168       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
169       \forall i,j \f$. The result is stored into the calling matrix.
170
171       \throw GSL_error if dimensions mis-match.
172    */
173    void div_elements(const matrix& b);
174
175    ///
176    /// Check whether matrices are equal within a user defined
177    /// precision, set by \a precision.
178    ///
179    /// @return True if each element deviates less or equal than \a
180    /// d. If any matrix contain a NaN, false is always returned.
181    ///
182    bool equal(const matrix&, const double precision=0) const;
183
184    ///
185    /// @return A const pointer to the internal GSL matrix.
186    ///
187    const gsl_matrix* gsl_matrix_p(void) const;
188
189    ///
190    /// @return A pointer to the internal GSL matrix.
191    ///
192    gsl_matrix* gsl_matrix_p(void);
193
194    ///
195    /// @return True if all elements in the matrix is zero, false
196    /// othwerwise;
197    ///
198    bool isnull(void) const;
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    /// @return The maximum value of the matrix.
210    ///
211    double max(void) const;
212
213    ///
214    /// @return The minimum value of the matrix.
215    ///
216    double min(void) const;
217
218    ///
219    /// @return The indecies to the element with the minimum and
220    /// maximum values of the matrix, respectively. The lowest index
221    /// has precedence (searching in row-major order).
222    ///
223    void minmax_index(std::pair<size_t,size_t>& min,
224                      std::pair<size_t,size_t>& max) const;
225
226    ///
227    /// The argument is changed into a matrix containing 1's and 0's
228    /// only. 1 means element in matrix is valid and a zero means
229    /// element is a NaN.
230    ///
231    /// @return true if matrix contains at least one NaN.
232    ///
233    bool nan(matrix&) const;
234
235    /**
236       Multiply the elements of matrix \a b with the elements of the
237       calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
238       \f$. The result is stored into the calling matrix.
239
240       \throw GSL_error if dimensions mis-match.
241    */
242    void mul_elements(const matrix& b);
243
244    ///
245    /// @return The number of rows in the matrix.
246    ///
247    size_t rows(void) const;
248
249    /**
250       Multiply the elements of the calling matrix with a scalar \a d,
251       \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$. The result is
252       stored into the calling matrix.
253    */
254    void scale(const double d);
255
256    /**
257       \brief Set element values to values in \a mat.
258
259       This function is needed for assignment of viewed elements.
260
261       \see const matrix& operator=(const matrix&)
262
263       \throw GSL_error if dimensions mis-match.
264    */
265    void set(const matrix& mat);
266
267    ///
268    /// Set all elements to \a value.
269    ///
270    void set_all(const double value);
271
272    /**
273       \brief Set \a column values to values in \a vec.
274
275       \note No check on size is done.
276
277       \throw GSL_error if index is out of range or mis-match in
278       sizes.
279    */
280    void set_column(const size_t column, const vector& vec);
281
282    /**
283       \brief Set \a row values to values in \a vec.
284
285       \note No check on size is done.
286
287       \throw GSL_error if index is out of range or mis-match in
288       sizes.
289    */
290    void set_row(const size_t row, const vector& vec);
291
292    /**
293       Subtract the elements of matrix \a b from the elements of the
294       calling matrix ,\f$ a_{ij} = a_{ij} - b_{ij} \; \forall i,j
295       \f$. The result is stored into the calling matrix.
296
297       \throw GSL_error if dimensions mis-match.
298    */
299    void sub(const matrix& b);
300
301    /**
302       \brief Exchange the elements of the this and \a other by copying.
303
304       The two matrices must have the same size.
305
306       \throw GSL_error if either index is out of bounds.
307    */
308    void swap(matrix& other);
309
310    /**
311       \brief Swap columns \a i and \a j.
312
313       \throw GSL_error if either index is out of bounds.
314    */
315    void swap_columns(const size_t i,const size_t j);
316
317    /**
318       \brief Swap row \a i and column \a j.
319
320       The matrix must be square.
321
322       \throw GSL_error if either index is out of bounds, or if matrix
323       is not square.
324    */
325    void swap_rowcol(const size_t i,const size_t j);
326
327    /**
328       \brief Swap rows \a i and \a j.
329
330       \throw GSL_error if either index is out of bounds.
331    */
332    void swap_rows(const size_t i, const size_t j);
333
334    /**
335       \brief Transpose the matrix.
336
337       \throw GSL_error if memory allocation fails for the new
338       transposed matrix.
339    */
340    void transpose(void);
341
342    /**
343       \brief Element access operator.
344
345       \return Reference to the element position (\a row, \a column).
346
347       \throw If GSL range checks are enabled in the underlying GSL
348       library a GSL_error exception is thrown if either index is out
349       of range.
350    */
351    double& operator()(size_t row,size_t column);
352
353    /**
354       \brief Element access operator.
355
356       \return Const reference to the element position (\a row, \a
357       column).
358
359       \throw If GSL range checks are enabled in the underlying GSL
360       library a GSL_error exception is thrown if either index is out
361       of range.
362    */
363    const double& operator()(size_t row,size_t column) const;
364
365    ///
366    /// @brief Comparison operator.
367    ///
368    /// @return True if all elements are equal otherwise False.
369    ///
370    /// @see equal
371    ///
372    bool operator==(const matrix& other) const;
373
374    ///
375    /// @brief Comparison operator.
376    ///
377    /// @return False if all elements are equal otherwise True.
378    ///
379    /// @see equal
380    ///
381    bool operator!=(const matrix& other) const;
382
383    /**
384       \brief The assignment operator.
385
386       There is no requirements on dimensions, i.e. the matrix is
387       remapped in memory if necessary. This implies that in general
388       views cannot be assigned using this operator. Views will be
389       mutated into normal matrices. The only exception to this
390       behaviour on views is when self-assignemnt is done, since
391       self-assignment is ignored.
392
393       \return A const reference to the resulting matrix.
394
395       \see void set(const matrix&).
396
397       \throw A GSL_error is indirectly thrown if memory allocation
398       fails, or if dimensions mis-match.
399    */
400    const matrix& operator=(const matrix& other);
401
402    /**
403       \brief Add and assign operator.
404
405       \return A const reference to the resulting matrix.
406
407       \throw GSL_error if dimensions mis-match.
408    */
409    const matrix& operator+=(const matrix&);
410
411    ///
412    /// @brief Add and assign operator.
413    ///
414    const matrix& operator+=(const double);
415
416    /**
417       \brief Subtract and assign operator.
418
419       \return A const reference to the resulting matrix.
420
421       \throw GSL_error if dimensions mis-match.
422    */
423    const matrix& operator-=(const matrix&);
424
425    /**
426       \brief Multiply and assigment operator.
427
428       \return Const reference to the resulting matrix.
429
430       \throw GSL_error if memory allocation fails.
431    */
432    const matrix& operator*=(const matrix&);
433
434    /**
435       \brief Multiply and assignment operator
436
437       \throw GSL_error if memory allocation fails.
438    */
439    const matrix& operator*=(const double);
440
441  private:
442
443    /**
444       \brief Create a new copy of the internal GSL matrix.
445
446       Necessary memory for the new GSL matrix is allocated and the
447       caller is responsible for freeing the allocated memory.
448
449       \return A pointer to a copy of the internal GSL matrix.
450
451       \throw GSL_error if memory cannot be allocated for the new
452       copy, or if dimensions mis-match.
453    */
454    gsl_matrix* create_gsl_matrix_copy(void) const;
455
456    // blas_result_ is used to temporarily store result in BLAS calls
457    // and when not NULL it should always have the same size as
458    // m_. Memory is not allocated for blas_result_ until it is
459    // needed.
460    gsl_matrix* blas_result_;
461    gsl_matrix* m_;
462    gsl_matrix_view* view_;
463  };
464
465  ///
466  /// The output operator for the matrix class.
467  ///
468  std::ostream& operator<< (std::ostream& s, const matrix&);
469
470  /**
471     @brief vector matrix multiplication
472   */
473  vector operator*(const matrix&, const vector&);
474
475  /**
476     @brief matrix vector multiplication
477   */
478  vector operator*(const vector&, const matrix&);
479
480}}} // of namespace utility, yat, and theplu
481
482#endif
Note: See TracBrowser for help on using the repository browser.