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

Last change on this file since 3999 was 3999, checked in by Peter, 20 months ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.2 KB
Line 
1#ifndef _theplu_yat_utility_matrix_
2#define _theplu_yat_utility_matrix_
3
4// $Id: Matrix.h 3999 2020-10-08 23:22: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, Peter Johansson, Markus Ringnér
10  Copyright (C) 2007, 2008, 2009 Jari Häkkinen, Peter Johansson
11  Copyright (C) 2012, 2017, 2018, 2019, 2020 Peter Johansson
12
13  This file is part of the yat library, http://dev.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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
27*/
28
29#include <cassert>
30
31#include "config_public.h"
32
33// include these operations as they are expected when using Matrix
34#include "BLAS_level2.h"
35#include "BLAS_level3.h"
36
37#include "Exception.h"
38#include "MatrixExpression.h"
39#include "StrideIterator.h"
40#include "Vector.h"
41#include "VectorConstView.h"
42#include "VectorExpression.h"
43#include "VectorView.h"
44#include "yat_assert.h"
45
46#include <gsl/gsl_matrix.h>
47
48#include <cstddef> // size_t
49#include <iosfwd>
50
51namespace theplu {
52namespace yat {
53namespace utility {
54
55  class VectorBase;
56
57  ///
58  /// @brief Interface to GSL matrix.
59  ///
60  /// For the time being 'double' is the only type supported.
61  ///
62  /// \par[File streams] Reading and writing vectors to file streams
63  /// are of course supported. These are implemented without using GSL
64  /// functionality, and thus binary read and write to streams are not
65  /// supported.
66  ///
67  /// @note All GSL matrix related functions are not implement but
68  /// most functionality defined for GSL matrices can be achieved with
69  /// this interface class. Most notable GSL functionality not
70  /// supported are no binary file support and views on arrays,
71  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
72  /// superdiagonals.
73  ///
74  class Matrix : public BasicMatrix<Matrix>
75  {
76  public:
77    /**
78       value_type is double
79
80       \since New in yat 0.5
81     */
82    typedef double value_type;
83
84    /**
85       reference type is double&
86
87       \since New in yat 0.5
88     */
89    typedef double& reference;
90
91    /**
92       const_reference type is const double&
93
94       \since New in yat 0.5
95     */
96    typedef const double& const_reference;
97
98    /**
99       Mutable iterator that iterates over all elements
100     */
101    typedef StrideIterator<double*> iterator;
102
103    /**
104       Read-only iterator that iterates over all elements
105     */
106    typedef StrideIterator<const double*> const_iterator;
107
108    /**
109       Mutable iterator that iterates over one column
110     */
111    typedef StrideIterator<double*> column_iterator;
112
113    /**
114       Read-only iterator that iterates over one column
115     */
116    typedef StrideIterator<const double*> const_column_iterator;
117
118    /**
119       Mutable iterator that iterates over one row
120     */
121    typedef StrideIterator<double*> row_iterator;
122
123    /**
124       Read-only iterator that iterates over one row
125     */
126    typedef StrideIterator<const double*> const_row_iterator;
127
128    /**
129       @brief The default constructor.
130
131       This constructor does not initialize underlying (essential)
132       structures.
133    */
134    Matrix(void);
135
136    /**
137       \brief Constructor allocating memory space for \a r times \a c
138       elements, and sets all elements to \a init_value.
139
140       If \a r is zero \a c must be zero and vice versa.
141
142       \throw GSL_error if memory allocation fails.
143    */
144    Matrix(const size_t& r, const size_t& c, double init_value=0);
145
146    /**
147       \brief The copy constructor.
148
149       \throw A GSL_error is indirectly thrown if memory allocation
150       fails.
151    */
152    Matrix(const Matrix&);
153
154    /**
155       Constructor from a matrix expression. A matrix expression is
156       result \c operator+, \c operator-, and operator*, or
157       combinations of them.
158
159       A typical usage looks like
160
161       \code
162       Matrix A( B * C + D)
163       \endcode
164
165       where B, C, and D are all instances of class Matrix.
166
167       Typically this constructor is not used, if rvalues are
168       enabled. \see Matrix(MatrixExpression<T>&&).
169
170       \since new in yat 0.15
171     */
172    template<class T>
173    Matrix(const MatrixExpression<T>& other);
174
175    /**
176       \brief The move constructor.
177
178       \since new in yat 0.15
179    */
180    Matrix(Matrix&&) noexcept;
181
182    /**
183       Move constructor from a matrix expression. A matrix expression is
184       result \c operator+, \c operator-, and operator*, or
185       combinations of them.
186
187       A typical usage looks like
188
189       \code
190       Matrix A( B * C + D)
191       \endcode
192
193       where B, C, and D are all instances of class Matrix.
194
195       \since new in yat 0.15
196    */
197    template<class T>
198    Matrix(MatrixExpression<T>&& other)
199      : blas_result_(NULL), m_(NULL)
200    {
201      other.move(m_);
202    }
203
204
205    /**
206       \brief The istream constructor.
207
208       The std::istream will be interpreted as outlined here:
209
210       Missing values, i.e. empty elements, are treated as NaN values
211       (std::numeric_limits<double>::quiet_NaN() to be specific).
212
213       Matrix rows are separated with the new line character.
214
215       Column element separation has two modes depending on the value
216       of \a sep.
217
218       - If \a sep is the default '\\0' value then column elements are
219       separated with white space characters except the new line
220       character. Multiple sequential white space characters are
221       treated as one separator.
222
223       - Setting \a sep to something else than the default value will
224       change the behaviour to use the \a sep character as the
225       separator between column elements. Multiple sequential \a sep
226       characters will be treated as separating elements with missing
227       values.
228
229       End of input is the end of file marker and this treatment
230       cannot be redefined using the provided API.
231
232       \throw GSL_error if memory allocation fails, IO_error if
233       unexpected input is found in the input stream.
234    */
235    explicit Matrix(std::istream &, char sep='\0');
236
237    /**
238       \brief The destructor.
239    */
240    ~Matrix(void);
241
242    ///
243    /// Set all elements to \a value.
244    ///
245    void all(const double value);
246
247    /**
248       Iterator iterates along a row. When end of row is reached it
249       jumps to beginning of next row.
250
251       \return iterator pointing to upper-left element.
252     */
253    iterator begin(void);
254
255    /**
256       Iterator iterates along a row. When end of row is reached it
257       jumps to beginning of next row.
258
259       \return const_iterator pointing to upper-left element.
260     */
261    const_iterator begin(void) const;
262
263    /**
264       Iterator iterates along a column.
265
266       \return iterator pointing to first element of column \a i.
267     */
268    iterator begin_column(size_t i);
269
270    /**
271       Iterator iterates along a column.
272
273       \return const_iterator pointing to first element of column \a i.
274     */
275    const_iterator begin_column(size_t i) const;
276
277    /**
278       Iterator iterates along a row.
279
280       \return iterator pointing to first element of row \a i.
281     */
282    iterator begin_row(size_t i);
283
284    /**
285       Iterator iterates along a row.
286
287       \return const_iterator pointing to first element of row \a i.
288     */
289    const_iterator begin_row(size_t i) const;
290
291    /**
292       \return Vector view of column \a i
293     */
294    VectorView column_view(size_t i);
295
296    /**
297       \return const vector view of column \a i
298     */
299    const VectorConstView column_const_view(size_t) const;
300
301    ///
302    /// @return The number of columns in the matrix.
303    ///
304    size_t columns(void) const;
305
306    /**
307       Elementwise division of the elements of the calling matrix by
308       the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
309       \forall i,j \f$. The result is stored into the calling matrix.
310
311       \throw GSL_error if dimensions mismatch.
312    */
313    void div(const Matrix& b);
314
315    /**
316       \return iterator pointing to end of matrix
317     */
318    iterator end(void);
319
320    /**
321       \return const_iterator pointing to end of matrix
322     */
323    const_iterator end(void) const;
324
325    /**
326       \return iterator pointing to end of column \a i
327     */
328    iterator end_column(size_t i);
329
330    /**
331       \return const_iterator pointing to end of column \a i
332     */
333    const_iterator end_column(size_t i) const;
334
335    /**
336       \return iterator pointing to end of row \a i
337     */
338    iterator end_row(size_t i);
339
340    /**
341       \return const_iterator pointing to end of row \a i
342     */
343    const_iterator end_row(size_t i) const;
344
345    /**
346       \brief Check whether matrices are equal within a user defined
347       precision, set by \a precision.
348
349       \return True if each element deviates less or equal than \a
350       d. If any matrix contain a NaN, false is always returned.
351
352       \see operator== and operator!=
353    */
354    bool equal(const Matrix&, const double precision=0) const;
355
356    ///
357    /// @return A const pointer to the internal GSL matrix.
358    ///
359    const gsl_matrix* gsl_matrix_p(void) const;
360
361    ///
362    /// @return A pointer to the internal GSL matrix.
363    ///
364    gsl_matrix* gsl_matrix_p(void);
365
366    /**
367       Multiply the elements of Matrix \a b with the elements of the
368       calling Matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
369       \f$. The result is stored into the calling Matrix.
370
371       \throw GSL_error if dimensions mismatch.
372    */
373    void mul(const Matrix& b);
374
375    /**
376       \brief Resize Matrix
377
378       All elements are set to @a init_value.
379
380       If \a r is zero \a c must be zero and vice versa.
381
382       \note underlying GSL matrix is destroyed and views and
383       iterators are invalidated
384    */
385    void resize(size_t r, size_t c, double init_value=0);
386
387    ///
388    /// @return The number of rows in the matrix.
389    ///
390    size_t rows(void) const;
391
392    /**
393       \return Vector view of \a i
394     */
395    VectorView row_view(size_t);
396
397    /**
398       \return const vector view of \a i
399     */
400    const VectorConstView row_const_view(size_t) const;
401
402    /**
403       \brief Swap columns \a i and \a j.
404
405       \throw GSL_error if either index is out of bounds.
406    */
407    void swap_columns(const size_t i,const size_t j);
408
409    /**
410       \brief Swap row \a i and column \a j.
411
412       The Matrix must be square.
413
414       \throw GSL_error if either index is out of bounds, or if matrix
415       is not square.
416    */
417    void swap_rowcol(const size_t i,const size_t j);
418
419    /**
420       \brief Swap rows \a i and \a j.
421
422       \throw GSL_error if either index is out of bounds.
423    */
424    void swap_rows(const size_t i, const size_t j);
425
426    /**
427       \brief Transpose the matrix.
428
429       \note Unless matrix is square, views and iterators are
430       invalidated.
431
432       \throw GSL_error if memory allocation fails for the new
433       transposed matrix.
434    */
435    void transpose(void);
436
437    /**
438       \brief Element access operator.
439
440       \return Reference to the element position (\a row, \a column).
441
442       \throw If GSL range checks are enabled in the underlying GSL
443       library a GSL_error exception is thrown if either index is out
444       of range.
445    */
446    double& operator()(size_t row,size_t column);
447
448    /**
449       \brief Element access operator.
450
451       \return Const reference to the element position (\a row, \a
452       column).
453
454       \throw If GSL range checks are enabled in the underlying GSL
455       library a GSL_error exception is thrown if either index is out
456       of range.
457    */
458    const double& operator()(size_t row,size_t column) const;
459
460    /**
461       \brief Comparison operator. Takes squared time.
462
463       Checks are performed with exact matching, i.e., rounding off
464       effects may destroy comparison. Use the equal function for
465       comparing elements within a user defined precision.
466
467       \return True if all elements are equal otherwise false.
468
469       \see equal
470    */
471    bool operator==(const Matrix& other) const;
472
473    /**
474       \brief Comparison operator. Takes squared time.
475
476       Checks are performed with exact matching, i.e., rounding off
477       effects may destroy comparison. Use the equal function for
478       comparing elements within a user defined precision.
479
480       \return False if all elements are equal otherwise true.
481
482       \see equal
483    */
484    bool operator!=(const Matrix& other) const;
485
486    /**
487       \brief The assignment operator.
488
489       \note If lhs needs to be resized, views and iterators are
490       invalidated.
491
492       \return A const reference to the resulting Matrix.
493    */
494    const Matrix& operator=(const Matrix& other);
495
496    /**
497       Assignment from a matrix expression. A matrix expression is the
498       result of \c operator+, \c operator-, and operator*, or
499       combinations of them.
500
501       A typical usage looks like
502
503       \code
504       A = B * C + D;
505       \endcode
506
507       where A, B, C, and D are all instances of class Matrix.
508
509       Typically this constructor is not used, if rvalues are
510       enabled. \see Matrix(MatrixExpression<T>&&).
511
512       \since new in yat 0.15
513     */
514    template<class T>
515    Matrix& operator=(const MatrixExpression<T>& rhs);
516
517
518    /**
519       \brief Move assignment operator
520
521       \since new in yat 0.15
522     */
523    Matrix& operator=(Matrix&& other);
524
525    /**
526       Move assignment from a matrix expression. A matrix expression
527       is the result of \c operator+, \c operator-, and operator*, or
528       combinations of them.
529
530       A typical usage looks like
531
532       \code
533       A = B * C + D;
534       \endcode
535
536       where B, C, and D are all instances of class Matrix.
537
538       Typically this constructor is not used, if rvalues are
539       enabled. \see Matrix(MatrixExpression<T>&&).
540
541       \since new in yat 0.15
542     */
543    template<class T>
544    Matrix& operator=(MatrixExpression<T>&& rhs)
545    {
546      // This function steals data from rhs into this, and give the old
547      // m_ in return to rhs so destructor of rhs takes care of
548      // deallocation.
549
550      rhs.move(m_);
551      return *this;
552    }
553
554
555    /**
556       \brief Add and assign operator.
557
558       Elementwise addition of the elements of Matrix \a b to the left
559       hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
560       \f$.
561
562       \return A const reference to the resulting Matrix.
563
564       \throw GSL_error if dimensions mismatch.
565    */
566    const Matrix& operator+=(const Matrix& b);
567
568    /**
569       \brief Addition and assign operator.
570
571       \return A reference to resulting Matrix
572    */
573    template<class T>
574    Matrix& operator+=(const MatrixExpression<T>& rhs)
575    {
576      *this = *this + rhs;
577      return *this;
578    }
579
580    /**
581       \brief Add and assign operator
582
583       Add the scalar value \a d to the left hand side Matrix, \f$
584       a_{ij} = a_{ij} + d \; \forall i,j \f$.
585    */
586    const Matrix& operator+=(const double d);
587
588    /**
589       \brief Subtract and assign operator.
590
591       Elementwise subtraction of the elements of Matrix \a b to the
592       left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
593       i,j \f$.
594
595       \return A const reference to the resulting Matrix.
596
597       \throw GSL_error if dimensions mismatch.
598    */
599    const Matrix& operator-=(const Matrix&);
600
601    /**
602       \brief Subtraction and assign operator.
603
604       \return A reference to the resulting VectorMutable.
605    */
606    template<class T>
607    Matrix& operator-=(const MatrixExpression<T>& rhs)
608    {
609      *this = *this - rhs;
610      return *this;
611    }
612
613    /**
614       \brief Subtract and assign operator
615
616       Subtract the scalar value \a d to the left hand side Matrix,
617       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
618    */
619    const Matrix& operator-=(const double d);
620
621    /**
622       \brief Multiply and assignment operator.
623
624       \return Const reference to the resulting Matrix.
625
626       \note invalidates views and iterators
627
628       \throw GSL_error if memory allocation fails.
629    */
630    const Matrix& operator*=(const Matrix&);
631
632    /**
633       \brief Multiply and assignment operator
634
635       Multiply the elements of the left hand side Matrix with a
636       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
637
638       \throw GSL_error if memory allocation fails.
639    */
640    const Matrix& operator*=(double d);
641
642  private:
643
644    /**
645       \brief Create a new copy of the internal GSL matrix.
646
647       Necessary memory for the new GSL matrix is allocated and the
648       caller is responsible for freeing the allocated memory.
649
650       \return A pointer to a copy of the internal GSL matrix.
651
652       \throw GSL_error if memory cannot be allocated for the new
653       copy, or if dimensions mismatch.
654    */
655    gsl_matrix* create_gsl_matrix_copy(void) const;
656
657    /**
658       \brief Clear all dynamically allocated memory.
659
660       Internal utility function.
661    */
662    void delete_allocated_memory(void);
663
664    // blas_result_ is used to temporarily store result in BLAS calls
665    // and when not NULL it should always have the same size as
666    // m_. Memory is not allocated for blas_result_ until it is
667    // needed.
668    gsl_matrix* blas_result_;
669    gsl_matrix* m_;
670  };
671
672  /**
673     Calculate inverse of a (square) Matrix using SVD
674
675     \since New in yat 0.15
676
677     \relates Matrix
678   */
679  void inverse_svd(const Matrix& input, Matrix& result);
680
681  /**
682     \brief Check if all elements of the Matrix are zero.
683
684     \return True if all elements in the Matrix is zero, false
685     otherwise.
686
687     \relates Matrix
688  */
689  bool isnull(const Matrix&);
690
691  /**
692     \brief Get the maximum value of the Matrix.
693
694     \return The maximum value of the Matrix.
695
696     \relates Matrix
697  */
698  double max(const Matrix&);
699
700  /**
701     \brief Get the minimum value of the Matrix.
702
703     \return The minimum value of the Matrix.
704
705     \relates Matrix
706  */
707  double min(const Matrix&);
708
709  /**
710     \brief Locate the maximum and minimum element in the Matrix.
711
712     \return The indices to the element with the minimum and maximum
713     values of the Matrix, respectively.
714
715     \note Lower index has precedence (searching in row-major order).
716
717     \relates Matrix
718  */
719  void minmax_index(const Matrix&,
720                    std::pair<size_t,size_t>& min,
721                    std::pair<size_t,size_t>& max);
722
723  /**
724     \brief Create a Matrix \a flag indicating NaN's in another Matrix
725     \a templat.
726
727     The \a flag Matrix is changed to contain 1's and 0's only. A 1
728     means that the corresponding element in the \a templat Matrix is
729     valid and a zero means that the corresponding element is a NaN.
730
731     \note Space for Matrix \a flag is reallocated to fit the size of
732     Matrix \a templat if sizes mismatch.
733
734     \return True if the \a templat Matrix contains at least one NaN.
735
736     \relates Matrix
737  */
738  bool nan(const Matrix& templat, Matrix& flag);
739
740  /**
741     \brief Exchange all elements between the matrices by copying.
742
743     This function swaps element by element and matrices must have the
744     same size. References and iterators are not invalidated.
745
746     If you use c++11, do not care about iterators and references,
747     std::swap is typically faster (only swapping a couple of pointers)
748
749     \throw GSL_error if sizes are not equal.
750
751     \relates Matrix
752  */
753  void swap(Matrix&, Matrix&);
754
755  /**
756     \brief The output operator for the Matrix class.
757
758     \relates Matrix
759  */
760  std::ostream& operator<< (std::ostream& s, const Matrix&);
761
762  /**
763     \brief Trace of matrix
764
765     \relates Matrix
766
767     \return sum of diagonal elements
768
769     \since New in yat 0.9
770  */
771  double trace(const Matrix&);
772
773  /////  template implemantions
774  template<class T>
775  Matrix::Matrix(const MatrixExpression<T>& other)
776    : blas_result_(NULL), m_(NULL)
777  {
778    // In principle, we could implement this as a move constructor
779    // and steal the data (rather than copy), but it's a bit hackish
780    // to workaround the constness, so users who are eager for speed
781    // should enable rvalues and get access to the faster move
782    // constructor below.
783
784    try {
785      detail::copy(m_, other.gsl_matrix_p());
786    }
787    catch (GSL_error& e) {
788      detail::deallocate(m_);
789      throw e;
790    }
791  }
792
793
794  template<class T>
795  Matrix& Matrix::operator=(const MatrixExpression<T>& rhs)
796  {
797    // Be careful because *this might be embedded in rhs.
798    rhs.get(blas_result_);
799    std::swap(m_, blas_result_);
800    return *this;
801  }
802
803
804}}} // of namespace utility, yat, and theplu
805
806#endif
Note: See TracBrowser for help on using the repository browser.