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

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

Iterators for KernelLookup? - refs #267

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