yat/utility/MatrixMutable.h

Code
Comments
Other
Rev Date Author Line
4119 07 Nov 21 peter 1 #ifndef _theplu_yat_utility_matrix_mutable_
4119 07 Nov 21 peter 2 #define _theplu_yat_utility_matrix_mutable_
616 31 Aug 06 jari 3
42 26 Feb 04 jari 4 // $Id$
12 19 Jun 03 daniel 5
570 05 Apr 06 jari 6 /*
570 05 Apr 06 jari 7   Copyright (C) 2003 Daniel Dalevi, Peter Johansson
2119 12 Dec 09 peter 8   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 9   Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 10   Copyright (C) 2007, 2008, 2009 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 11   Copyright (C) 2012, 2017, 2018, 2020, 2021, 2022 Peter Johansson
570 05 Apr 06 jari 12
4119 07 Nov 21 peter 13   This file is part of the yat library, https://dev.thep.lu.se/yat
570 05 Apr 06 jari 14
675 10 Oct 06 jari 15   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 16   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 17   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 18   License, or (at your option) any later version.
570 05 Apr 06 jari 19
675 10 Oct 06 jari 20   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 21   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 22   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
570 05 Apr 06 jari 23   General Public License for more details.
570 05 Apr 06 jari 24
570 05 Apr 06 jari 25   You should have received a copy of the GNU General Public License
4119 07 Nov 21 peter 26   along with yat. If not, see <https://www.gnu.org/licenses/>.
570 05 Apr 06 jari 27 */
570 05 Apr 06 jari 28
1064 10 Feb 08 peter 29 #include "Exception.h"
4119 07 Nov 21 peter 30 #include "MatrixBase.h"
3605 27 Jan 17 peter 31 #include "MatrixExpression.h"
1064 10 Feb 08 peter 32 #include "StrideIterator.h"
1120 21 Feb 08 peter 33 #include "Vector.h"
3605 27 Jan 17 peter 34 #include "VectorExpression.h"
1015 01 Feb 08 peter 35 #include "VectorView.h"
3605 27 Jan 17 peter 36 #include "yat_assert.h"
12 19 Jun 03 daniel 37
12 19 Jun 03 daniel 38 #include <gsl/gsl_matrix.h>
12 19 Jun 03 daniel 39
1792 11 Feb 09 peter 40 #include <cstddef> // size_t
1792 11 Feb 09 peter 41 #include <iosfwd>
1792 11 Feb 09 peter 42
42 26 Feb 04 jari 43 namespace theplu {
680 11 Oct 06 jari 44 namespace yat {
616 31 Aug 06 jari 45 namespace utility {
12 19 Jun 03 daniel 46
1017 01 Feb 08 peter 47   class VectorBase;
43 26 Feb 04 jari 48
4129 19 Jan 22 peter 49   namespace detail {
4129 19 Jan 22 peter 50     class Mover;
4129 19 Jan 22 peter 51   }
4129 19 Jan 22 peter 52
4119 07 Nov 21 peter 53   /**
4119 07 Nov 21 peter 54      Interface class for mutable matrices including Matrix and
4119 07 Nov 21 peter 55      view. The interface allows altering elements of the matrix data
4119 07 Nov 21 peter 56      but not alter the dimensions.
4119 07 Nov 21 peter 57    */
4119 07 Nov 21 peter 58   class MatrixMutable : public MatrixBase
420 02 Dec 05 jari 59   {
420 02 Dec 05 jari 60   public:
1064 10 Feb 08 peter 61     /**
1547 03 Oct 08 peter 62        reference type is double&
1547 03 Oct 08 peter 63
1547 03 Oct 08 peter 64        \since New in yat 0.5
1547 03 Oct 08 peter 65      */
1547 03 Oct 08 peter 66     typedef double& reference;
1547 03 Oct 08 peter 67
1547 03 Oct 08 peter 68     /**
1125 22 Feb 08 peter 69        Mutable iterator that iterates over all elements
1064 10 Feb 08 peter 70      */
4136 21 Jan 22 peter 71     typedef Container2DIterator<MatrixMutable, double, double&> iterator;
12 19 Jun 03 daniel 72
614 31 Aug 06 jari 73     /**
1125 22 Feb 08 peter 74        Mutable iterator that iterates over one column
1103 18 Feb 08 peter 75      */
1103 18 Feb 08 peter 76     typedef StrideIterator<double*> column_iterator;
1103 18 Feb 08 peter 77
1103 18 Feb 08 peter 78     /**
1125 22 Feb 08 peter 79        Mutable iterator that iterates over one row
1103 18 Feb 08 peter 80      */
1103 18 Feb 08 peter 81     typedef StrideIterator<double*> row_iterator;
1103 18 Feb 08 peter 82
1103 18 Feb 08 peter 83     /**
4124 13 Jan 22 peter 84      */
4119 07 Nov 21 peter 85     MatrixMutable(void);
12 19 Jun 03 daniel 86
754 17 Feb 07 jari 87     /**
754 17 Feb 07 jari 88        \brief The copy constructor.
754 17 Feb 07 jari 89     */
4119 07 Nov 21 peter 90     MatrixMutable(const MatrixMutable&) = delete;
12 19 Jun 03 daniel 91
3605 27 Jan 17 peter 92     /**
754 17 Feb 07 jari 93        \brief The destructor.
754 17 Feb 07 jari 94     */
4119 07 Nov 21 peter 95     virtual ~MatrixMutable(void);
17 08 Jul 03 peter 96
1009 01 Feb 08 peter 97     ///
1009 01 Feb 08 peter 98     /// Set all elements to \a value.
1009 01 Feb 08 peter 99     ///
1009 01 Feb 08 peter 100     void all(const double value);
1009 01 Feb 08 peter 101
4119 07 Nov 21 peter 102     using MatrixBase::begin;
4119 07 Nov 21 peter 103     using MatrixBase::begin_column;
4119 07 Nov 21 peter 104     using MatrixBase::begin_row;
4119 07 Nov 21 peter 105
792 11 Mar 07 jari 106     /**
1065 10 Feb 08 peter 107        Iterator iterates along a row. When end of row is reached it
1065 10 Feb 08 peter 108        jumps to beginning of next row.
1065 10 Feb 08 peter 109
1065 10 Feb 08 peter 110        \return iterator pointing to upper-left element.
1064 10 Feb 08 peter 111      */
1064 10 Feb 08 peter 112     iterator begin(void);
1064 10 Feb 08 peter 113
1064 10 Feb 08 peter 114     /**
1065 10 Feb 08 peter 115        Iterator iterates along a column.
1065 10 Feb 08 peter 116
1065 10 Feb 08 peter 117        \return iterator pointing to first element of column \a i.
1064 10 Feb 08 peter 118      */
4133 21 Jan 22 peter 119     column_iterator begin_column(size_t i);
1064 10 Feb 08 peter 120
1064 10 Feb 08 peter 121     /**
1065 10 Feb 08 peter 122        Iterator iterates along a row.
1065 10 Feb 08 peter 123
1065 10 Feb 08 peter 124        \return iterator pointing to first element of row \a i.
1064 10 Feb 08 peter 125      */
4133 21 Jan 22 peter 126     row_iterator begin_row(size_t i);
1064 10 Feb 08 peter 127
1064 10 Feb 08 peter 128     /**
1065 10 Feb 08 peter 129        \return Vector view of column \a i
1009 01 Feb 08 peter 130      */
1065 10 Feb 08 peter 131     VectorView column_view(size_t i);
1009 01 Feb 08 peter 132
1009 01 Feb 08 peter 133     /**
1822 24 Feb 09 jari 134        Elementwise division of the elements of the calling matrix by
755 18 Feb 07 jari 135        the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
755 18 Feb 07 jari 136        \forall i,j \f$. The result is stored into the calling matrix.
420 02 Dec 05 jari 137
1822 24 Feb 09 jari 138        \throw GSL_error if dimensions mismatch.
755 18 Feb 07 jari 139     */
4119 07 Nov 21 peter 140     void div(const MatrixBase& b);
755 18 Feb 07 jari 141
4119 07 Nov 21 peter 142     using MatrixBase::end;
4119 07 Nov 21 peter 143     using MatrixBase::end_column;
4119 07 Nov 21 peter 144     using MatrixBase::end_row;
4119 07 Nov 21 peter 145
788 10 Mar 07 jari 146     /**
1065 10 Feb 08 peter 147        \return iterator pointing to end of matrix
1064 10 Feb 08 peter 148      */
1064 10 Feb 08 peter 149     iterator end(void);
1064 10 Feb 08 peter 150
1064 10 Feb 08 peter 151     /**
1065 10 Feb 08 peter 152        \return iterator pointing to end of column \a i
1064 10 Feb 08 peter 153      */
4133 21 Jan 22 peter 154     column_iterator end_column(size_t i);
1064 10 Feb 08 peter 155
1064 10 Feb 08 peter 156     /**
1065 10 Feb 08 peter 157        \return iterator pointing to end of row \a i
1064 10 Feb 08 peter 158      */
4133 21 Jan 22 peter 159     row_iterator end_row(size_t i);
1064 10 Feb 08 peter 160
4119 07 Nov 21 peter 161     using MatrixBase::gsl_matrix_p;
1064 10 Feb 08 peter 162
42 26 Feb 04 jari 163     ///
42 26 Feb 04 jari 164     /// @return A pointer to the internal GSL matrix.
42 26 Feb 04 jari 165     ///
4129 19 Jan 22 peter 166     virtual gsl_matrix* gsl_matrix_p(void)=0;
12 19 Jun 03 daniel 167
755 18 Feb 07 jari 168     /**
4124 13 Jan 22 peter 169        @return A const pointer to the internal GSL matrix.
4124 13 Jan 22 peter 170     */
4129 19 Jan 22 peter 171     virtual const gsl_matrix* gsl_matrix_p(void) const=0;
4124 13 Jan 22 peter 172
4124 13 Jan 22 peter 173     /**
1121 22 Feb 08 peter 174        Multiply the elements of Matrix \a b with the elements of the
1121 22 Feb 08 peter 175        calling Matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
1121 22 Feb 08 peter 176        \f$. The result is stored into the calling Matrix.
420 02 Dec 05 jari 177
1822 24 Feb 09 jari 178        \throw GSL_error if dimensions mismatch.
755 18 Feb 07 jari 179     */
4119 07 Nov 21 peter 180     void mul(const MatrixBase& b);
755 18 Feb 07 jari 181
811 16 Mar 07 peter 182     /**
1065 10 Feb 08 peter 183        \return Vector view of \a i
1009 01 Feb 08 peter 184      */
1028 03 Feb 08 peter 185     VectorView row_view(size_t);
131 09 Aug 04 peter 186
754 17 Feb 07 jari 187     /**
755 18 Feb 07 jari 188        \brief Swap columns \a i and \a j.
12 19 Jun 03 daniel 189
755 18 Feb 07 jari 190        \throw GSL_error if either index is out of bounds.
755 18 Feb 07 jari 191     */
755 18 Feb 07 jari 192     void swap_columns(const size_t i,const size_t j);
755 18 Feb 07 jari 193
754 17 Feb 07 jari 194     /**
755 18 Feb 07 jari 195        \brief Swap row \a i and column \a j.
755 18 Feb 07 jari 196
1121 22 Feb 08 peter 197        The Matrix must be square.
762 20 Feb 07 jari 198
755 18 Feb 07 jari 199        \throw GSL_error if either index is out of bounds, or if matrix
755 18 Feb 07 jari 200        is not square.
755 18 Feb 07 jari 201     */
755 18 Feb 07 jari 202     void swap_rowcol(const size_t i,const size_t j);
755 18 Feb 07 jari 203
755 18 Feb 07 jari 204     /**
755 18 Feb 07 jari 205        \brief Swap rows \a i and \a j.
755 18 Feb 07 jari 206
755 18 Feb 07 jari 207        \throw GSL_error if either index is out of bounds.
755 18 Feb 07 jari 208     */
755 18 Feb 07 jari 209     void swap_rows(const size_t i, const size_t j);
755 18 Feb 07 jari 210
755 18 Feb 07 jari 211     /**
755 18 Feb 07 jari 212        \brief Element access operator.
755 18 Feb 07 jari 213
755 18 Feb 07 jari 214        \return Reference to the element position (\a row, \a column).
755 18 Feb 07 jari 215
755 18 Feb 07 jari 216        \throw If GSL range checks are enabled in the underlying GSL
755 18 Feb 07 jari 217        library a GSL_error exception is thrown if either index is out
755 18 Feb 07 jari 218        of range.
755 18 Feb 07 jari 219     */
716 25 Dec 06 jari 220     double& operator()(size_t row,size_t column);
390 13 Sep 05 peter 221
4119 07 Nov 21 peter 222     // to allow overload from base class
4119 07 Nov 21 peter 223     using MatrixBase::operator();
755 18 Feb 07 jari 224
788 10 Mar 07 jari 225     /**
754 17 Feb 07 jari 226        \brief The assignment operator.
754 17 Feb 07 jari 227
2689 27 Feb 12 peter 228        \note If lhs needs to be resized, views and iterators are
2689 27 Feb 12 peter 229        invalidated.
2689 27 Feb 12 peter 230
1121 22 Feb 08 peter 231        \return A const reference to the resulting Matrix.
754 17 Feb 07 jari 232     */
4119 07 Nov 21 peter 233     MatrixMutable& operator=(const MatrixMutable& other);
12 19 Jun 03 daniel 234
3605 27 Jan 17 peter 235     /**
4119 07 Nov 21 peter 236        \brief The assignment operator.
4119 07 Nov 21 peter 237      */
4119 07 Nov 21 peter 238     MatrixMutable& operator=(const MatrixBase& other);
4119 07 Nov 21 peter 239
4119 07 Nov 21 peter 240     /**
3605 27 Jan 17 peter 241        Assignment from a matrix expression. A matrix expression is the
3605 27 Jan 17 peter 242        result of \c operator+, \c operator-, and operator*, or
3605 27 Jan 17 peter 243        combinations of them.
3605 27 Jan 17 peter 244
3605 27 Jan 17 peter 245        A typical usage looks like
3605 27 Jan 17 peter 246
3605 27 Jan 17 peter 247        \code
3605 27 Jan 17 peter 248        A = B * C + D;
3605 27 Jan 17 peter 249        \endcode
3605 27 Jan 17 peter 250
3605 27 Jan 17 peter 251        where A, B, C, and D are all instances of class Matrix.
3605 27 Jan 17 peter 252
4118 24 Oct 21 peter 253        Typically MatrixExpression only exists as rvalue, and this
4118 24 Oct 21 peter 254        operator is not called but the rvalue version,
4118 24 Oct 21 peter 255        operator=(MatrixExpression<T>&&)
3605 27 Jan 17 peter 256
4118 24 Oct 21 peter 257        Invalidates, references, iterators and views.
4118 24 Oct 21 peter 258
3605 27 Jan 17 peter 259        \since new in yat 0.15
3605 27 Jan 17 peter 260      */
3605 27 Jan 17 peter 261     template<class T>
4124 13 Jan 22 peter 262     MatrixMutable& operator=(const MatrixExpression<T>& rhs)
4124 13 Jan 22 peter 263     {
4124 13 Jan 22 peter 264       rhs.move(blas_result_);
4129 19 Jan 22 peter 265       move_assign(std::move(blas_result_));
4124 13 Jan 22 peter 266       return *this;
4124 13 Jan 22 peter 267     }
3605 27 Jan 17 peter 268
3938 16 Jul 20 peter 269
754 17 Feb 07 jari 270     /**
3585 19 Jan 17 peter 271        \brief Move assignment operator
4119 07 Nov 21 peter 272      */
4119 07 Nov 21 peter 273     MatrixMutable& operator=(MatrixMutable&& other);
3596 21 Jan 17 peter 274
4119 07 Nov 21 peter 275     /**
3605 27 Jan 17 peter 276        Move assignment from a matrix expression. A matrix expression
3605 27 Jan 17 peter 277        is the result of \c operator+, \c operator-, and operator*, or
3605 27 Jan 17 peter 278        combinations of them.
3605 27 Jan 17 peter 279
3605 27 Jan 17 peter 280        A typical usage looks like
3605 27 Jan 17 peter 281
3605 27 Jan 17 peter 282        \code
3605 27 Jan 17 peter 283        A = B * C + D;
3605 27 Jan 17 peter 284        \endcode
3605 27 Jan 17 peter 285
3605 27 Jan 17 peter 286        where B, C, and D are all instances of class Matrix.
3605 27 Jan 17 peter 287
4118 24 Oct 21 peter 288        Invalidates references, iterators, and views.
3605 27 Jan 17 peter 289
3605 27 Jan 17 peter 290        \since new in yat 0.15
3605 27 Jan 17 peter 291      */
3605 27 Jan 17 peter 292     template<class T>
4119 07 Nov 21 peter 293     MatrixMutable& operator=(MatrixExpression<T>&& rhs)
3605 27 Jan 17 peter 294     {
3653 28 Jun 17 peter 295       // This function steals data from rhs into this, and give the old
3605 27 Jan 17 peter 296       // m_ in return to rhs so destructor of rhs takes care of
3605 27 Jan 17 peter 297       // deallocation.
3605 27 Jan 17 peter 298
4124 13 Jan 22 peter 299       rhs.move(blas_result_);
4129 19 Jan 22 peter 300       move_assign(std::move(blas_result_));
3605 27 Jan 17 peter 301       return *this;
3605 27 Jan 17 peter 302     }
3605 27 Jan 17 peter 303
3585 19 Jan 17 peter 304     /**
754 17 Feb 07 jari 305        \brief Add and assign operator.
754 17 Feb 07 jari 306
1121 22 Feb 08 peter 307        Elementwise addition of the elements of Matrix \a b to the left
1121 22 Feb 08 peter 308        hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
773 01 Mar 07 jari 309        \f$.
773 01 Mar 07 jari 310
1121 22 Feb 08 peter 311        \return A const reference to the resulting Matrix.
754 17 Feb 07 jari 312
1822 24 Feb 09 jari 313        \throw GSL_error if dimensions mismatch.
754 17 Feb 07 jari 314     */
4119 07 Nov 21 peter 315     const MatrixMutable& operator+=(const MatrixBase& b);
420 02 Dec 05 jari 316
773 01 Mar 07 jari 317     /**
3609 27 Jan 17 peter 318        \brief Addition and assign operator.
3609 27 Jan 17 peter 319
3903 07 May 20 peter 320        \return A reference to resulting Matrix
3609 27 Jan 17 peter 321     */
3609 27 Jan 17 peter 322     template<class T>
4119 07 Nov 21 peter 323     MatrixMutable& operator+=(const MatrixExpression<T>& rhs)
3609 27 Jan 17 peter 324     {
3609 27 Jan 17 peter 325       *this = *this + rhs;
3609 27 Jan 17 peter 326       return *this;
3609 27 Jan 17 peter 327     }
3609 27 Jan 17 peter 328
3609 27 Jan 17 peter 329     /**
773 01 Mar 07 jari 330        \brief Add and assign operator
420 02 Dec 05 jari 331
1121 22 Feb 08 peter 332        Add the scalar value \a d to the left hand side Matrix, \f$
773 01 Mar 07 jari 333        a_{ij} = a_{ij} + d \; \forall i,j \f$.
773 01 Mar 07 jari 334     */
4119 07 Nov 21 peter 335     const MatrixMutable& operator+=(const double d);
773 01 Mar 07 jari 336
754 17 Feb 07 jari 337     /**
754 17 Feb 07 jari 338        \brief Subtract and assign operator.
754 17 Feb 07 jari 339
1121 22 Feb 08 peter 340        Elementwise subtraction of the elements of Matrix \a b to the
1121 22 Feb 08 peter 341        left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
773 01 Mar 07 jari 342        i,j \f$.
773 01 Mar 07 jari 343
1121 22 Feb 08 peter 344        \return A const reference to the resulting Matrix.
754 17 Feb 07 jari 345
1822 24 Feb 09 jari 346        \throw GSL_error if dimensions mismatch.
754 17 Feb 07 jari 347     */
4119 07 Nov 21 peter 348     const MatrixMutable& operator-=(const MatrixBase&);
90 29 May 04 jari 349
754 17 Feb 07 jari 350     /**
3609 27 Jan 17 peter 351        \brief Subtraction and assign operator.
3609 27 Jan 17 peter 352
3609 27 Jan 17 peter 353        \return A reference to the resulting VectorMutable.
3609 27 Jan 17 peter 354     */
3609 27 Jan 17 peter 355     template<class T>
4119 07 Nov 21 peter 356     MatrixMutable& operator-=(const MatrixExpression<T>& rhs)
3609 27 Jan 17 peter 357     {
3609 27 Jan 17 peter 358       *this = *this - rhs;
3609 27 Jan 17 peter 359       return *this;
3609 27 Jan 17 peter 360     }
3609 27 Jan 17 peter 361
3609 27 Jan 17 peter 362     /**
773 01 Mar 07 jari 363        \brief Subtract and assign operator
773 01 Mar 07 jari 364
1121 22 Feb 08 peter 365        Subtract the scalar value \a d to the left hand side Matrix,
773 01 Mar 07 jari 366        \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
773 01 Mar 07 jari 367     */
4119 07 Nov 21 peter 368     const MatrixMutable& operator-=(const double d);
773 01 Mar 07 jari 369
773 01 Mar 07 jari 370     /**
754 17 Feb 07 jari 371        \brief Multiply and assignment operator
754 17 Feb 07 jari 372
1121 22 Feb 08 peter 373        Multiply the elements of the left hand side Matrix with a
773 01 Mar 07 jari 374        scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
773 01 Mar 07 jari 375
754 17 Feb 07 jari 376        \throw GSL_error if memory allocation fails.
754 17 Feb 07 jari 377     */
4119 07 Nov 21 peter 378     const MatrixMutable& operator*=(double d);
12 19 Jun 03 daniel 379
4119 07 Nov 21 peter 380   protected:
754 17 Feb 07 jari 381     /**
4129 19 Jan 22 peter 382        blas_result_ is used to temporarily store result in BLAS calls.
4129 19 Jan 22 peter 383        Memory is not allocated for blas_result_ until it is needed.
4129 19 Jan 22 peter 384     */
4129 19 Jan 22 peter 385     gsl_matrix* blas_result_;
4129 19 Jan 22 peter 386
4129 19 Jan 22 peter 387     /// Behaves like operator=(const MatrixBase&)
4124 13 Jan 22 peter 388     void copy_assign(const MatrixBase& other);
4124 13 Jan 22 peter 389
4129 19 Jan 22 peter 390     /// Behaves like operator=(const MatrixBase&)
4129 19 Jan 22 peter 391     virtual void copy_assign(const gsl_matrix* rhs)=0;
4129 19 Jan 22 peter 392
4124 13 Jan 22 peter 393     /**
4129 19 Jan 22 peter 394        Behaves like operator=(MatrixMutable&&)
4124 13 Jan 22 peter 395      */
4129 19 Jan 22 peter 396     virtual void move_assign(MatrixMutable&& other);
4124 13 Jan 22 peter 397
4124 13 Jan 22 peter 398     /**
4129 19 Jan 22 peter 399        Behaves like operator=(MatrixMutable&&)
4124 13 Jan 22 peter 400      */
4129 19 Jan 22 peter 401     virtual void move_assign(gsl_matrix*&& rhs)=0;
4119 07 Nov 21 peter 402   private:
4129 19 Jan 22 peter 403     friend class detail::Mover;
4129 19 Jan 22 peter 404     virtual void visit(detail::Mover& mover);
42 26 Feb 04 jari 405   };
12 19 Jun 03 daniel 406
774 01 Mar 07 jari 407   /**
774 01 Mar 07 jari 408      \brief Exchange all elements between the matrices by copying.
774 01 Mar 07 jari 409
3735 08 May 18 peter 410      This function swaps element by element and matrices must have the
3735 08 May 18 peter 411      same size. References and iterators are not invalidated.
774 01 Mar 07 jari 412
4106 27 Sep 21 peter 413      If you do not care about iterators and references, std::swap is
4106 27 Sep 21 peter 414      typically faster (only swapping a couple of pointers).
3735 08 May 18 peter 415
1437 25 Aug 08 peter 416      \throw GSL_error if sizes are not equal.
1887 31 Mar 09 peter 417
4129 19 Jan 22 peter 418      \relates MatrixMutable
774 01 Mar 07 jari 419   */
4119 07 Nov 21 peter 420   void swap(MatrixMutable&, MatrixMutable&);
774 01 Mar 07 jari 421
687 16 Oct 06 jari 422 }}} // of namespace utility, yat, and theplu
12 19 Jun 03 daniel 423
420 02 Dec 05 jari 424 #endif