yat/utility/BasicMatrix.h

Code
Comments
Other
Rev Date Author Line
3605 27 Jan 17 peter 1 #ifndef _theplu_yat_utility_basic_matrix
3605 27 Jan 17 peter 2 #define _theplu_yat_utility_basic_matrix
3605 27 Jan 17 peter 3
3605 27 Jan 17 peter 4 // $Id$
3605 27 Jan 17 peter 5
3605 27 Jan 17 peter 6 /*
4207 26 Aug 22 peter 7   Copyright (C) 2017, 2022 Peter Johansson
3605 27 Jan 17 peter 8
3605 27 Jan 17 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3605 27 Jan 17 peter 10
3605 27 Jan 17 peter 11   The yat library is free software; you can redistribute it and/or
3605 27 Jan 17 peter 12   modify it under the terms of the GNU General Public License as
3605 27 Jan 17 peter 13   published by the Free Software Foundation; either version 3 of the
3605 27 Jan 17 peter 14   License, or (at your option) any later version.
3605 27 Jan 17 peter 15
3605 27 Jan 17 peter 16   The yat library is distributed in the hope that it will be useful,
3605 27 Jan 17 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3605 27 Jan 17 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3605 27 Jan 17 peter 19   General Public License for more details.
3605 27 Jan 17 peter 20
3605 27 Jan 17 peter 21   You should have received a copy of the GNU General Public License
3605 27 Jan 17 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3605 27 Jan 17 peter 23 */
3605 27 Jan 17 peter 24
3605 27 Jan 17 peter 25 #include <gsl/gsl_matrix.h>
3605 27 Jan 17 peter 26
3605 27 Jan 17 peter 27 #include <cstddef>
3605 27 Jan 17 peter 28
3605 27 Jan 17 peter 29 namespace theplu {
3605 27 Jan 17 peter 30 namespace yat {
3605 27 Jan 17 peter 31 namespace utility {
3605 27 Jan 17 peter 32
3605 27 Jan 17 peter 33   /**
3605 27 Jan 17 peter 34      Something that looks like a Matrix. Is either a Matrix or a
3605 27 Jan 17 peter 35      MatrixExpression (see below)
3605 27 Jan 17 peter 36   */
3605 27 Jan 17 peter 37   template<class Derived>
3605 27 Jan 17 peter 38   class BasicMatrix
3605 27 Jan 17 peter 39   {
3605 27 Jan 17 peter 40   public:
3605 27 Jan 17 peter 41     /**
3605 27 Jan 17 peter 42        Class that inherits from BasicMatrix
3605 27 Jan 17 peter 43      */
3605 27 Jan 17 peter 44     typedef Derived derived_type;
3605 27 Jan 17 peter 45
3605 27 Jan 17 peter 46     /**
3605 27 Jan 17 peter 47        \return element in row \a row and column \a column
3605 27 Jan 17 peter 48      */
3605 27 Jan 17 peter 49     double operator()(size_t row, size_t column) const
3605 27 Jan 17 peter 50     {  return (*static_cast<const Derived*>(this))(row, column);  }
3605 27 Jan 17 peter 51
3605 27 Jan 17 peter 52     /**
3605 27 Jan 17 peter 53        \return underlying data structure
3605 27 Jan 17 peter 54      */
3605 27 Jan 17 peter 55     const gsl_matrix* gsl_matrix_p(void) const
3605 27 Jan 17 peter 56     {  return static_cast<const Derived*>(this)->gsl_matrix_p();  }
3605 27 Jan 17 peter 57
3605 27 Jan 17 peter 58     /**
3605 27 Jan 17 peter 59        \return number of rows
3605 27 Jan 17 peter 60      */
3605 27 Jan 17 peter 61     size_t rows(void) const
3605 27 Jan 17 peter 62     {  return static_cast<const Derived*>(this)->rows();  }
3605 27 Jan 17 peter 63
3605 27 Jan 17 peter 64     /**
3605 27 Jan 17 peter 65        \return number of columns
3605 27 Jan 17 peter 66      */
3605 27 Jan 17 peter 67     size_t columns(void) const
3605 27 Jan 17 peter 68     { return static_cast<const Derived*>(this)->columns();}
3605 27 Jan 17 peter 69
3605 27 Jan 17 peter 70   };
3605 27 Jan 17 peter 71
3651 28 Jun 17 peter 72
3651 28 Jun 17 peter 73   // Some low-level convenience functions to call gsl functionality,
3651 28 Jun 17 peter 74   // which checks if error occurs and throw an exception if error is
3651 28 Jun 17 peter 75   // detected.
3651 28 Jun 17 peter 76   //
3651 28 Jun 17 peter 77   // Functions are placed here so they can be used in both Matrix and
3651 28 Jun 17 peter 78   // MatrixExpression
3651 28 Jun 17 peter 79   namespace detail
3651 28 Jun 17 peter 80   {
3651 28 Jun 17 peter 81     // allocate memory
3651 28 Jun 17 peter 82     // throw if allocation fails
3651 28 Jun 17 peter 83     gsl_matrix* allocate_matrix(size_t row, size_t col);
3651 28 Jun 17 peter 84
3651 28 Jun 17 peter 85     // copy data in src to dest
3651 28 Jun 17 peter 86     // dimensions must agree
3651 28 Jun 17 peter 87     void assign(gsl_matrix* dest, const gsl_matrix* src);
3651 28 Jun 17 peter 88
3651 28 Jun 17 peter 89     // if \a p is NULL, return 0; otherwise return size2
3651 28 Jun 17 peter 90     size_t columns(const gsl_matrix* p);
3651 28 Jun 17 peter 91
3651 28 Jun 17 peter 92     // copy data in src to dest
3651 28 Jun 17 peter 93     // if dimensions do not agree, dest is reallocated
3651 28 Jun 17 peter 94     void copy(gsl_matrix*& dest, const gsl_matrix* src);
3651 28 Jun 17 peter 95
4129 19 Jan 22 peter 96     // return one passed end of underlying data
4129 19 Jan 22 peter 97     const double* end(const gsl_matrix* p);
4129 19 Jan 22 peter 98
4124 13 Jan 22 peter 99     // If both own their data, swap pointers; otherwise call copy
4129 19 Jan 22 peter 100     // can practically only be called in Matrix, so move it there
4124 13 Jan 22 peter 101     void move_if_owner(gsl_matrix*& dest, gsl_matrix*& src);
4124 13 Jan 22 peter 102
3651 28 Jun 17 peter 103     // create a copy of p
3651 28 Jun 17 peter 104     gsl_matrix* create_copy(const gsl_matrix* p);
3651 28 Jun 17 peter 105
3651 28 Jun 17 peter 106     // free pointer and set to NULL
3651 28 Jun 17 peter 107     void deallocate(gsl_matrix*& p);
3651 28 Jun 17 peter 108
3651 28 Jun 17 peter 109     // Similar to realloc in std. Allocate memory for a pointer of
3651 28 Jun 17 peter 110     // size row and col. If passed pointer \a p is already row x col,
3651 28 Jun 17 peter 111     // p is returned, otherwise p is deleted and a new gsl_matrix is
3651 28 Jun 17 peter 112     // allocated.
3651 28 Jun 17 peter 113     void reallocate(gsl_matrix*& p, size_t row, size_t col);
3651 28 Jun 17 peter 114
3651 28 Jun 17 peter 115     // if \a p is NULL, return 0; otherwise return size1
3651 28 Jun 17 peter 116     size_t rows(const gsl_matrix* p);
4129 19 Jan 22 peter 117
4129 19 Jan 22 peter 118     /**
4129 19 Jan 22 peter 119        \return true if underlying ranges [data, data+size) overlap
4129 19 Jan 22 peter 120      */
4129 19 Jan 22 peter 121     bool overlap(const gsl_matrix* lhs, const gsl_matrix* rhs);
3651 28 Jun 17 peter 122   }
3651 28 Jun 17 peter 123
3605 27 Jan 17 peter 124 }}} // of namespace utility, yat, and theplu
3605 27 Jan 17 peter 125
3605 27 Jan 17 peter 126 #endif