source: trunk/yat/utility/Matrix.h @ 1487

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

Addresses #436. GPL license copy reference should also be updated.

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