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

Last change on this file since 1015 was 1015, checked in by Peter, 13 years ago

changing class names vectorView => VectorView? and vectorBase => VectorBase?. Refs #256

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