yat/utility/Matrix.h

Code
Comments
Other
Rev Date Author Line
680 11 Oct 06 jari 1 #ifndef _theplu_yat_utility_matrix_
680 11 Oct 06 jari 2 #define _theplu_yat_utility_matrix_
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
1437 25 Aug 08 peter 13   This file is part of the yat library, http://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
1487 10 Sep 08 jari 26   along with yat. If not, see <http://www.gnu.org/licenses/>.
570 05 Apr 06 jari 27 */
570 05 Apr 06 jari 28
4129 19 Jan 22 peter 29 #include <cassert>
4129 19 Jan 22 peter 30
3585 19 Jan 17 peter 31 #include "config_public.h"
3585 19 Jan 17 peter 32
1064 10 Feb 08 peter 33 #include "Exception.h"
3605 27 Jan 17 peter 34 #include "MatrixExpression.h"
4119 07 Nov 21 peter 35 #include "MatrixMutable.h"
1064 10 Feb 08 peter 36 #include "StrideIterator.h"
1120 21 Feb 08 peter 37 #include "Vector.h"
1027 02 Feb 08 peter 38 #include "VectorConstView.h"
3605 27 Jan 17 peter 39 #include "VectorExpression.h"
1015 01 Feb 08 peter 40 #include "VectorView.h"
3605 27 Jan 17 peter 41 #include "yat_assert.h"
12 19 Jun 03 daniel 42
12 19 Jun 03 daniel 43 #include <gsl/gsl_matrix.h>
12 19 Jun 03 daniel 44
1792 11 Feb 09 peter 45 #include <cstddef> // size_t
1792 11 Feb 09 peter 46 #include <iosfwd>
1792 11 Feb 09 peter 47
42 26 Feb 04 jari 48 namespace theplu {
680 11 Oct 06 jari 49 namespace yat {
616 31 Aug 06 jari 50 namespace utility {
12 19 Jun 03 daniel 51
1017 01 Feb 08 peter 52   class VectorBase;
4129 19 Jan 22 peter 53   class Matrix;
43 26 Feb 04 jari 54
4129 19 Jan 22 peter 55   /// \cond IGNORE_DOXYGEN
4129 19 Jan 22 peter 56
4129 19 Jan 22 peter 57   namespace detail {
4129 19 Jan 22 peter 58
4129 19 Jan 22 peter 59     /*
4129 19 Jan 22 peter 60       Class involved in move assignment of MatrixMutable to determine
4129 19 Jan 22 peter 61       when both lhs and rhs are a Matrix. When the behavoiur only
4129 19 Jan 22 peter 62       depends on the type of one object, dynamic polymorphism with a
4129 19 Jan 22 peter 63       virtual function suffices; when the behaviour depends on two
4129 19 Jan 22 peter 64       type of two objects, we need to simulate double dispatch and
4129 19 Jan 22 peter 65       this class is helping the Matrix and MatrixMutable to implement
4129 19 Jan 22 peter 66       that.
4129 19 Jan 22 peter 67
4129 19 Jan 22 peter 68       In brief, it's constructed inside lhs (being a Matrix), taking
4129 19 Jan 22 peter 69       and storing a reference to lhs. It then calls operator() taking
4129 19 Jan 22 peter 70       rhs. It then uses the Visitor pattern calling rhs::visit(Mover),
4129 19 Jan 22 peter 71       where visit is a virtual function having a special implemntation
4129 19 Jan 22 peter 72       for Matrix where it's using fast move assignment, whereas for
4129 19 Jan 22 peter 73       other rhs types move assignment can't be used (because the rhs
4129 19 Jan 22 peter 74       doesn't own its data) and we use copy assignment instead.
4129 19 Jan 22 peter 75      */
4129 19 Jan 22 peter 76     class Mover
4129 19 Jan 22 peter 77     {
4129 19 Jan 22 peter 78     public:
4129 19 Jan 22 peter 79       Mover(Matrix& lhs);
4129 19 Jan 22 peter 80       void operator()(MatrixMutable& rhs);
4129 19 Jan 22 peter 81       Matrix& lhs(void);
4129 19 Jan 22 peter 82       void copy_assign(const MatrixMutable& rhs);
4129 19 Jan 22 peter 83     private:
4129 19 Jan 22 peter 84       Matrix& lhs_;
4129 19 Jan 22 peter 85     };
4129 19 Jan 22 peter 86   }
4129 19 Jan 22 peter 87   /// \endcond
4129 19 Jan 22 peter 88
42 26 Feb 04 jari 89   ///
767 22 Feb 07 peter 90   /// @brief Interface to GSL matrix.
42 26 Feb 04 jari 91   ///
767 22 Feb 07 peter 92   /// For the time being 'double' is the only type supported.
767 22 Feb 07 peter 93   ///
420 02 Dec 05 jari 94   /// \par[File streams] Reading and writing vectors to file streams
420 02 Dec 05 jari 95   /// are of course supported. These are implemented without using GSL
420 02 Dec 05 jari 96   /// functionality, and thus binary read and write to streams are not
420 02 Dec 05 jari 97   /// supported.
420 02 Dec 05 jari 98   ///
420 02 Dec 05 jari 99   /// @note All GSL matrix related functions are not implement but
420 02 Dec 05 jari 100   /// most functionality defined for GSL matrices can be achieved with
420 02 Dec 05 jari 101   /// this interface class. Most notable GSL functionality not
4141 01 Feb 22 peter 102   /// supported are no binary file support.
420 02 Dec 05 jari 103   ///
4119 07 Nov 21 peter 104   class Matrix : public MatrixMutable
420 02 Dec 05 jari 105   {
420 02 Dec 05 jari 106   public:
1064 10 Feb 08 peter 107     /**
614 31 Aug 06 jari 108        @brief The default constructor.
2689 27 Feb 12 peter 109
1822 24 Feb 09 jari 110        This constructor does not initialize underlying (essential)
614 31 Aug 06 jari 111        structures.
614 31 Aug 06 jari 112     */
1121 22 Feb 08 peter 113     Matrix(void);
12 19 Jun 03 daniel 114
754 17 Feb 07 jari 115     /**
754 17 Feb 07 jari 116        \brief Constructor allocating memory space for \a r times \a c
754 17 Feb 07 jari 117        elements, and sets all elements to \a init_value.
754 17 Feb 07 jari 118
1890 02 Apr 09 peter 119        If \a r is zero \a c must be zero and vice versa.
1890 02 Apr 09 peter 120
754 17 Feb 07 jari 121        \throw GSL_error if memory allocation fails.
754 17 Feb 07 jari 122     */
1121 22 Feb 08 peter 123     Matrix(const size_t& r, const size_t& c, double init_value=0);
420 02 Dec 05 jari 124
754 17 Feb 07 jari 125     /**
754 17 Feb 07 jari 126        \brief The copy constructor.
754 17 Feb 07 jari 127
754 17 Feb 07 jari 128        \throw A GSL_error is indirectly thrown if memory allocation
754 17 Feb 07 jari 129        fails.
754 17 Feb 07 jari 130     */
1121 22 Feb 08 peter 131     Matrix(const Matrix&);
12 19 Jun 03 daniel 132
3605 27 Jan 17 peter 133     /**
4129 19 Jan 22 peter 134        Copy from the base class.
4119 07 Nov 21 peter 135
4119 07 Nov 21 peter 136        \since New in yat 0.20
4119 07 Nov 21 peter 137      */
4124 13 Jan 22 peter 138     explicit Matrix(const MatrixBase&);
4119 07 Nov 21 peter 139
4119 07 Nov 21 peter 140     /**
4174 18 May 22 peter 141        \brief Destructor
4174 18 May 22 peter 142      */
4174 18 May 22 peter 143     virtual ~Matrix(void);
4174 18 May 22 peter 144
4174 18 May 22 peter 145     /**
3605 27 Jan 17 peter 146        Constructor from a matrix expression. A matrix expression is
3605 27 Jan 17 peter 147        result \c operator+, \c operator-, and operator*, or
3605 27 Jan 17 peter 148        combinations of them.
3605 27 Jan 17 peter 149
3605 27 Jan 17 peter 150        A typical usage looks like
3605 27 Jan 17 peter 151
3605 27 Jan 17 peter 152        \code
3605 27 Jan 17 peter 153        Matrix A( B * C + D)
3605 27 Jan 17 peter 154        \endcode
3605 27 Jan 17 peter 155
3605 27 Jan 17 peter 156        where B, C, and D are all instances of class Matrix.
3605 27 Jan 17 peter 157
4118 24 Oct 21 peter 158        Typically MatrixExpression only exists as rvalue, and this
4118 24 Oct 21 peter 159        operator is not called but the rvalue version,
4118 24 Oct 21 peter 160        operator=(MatrixExpression<T>&&)
3605 27 Jan 17 peter 161
4118 24 Oct 21 peter 162        Invalidates, references, iterators and views.
4118 24 Oct 21 peter 163
3605 27 Jan 17 peter 164        \since new in yat 0.15
3605 27 Jan 17 peter 165      */
3605 27 Jan 17 peter 166     template<class T>
3651 28 Jun 17 peter 167     Matrix(const MatrixExpression<T>& other);
3605 27 Jan 17 peter 168
754 17 Feb 07 jari 169     /**
3585 19 Jan 17 peter 170        \brief The move constructor.
3596 21 Jan 17 peter 171
4118 24 Oct 21 peter 172        Invalidates references, iterators, and views.
4118 24 Oct 21 peter 173
3596 21 Jan 17 peter 174        \since new in yat 0.15
3585 19 Jan 17 peter 175     */
3691 14 Sep 17 peter 176     Matrix(Matrix&&) noexcept;
3605 27 Jan 17 peter 177
3605 27 Jan 17 peter 178     /**
4119 07 Nov 21 peter 179        Not implemented yet
4119 07 Nov 21 peter 180
4119 07 Nov 21 peter 181        \since New in yat 0.20
4119 07 Nov 21 peter 182      */
4124 13 Jan 22 peter 183     Matrix(MatrixMutable&&);
4119 07 Nov 21 peter 184
4119 07 Nov 21 peter 185     /**
3605 27 Jan 17 peter 186        Move constructor from a matrix expression. A matrix expression is
3605 27 Jan 17 peter 187        result \c operator+, \c operator-, and operator*, or
3605 27 Jan 17 peter 188        combinations of them.
3605 27 Jan 17 peter 189
3605 27 Jan 17 peter 190        A typical usage looks like
3605 27 Jan 17 peter 191
3605 27 Jan 17 peter 192        \code
3605 27 Jan 17 peter 193        Matrix A( B * C + D)
3605 27 Jan 17 peter 194        \endcode
3605 27 Jan 17 peter 195
3605 27 Jan 17 peter 196        where B, C, and D are all instances of class Matrix.
3605 27 Jan 17 peter 197
4118 24 Oct 21 peter 198        Invalidates, references, iterators and views.
4118 24 Oct 21 peter 199
3605 27 Jan 17 peter 200        \since new in yat 0.15
3605 27 Jan 17 peter 201     */
3605 27 Jan 17 peter 202     template<class T>
3605 27 Jan 17 peter 203     Matrix(MatrixExpression<T>&& other)
4129 19 Jan 22 peter 204       : m_(nullptr)
3605 27 Jan 17 peter 205     {
3653 28 Jun 17 peter 206       other.move(m_);
3605 27 Jan 17 peter 207     }
3605 27 Jan 17 peter 208
3585 19 Jan 17 peter 209
3585 19 Jan 17 peter 210     /**
754 17 Feb 07 jari 211        \brief The istream constructor.
754 17 Feb 07 jari 212
1823 24 Feb 09 jari 213        The std::istream will be interpreted as outlined here:
754 17 Feb 07 jari 214
1823 24 Feb 09 jari 215        Missing values, i.e. empty elements, are treated as NaN values
1823 24 Feb 09 jari 216        (std::numeric_limits<double>::quiet_NaN() to be specific).
1823 24 Feb 09 jari 217
1823 24 Feb 09 jari 218        Matrix rows are separated with the new line character.
1823 24 Feb 09 jari 219
1823 24 Feb 09 jari 220        Column element separation has two modes depending on the value
1823 24 Feb 09 jari 221        of \a sep.
1823 24 Feb 09 jari 222
1823 24 Feb 09 jari 223        - If \a sep is the default '\\0' value then column elements are
1823 24 Feb 09 jari 224        separated with white space characters except the new line
1823 24 Feb 09 jari 225        character. Multiple sequential white space characters are
1823 24 Feb 09 jari 226        treated as one separator.
1823 24 Feb 09 jari 227
1823 24 Feb 09 jari 228        - Setting \a sep to something else than the default value will
1823 24 Feb 09 jari 229        change the behaviour to use the \a sep character as the
1823 24 Feb 09 jari 230        separator between column elements. Multiple sequential \a sep
1823 24 Feb 09 jari 231        characters will be treated as separating elements with missing
1823 24 Feb 09 jari 232        values.
1823 24 Feb 09 jari 233
1823 24 Feb 09 jari 234        End of input is the end of file marker and this treatment
1823 24 Feb 09 jari 235        cannot be redefined using the provided API.
1823 24 Feb 09 jari 236
759 19 Feb 07 jari 237        \throw GSL_error if memory allocation fails, IO_error if
759 19 Feb 07 jari 238        unexpected input is found in the input stream.
754 17 Feb 07 jari 239     */
3845 13 Sep 19 peter 240     explicit Matrix(std::istream &, char sep='\0');
12 19 Jun 03 daniel 241
754 17 Feb 07 jari 242     /**
4129 19 Jan 22 peter 243        @return A pointer to the internal GSL matrix.
4129 19 Jan 22 peter 244      */
4129 19 Jan 22 peter 245     gsl_matrix* gsl_matrix_p(void);
4129 19 Jan 22 peter 246
4129 19 Jan 22 peter 247     /**
4129 19 Jan 22 peter 248        @return A const pointer to the internal GSL matrix.
4129 19 Jan 22 peter 249      */
4129 19 Jan 22 peter 250     const gsl_matrix* gsl_matrix_p(void) const;
4129 19 Jan 22 peter 251
4129 19 Jan 22 peter 252     /**
1121 22 Feb 08 peter 253        \brief Resize Matrix
2689 27 Feb 12 peter 254
811 16 Mar 07 peter 255        All elements are set to @a init_value.
811 16 Mar 07 peter 256
1890 02 Apr 09 peter 257        If \a r is zero \a c must be zero and vice versa.
1890 02 Apr 09 peter 258
2689 27 Feb 12 peter 259        \note underlying GSL matrix is destroyed and views and
2689 27 Feb 12 peter 260        iterators are invalidated
811 16 Mar 07 peter 261     */
1890 02 Apr 09 peter 262     void resize(size_t r, size_t c, double init_value=0);
808 15 Mar 07 peter 263
754 17 Feb 07 jari 264     /**
754 17 Feb 07 jari 265        \brief Transpose the matrix.
754 17 Feb 07 jari 266
2689 27 Feb 12 peter 267        \note Unless matrix is square, views and iterators are
2689 27 Feb 12 peter 268        invalidated.
2689 27 Feb 12 peter 269
754 17 Feb 07 jari 270        \throw GSL_error if memory allocation fails for the new
754 17 Feb 07 jari 271        transposed matrix.
754 17 Feb 07 jari 272     */
420 02 Dec 05 jari 273     void transpose(void);
12 19 Jun 03 daniel 274
755 18 Feb 07 jari 275     /**
754 17 Feb 07 jari 276        \brief The assignment operator.
754 17 Feb 07 jari 277
4129 19 Jan 22 peter 278        \note views and iterators are invalidated.
2689 27 Feb 12 peter 279
1121 22 Feb 08 peter 280        \return A const reference to the resulting Matrix.
754 17 Feb 07 jari 281     */
1121 22 Feb 08 peter 282     const Matrix& operator=(const Matrix& other);
4129 19 Jan 22 peter 283     using MatrixMutable::operator=; // access to =(const MatrixBase&)
12 19 Jun 03 daniel 284
3605 27 Jan 17 peter 285     /**
4129 19 Jan 22 peter 286        \brief Move assignment operatorOnly reason they mention those 1e-5 peaks
4119 07 Nov 21 peter 287
4119 07 Nov 21 peter 288
3605 27 Jan 17 peter 289
4118 24 Oct 21 peter 290        Invalidates references, iterators, and views.
4118 24 Oct 21 peter 291
3596 21 Jan 17 peter 292        \since new in yat 0.15
3585 19 Jan 17 peter 293      */
3585 19 Jan 17 peter 294     Matrix& operator=(Matrix&& other);
3605 27 Jan 17 peter 295
4129 19 Jan 22 peter 296     // to get access to operator*=(double)
4119 07 Nov 21 peter 297     using MatrixMutable::operator*=;
754 17 Feb 07 jari 298
773 01 Mar 07 jari 299     /**
1822 24 Feb 09 jari 300        \brief Multiply and assignment operator.
754 17 Feb 07 jari 301
1121 22 Feb 08 peter 302        \return Const reference to the resulting Matrix.
754 17 Feb 07 jari 303
2689 27 Feb 12 peter 304        \note invalidates views and iterators
2689 27 Feb 12 peter 305
754 17 Feb 07 jari 306        \throw GSL_error if memory allocation fails.
754 17 Feb 07 jari 307     */
1121 22 Feb 08 peter 308     const Matrix& operator*=(const Matrix&);
420 02 Dec 05 jari 309
4129 19 Jan 22 peter 310   protected:
4129 19 Jan 22 peter 311     /**
4129 19 Jan 22 peter 312        Allocate necessary memory in underlying GSL matrix and copy
4129 19 Jan 22 peter 313        data from \c rhs. It is safe to copy from an \c rhs that
4129 19 Jan 22 peter 314        overlap in underlying data structure.
4129 19 Jan 22 peter 315      */
4129 19 Jan 22 peter 316     void copy_assign(const gsl_matrix* rhs);
4129 19 Jan 22 peter 317     using MatrixMutable::copy_assign;
4129 19 Jan 22 peter 318
4129 19 Jan 22 peter 319     /**
4129 19 Jan 22 peter 320        If rhs is a Matrix, move assignment is used, i.e., the
4129 19 Jan 22 peter 321        operation is cheap (constant time) only swapping a two
4129 19 Jan 22 peter 322        pointers. If \c rhs is a MatrixView, which doesn't own its
4129 19 Jan 22 peter 323        data, data is aligned by copying instead. Dimensions of lhs and
4129 19 Jan 22 peter 324        rhs do not need to be the same.
4129 19 Jan 22 peter 325      */
4129 19 Jan 22 peter 326     void move_assign(MatrixMutable&& rhs);
4129 19 Jan 22 peter 327
4129 19 Jan 22 peter 328     /**
4129 19 Jan 22 peter 329        Move assign data in \c rhs into this. Currently, the passed \c
4129 19 Jan 22 peter 330        rhs must own its data and move assignment is therefore always
4129 19 Jan 22 peter 331        used. Dimensions of lhs and rhs do not need to be the same.
4129 19 Jan 22 peter 332      */
4129 19 Jan 22 peter 333     void move_assign(gsl_matrix*&& rhs);
4129 19 Jan 22 peter 334
420 02 Dec 05 jari 335   private:
4129 19 Jan 22 peter 336     gsl_matrix* m_;
4129 19 Jan 22 peter 337
754 17 Feb 07 jari 338     /**
754 17 Feb 07 jari 339        \brief Create a new copy of the internal GSL matrix.
754 17 Feb 07 jari 340
754 17 Feb 07 jari 341        Necessary memory for the new GSL matrix is allocated and the
754 17 Feb 07 jari 342        caller is responsible for freeing the allocated memory.
754 17 Feb 07 jari 343
754 17 Feb 07 jari 344        \return A pointer to a copy of the internal GSL matrix.
754 17 Feb 07 jari 345
754 17 Feb 07 jari 346        \throw GSL_error if memory cannot be allocated for the new
1822 24 Feb 09 jari 347        copy, or if dimensions mismatch.
754 17 Feb 07 jari 348     */
420 02 Dec 05 jari 349     gsl_matrix* create_gsl_matrix_copy(void) const;
12 19 Jun 03 daniel 350
4129 19 Jan 22 peter 351     // Used in move_assignment, see docs for Mover class.
4129 19 Jan 22 peter 352     friend class detail::Mover;
4129 19 Jan 22 peter 353     void visit(detail::Mover& mover);
42 26 Feb 04 jari 354   };
12 19 Jun 03 daniel 355
774 01 Mar 07 jari 356   /**
3657 13 Jul 17 peter 357      Calculate inverse of a (square) Matrix using SVD
3657 13 Jul 17 peter 358
3657 13 Jul 17 peter 359      \since New in yat 0.15
3657 13 Jul 17 peter 360
4118 24 Oct 21 peter 361      \note Space for Matrix \a result reallocated and references,
4118 24 Oct 21 peter 362      iterators and views are invalidated.
4118 24 Oct 21 peter 363
3657 13 Jul 17 peter 364      \relates Matrix
3657 13 Jul 17 peter 365    */
4140 29 Jan 22 peter 366   void inverse_svd(const MatrixBase& input, Matrix& result);
3657 13 Jul 17 peter 367
3657 13 Jul 17 peter 368   /**
1121 22 Feb 08 peter 369      \brief Create a Matrix \a flag indicating NaN's in another Matrix
774 01 Mar 07 jari 370      \a templat.
774 01 Mar 07 jari 371
1121 22 Feb 08 peter 372      The \a flag Matrix is changed to contain 1's and 0's only. A 1
1121 22 Feb 08 peter 373      means that the corresponding element in the \a templat Matrix is
774 01 Mar 07 jari 374      valid and a zero means that the corresponding element is a NaN.
774 01 Mar 07 jari 375
4118 24 Oct 21 peter 376      \note If sizes mismatch, space for Matrix \a flag is reallocated
4118 24 Oct 21 peter 377      and references, iterators, and views are invalidated.
774 01 Mar 07 jari 378
1121 22 Feb 08 peter 379      \return True if the \a templat Matrix contains at least one NaN.
1887 31 Mar 09 peter 380
1887 31 Mar 09 peter 381      \relates Matrix
774 01 Mar 07 jari 382   */
4138 28 Jan 22 peter 383   bool nan(const MatrixBase& templat, Matrix& flag);
774 01 Mar 07 jari 384
3651 28 Jun 17 peter 385   /////  template implemantions
3651 28 Jun 17 peter 386   template<class T>
3651 28 Jun 17 peter 387   Matrix::Matrix(const MatrixExpression<T>& other)
3651 28 Jun 17 peter 388   {
3651 28 Jun 17 peter 389     // In principle, we could implement this as a move constructor
3651 28 Jun 17 peter 390     // and steal the data (rather than copy), but it's a bit hackish
3651 28 Jun 17 peter 391     // to workaround the constness, so users who are eager for speed
3651 28 Jun 17 peter 392     // should enable rvalues and get access to the faster move
3651 28 Jun 17 peter 393     // constructor below.
3651 28 Jun 17 peter 394
3651 28 Jun 17 peter 395     try {
3651 28 Jun 17 peter 396       detail::copy(m_, other.gsl_matrix_p());
3651 28 Jun 17 peter 397     }
3651 28 Jun 17 peter 398     catch (GSL_error& e) {
3651 28 Jun 17 peter 399       detail::deallocate(m_);
3651 28 Jun 17 peter 400       throw e;
3651 28 Jun 17 peter 401     }
3651 28 Jun 17 peter 402   }
3651 28 Jun 17 peter 403
687 16 Oct 06 jari 404 }}} // of namespace utility, yat, and theplu
12 19 Jun 03 daniel 405
420 02 Dec 05 jari 406 #endif