Changeset 4129


Ignore:
Timestamp:
Jan 19, 2022, 9:57:41 AM (20 months ago)
Author:
Peter
Message:

refs #202. Adding classes for Matrix View (and corresponding const version)

Location:
trunk
Files:
5 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/test/Makefile.am

    r4116 r4129  
    8989  test/matrix_lookup_weighted.test test/matrix.test \
    9090  test/matrix_expression.test \
     91  test/matrix_view.test \
    9192  test/matrix_weighted.test test/merge.test \
    9293  test/merge_iterator.test \

  • trunk/yat/utility/BasicMatrix.cc

    r4124 r4129  
    130130  void reallocate(gsl_matrix*& p, size_t row, size_t col)
    131131  {
     132    assert(p==nullptr || p->owner);
    132133    assert(static_cast<bool>(row) == static_cast<bool>(col));
    133134    if (row!=rows(p) || col!=columns(p)) {
    … …  
    147148
    148149
     150  const double* end(const gsl_matrix* p)
     151  {
     152    assert(p);
     153    assert(p->block);
     154    return p->data + p->block->size;
     155  }
     156
     157
     158  bool overlap(const gsl_matrix* lhs, const gsl_matrix* rhs)
     159  {
     160    if (!lhs || !rhs)
     161      return false;
     162
     163    return (rhs->data >= lhs->data && rhs->data < end(lhs)) ||
     164      (lhs->data >= rhs->data && lhs->data < end(rhs));
     165  }
     166
     167
    149168}}}} // of namespace detail, utility, yat and thep

  • trunk/yat/utility/BasicMatrix.h

    r4124 r4129  
    9494    void copy(gsl_matrix*& dest, const gsl_matrix* src);
    9595
     96    // return one passed end of underlying data
     97    const double* end(const gsl_matrix* p);
     98
    9699    // If both own their data, swap pointers; otherwise call copy
     100    // can practically only be called in Matrix, so move it there
    97101    void move_if_owner(gsl_matrix*& dest, gsl_matrix*& src);
    98102
    … …  
    111115    // if \a p is NULL, return 0; otherwise return size1
    112116    size_t rows(const gsl_matrix* p);
     117
     118    /**
     119       \return true if underlying ranges [data, data+size) overlap
     120     */
     121    bool overlap(const gsl_matrix* lhs, const gsl_matrix* rhs);
    113122  }
    114123

  • trunk/yat/utility/Makefile.am

    r4119 r4129  
    3737  yat/utility/Matrix.cc \
    3838  yat/utility/MatrixBase.cc \
     39  yat/utility/MatrixConstView.cc \
    3940  yat/utility/MatrixMutable.cc \
     41  yat/utility/MatrixView.cc \
    4042  yat/utility/MatrixWeighted.cc \
    4143  yat/utility/NNI.cc \
    … …  
    98100  $(srcdir)/yat/utility/Matrix.h \
    99101  $(srcdir)/yat/utility/MatrixBase.h \
     102  $(srcdir)/yat/utility/MatrixConstView.h \
     103  $(srcdir)/yat/utility/MatrixExpression.h \
    100104  $(srcdir)/yat/utility/MatrixMutable.h \
    101   $(srcdir)/yat/utility/MatrixExpression.h \
     105  $(srcdir)/yat/utility/MatrixView.h \
    102106  $(srcdir)/yat/utility/MatrixWeighted.h \
    103107  $(srcdir)/yat/utility/merge.h \

  • trunk/yat/utility/Matrix.cc

    r4124 r4129  
    5454
    5555  Matrix::Matrix(void)
     56    : m_(nullptr)
    5657  {
    5758  }
    … …  
    5960
    6061  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
     62    : m_(nullptr)
    6163  {
    6264    assert(static_cast<bool>(r) == static_cast<bool>(c));
    … …  
    6870
    6971  Matrix::Matrix(const Matrix& o)
     72    : m_(nullptr)
    7073  {
    7174    m_ = detail::create_copy(o.gsl_matrix_p());
    … …  
    7477
    7578  Matrix::Matrix(const MatrixBase& o)
     79    : m_(nullptr)
    7680  {
    7781    m_ = detail::create_copy(o.gsl_matrix_p());
    … …  
    8084
    8185  Matrix::Matrix(Matrix&& other) noexcept
     86    : m_(nullptr)
    8287  {
    8388    std::swap(m_, other.m_);
    … …  
    8691
    8792  Matrix::Matrix(MatrixMutable&& other)
    88   {
    89     move_assign(std::move(other));
     93    : m_(nullptr)
     94  {
     95    assert(0);
     96    //move_assign(std::move(other));
    9097  }
    9198
    … …  
    93100  // Constructor that gets data from istream
    94101  Matrix::Matrix(std::istream& is, char sep)
     102    : m_(nullptr)
    95103  {
    96104    if (!is.good())
    … …  
    129137      std::copy(data_matrix[i].begin(), data_matrix[i].end(), begin_row(i));
    130138    }
     139  }
     140
     141
     142  void Matrix::copy_assign(const gsl_matrix* rhs)
     143  {
     144    // copy via blas_result and then swap result to m_
     145    detail::copy(blas_result_, rhs);
     146    std::swap(blas_result_, m_);
     147  }
     148
     149
     150  gsl_matrix* Matrix::gsl_matrix_p(void)
     151  {
     152    return m_;
     153  }
     154
     155
     156  const gsl_matrix* Matrix::gsl_matrix_p(void) const
     157  {
     158    return m_;
     159  }
     160
     161
     162  void Matrix::move_assign(MatrixMutable&& rhs)
     163  {
     164    detail::Mover mover(*this);
     165    mover(rhs);
     166  }
     167
     168
     169  void Matrix::move_assign(gsl_matrix*&& rhs)
     170  {
     171    // if there is a use case for when rhs is a view, we could allow
     172    // that with an if statement, but ATM it seems a waste.
     173    assert(!rhs || rhs->owner);
     174    std::swap(m_, rhs);
    131175  }
    132176
    … …  
    171215
    172216
     217  void Matrix::visit(detail::Mover& mover)
     218  {
     219    std::swap(mover.lhs().m_, m_);
     220  }
     221
     222
    173223  const Matrix& Matrix::operator=( const Matrix& other )
    174224  {
    … …  
    178228
    179229
    180   const Matrix& Matrix::operator=( const MatrixBase& other )
    181   {
    182     copy_assign(other);
    183     return *this;
    184   }
    185 
    186 
    187   Matrix& Matrix::operator=(MatrixMutable&& other)
    188   {
    189     move_assign(std::move(other));
    190     return *this;
    191   }
    192 
    193 
    194230  Matrix& Matrix::operator=(Matrix&& other)
    195231  {
    196     move_assign(std::move(other));
     232    std::swap(m_, other.m_);
    197233    return *this;
    198234  }
    … …  
    244280
    245281
     282  namespace detail {
     283
     284    Mover::Mover(Matrix& m)
     285      : lhs_(m) {}
     286
     287
     288    void Mover::operator()(MatrixMutable& rhs)
     289    {
     290      rhs.visit(*this);
     291    }
     292
     293
     294    void Mover::copy_assign(const MatrixMutable& rhs)
     295    {
     296      lhs_.copy_assign(rhs);
     297    }
     298
     299
     300    Matrix& Mover::lhs(void)
     301    {
     302      return lhs_;
     303    }
     304  }
     305
    246306}}} // of namespace utility, yat and thep

  • trunk/yat/utility/Matrix.h

    r4126 r4129  
    2626  along with yat. If not, see <http://www.gnu.org/licenses/>.
    2727*/
     28
     29#include <cassert>
    2830
    2931#include "config_public.h"
    … …  
    4951
    5052  class VectorBase;
     53  class Matrix;
     54
     55  /// \cond IGNORE_DOXYGEN
     56
     57  namespace detail {
     58
     59    /*
     60      Class involved in move assignment of MatrixMutable to determine
     61      when both lhs and rhs are a Matrix. When the behavoiur only
     62      depends on the type of one object, dynamic polymorphism with a
     63      virtual function suffices; when the behaviour depends on two
     64      type of two objects, we need to simulate double dispatch and
     65      this class is helping the Matrix and MatrixMutable to implement
     66      that.
     67
     68      In brief, it's constructed inside lhs (being a Matrix), taking
     69      and storing a reference to lhs. It then calls operator() taking
     70      rhs. It then uses the Visitor pattern calling rhs::visit(Mover),
     71      where visit is a virtual function having a special implemntation
     72      for Matrix where it's using fast move assignment, whereas for
     73      other rhs types move assignment can't be used (because the rhs
     74      doesn't own its data) and we use copy assignment instead.
     75     */
     76    class Mover
     77    {
     78    public:
     79      Mover(Matrix& lhs);
     80      void operator()(MatrixMutable& rhs);
     81      Matrix& lhs(void);
     82      void copy_assign(const MatrixMutable& rhs);
     83    private:
     84      Matrix& lhs_;
     85    };
     86  }
     87  /// \endcond
    5188
    5289  ///
    … …  
    97134
    98135    /**
    99        Not implemented yet
     136       Copy from the base class.
    100137
    101138       \since New in yat 0.20
    … …  
    162199    template<class T>
    163200    Matrix(MatrixExpression<T>&& other)
     201      : m_(nullptr)
    164202    {
    165203      other.move(m_);
    … …  
    200238
    201239    /**
     240       @return A pointer to the internal GSL matrix.
     241     */
     242    gsl_matrix* gsl_matrix_p(void);
     243
     244    /**
     245       @return A const pointer to the internal GSL matrix.
     246     */
     247    const gsl_matrix* gsl_matrix_p(void) const;
     248
     249    /**
    202250       \brief Resize Matrix
    203251
    … …  
    225273       \brief The assignment operator.
    226274
    227        \note If lhs needs to be resized, views and iterators are
    228        invalidated.
     275       \note views and iterators are invalidated.
    229276
    230277       \return A const reference to the resulting Matrix.
    231278    */
    232279    const Matrix& operator=(const Matrix& other);
    233 
    234     /**
    235        Not implemented yet
    236 
    237        \since New in yat 0.20
    238      */
    239     const Matrix& operator=(const MatrixBase& other);
    240 
    241     /**
    242        Assignment from a matrix expression. A matrix expression is the
    243        result of \c operator+, \c operator-, and operator*, or
    244        combinations of them.
    245 
    246        A typical usage looks like
    247 
    248        \code
    249        A = B * C + D;
    250        \endcode
    251 
    252        where A, B, C, and D are all instances of class Matrix.
    253 
    254        Typically MatrixExpression only exists as rvalue, and this
    255        operator is not called but the rvalue version,
    256        operator=(MatrixExpression<T>&&)
    257 
    258        Invalidates, references, iterators and views.
     280    using MatrixMutable::operator=; // access to =(const MatrixBase&)
     281
     282    /**
     283       \brief Move assignment operatorOnly reason they mention those 1e-5 peaks
     284
     285
     286
     287       Invalidates references, iterators, and views.
    259288
    260289       \since new in yat 0.15
    261290     */
    262     template<class T>
    263     Matrix& operator=(const MatrixExpression<T>& rhs);
    264 
    265 
    266     /**
    267        \brief Move assignment operator
    268 
    269        Invalidates references, iterators, and views.
    270 
    271        \since new in yat 0.15
    272      */
    273291    Matrix& operator=(Matrix&& other);
    274292
    275     /**
    276        Not implemented yet
    277 
    278        \since New in yat 0.20
    279      */
    280     Matrix& operator=(MatrixMutable&& other);
    281 
    282     /**
    283        Move assignment from a matrix expression. A matrix expression
    284        is the result of \c operator+, \c operator-, and operator*, or
    285        combinations of them.
    286 
    287        A typical usage looks like
    288 
    289        \code
    290        A = B * C + D;
    291        \endcode
    292 
    293        where B, C, and D are all instances of class Matrix.
    294 
    295        Invalidates references, iterators, and views.
    296 
    297        \since new in yat 0.15
    298      */
    299     template<class T>
    300     Matrix& operator=(MatrixExpression<T>&& rhs)
    301     {
    302       // This function steals data from rhs into this, and give the old
    303       // m_ in return to rhs so destructor of rhs takes care of
    304       // deallocation.
    305       rhs.move(m_);
    306       return *this;
    307     }
    308 
     293    // to get access to operator*=(double)
    309294    using MatrixMutable::operator*=;
    310295
    … …  
    320305    const Matrix& operator*=(const Matrix&);
    321306
     307  protected:
     308    /**
     309       Allocate necessary memory in underlying GSL matrix and copy
     310       data from \c rhs. It is safe to copy from an \c rhs that
     311       overlap in underlying data structure.
     312     */
     313    void copy_assign(const gsl_matrix* rhs);
     314    using MatrixMutable::copy_assign;
     315
     316    /**
     317       If rhs is a Matrix, move assignment is used, i.e., the
     318       operation is cheap (constant time) only swapping a two
     319       pointers. If \c rhs is a MatrixView, which doesn't own its
     320       data, data is aligned by copying instead. Dimensions of lhs and
     321       rhs do not need to be the same.
     322     */
     323    void move_assign(MatrixMutable&& rhs);
     324
     325    /**
     326       Move assign data in \c rhs into this. Currently, the passed \c
     327       rhs must own its data and move assignment is therefore always
     328       used. Dimensions of lhs and rhs do not need to be the same.
     329     */
     330    void move_assign(gsl_matrix*&& rhs);
     331
    322332  private:
     333    gsl_matrix* m_;
     334
    323335    /**
    324336       \brief Create a new copy of the internal GSL matrix.
    … …  
    334346    gsl_matrix* create_gsl_matrix_copy(void) const;
    335347
     348    // Used in move_assignment, see docs for Mover class.
     349    friend class detail::Mover;
     350    void visit(detail::Mover& mover);
    336351  };
    337352
    … …  
    384399  }
    385400
    386 
    387   template<class T>
    388   Matrix& Matrix::operator=(const MatrixExpression<T>& rhs)
    389   {
    390     // Be careful because *this might be embedded in rhs.
    391     rhs.get(blas_result_);
    392     std::swap(m_, blas_result_);
    393     return *this;
    394   }
    395 
    396 
    397401}}} // of namespace utility, yat, and theplu
    398402

  • trunk/yat/utility/MatrixMutable.cc

    r4124 r4129  
    5454
    5555  MatrixMutable::MatrixMutable(void)
    56     : m_(nullptr), blas_result_(nullptr)
     56    : blas_result_(nullptr)
    5757  {
    5858  }
    … …  
    6161  MatrixMutable::~MatrixMutable(void)
    6262  {
    63     detail::deallocate(m_);
    6463    detail::deallocate(blas_result_);
    6564  }
    … …  
    9897
    9998
    100   void MatrixMutable::copy_assign(const MatrixBase& other)
    101   {
    102     if (this!=&other) {
    103       detail::copy(blas_result_, other.gsl_matrix_p());
    104       std::swap(blas_result_, m_);
    105     }
    106   }
    107 
    108 
    109   void MatrixMutable::delete_allocated_memory(void)
    110   {
    111     detail::deallocate(blas_result_);
    112   }
    113 
    114 
    11599  void MatrixMutable::div(const MatrixBase& other)
    116100  {
    … …  
    137121  {
    138122    return row_iterator(&(*this)(i,0)+columns(), 1);
    139   }
    140 
    141 
    142   gsl_matrix* MatrixMutable::gsl_matrix_p(void)
    143   {
    144     return m_;
    145   }
    146 
    147 
    148   const gsl_matrix* MatrixMutable::gsl_matrix_p(void) const
    149   {
    150     return m_;
    151   }
    152 
    153 
    154   void MatrixMutable::move_assign(MatrixMutable&& other)
    155   {
    156     detail::move_if_owner(m_, other.m_);
    157123  }
    158124
    … …  
    279245
    280246
     247  void MatrixMutable::copy_assign(const MatrixBase& other)
     248  {
     249    if (this != &other)
     250      copy_assign(other.gsl_matrix_p());
     251  }
     252
     253
     254  void MatrixMutable::move_assign(MatrixMutable&& other)
     255  {
     256    copy_assign(other);
     257  }
     258
     259
     260  void MatrixMutable::visit(detail::Mover& mover)
     261  {
     262    mover.copy_assign(*this);
     263  }
     264
     265
    281266  void swap(MatrixMutable& a, MatrixMutable& b)
    282267  {

  • trunk/yat/utility/MatrixMutable.h

    r4126 r4129  
    4747  class VectorBase;
    4848
     49  namespace detail {
     50    class Mover;
     51  }
     52
    4953  /**
    5054     Interface class for mutable matrices including Matrix and
    … …  
    160164    /// @return A pointer to the internal GSL matrix.
    161165    ///
    162     gsl_matrix* gsl_matrix_p(void);
     166    virtual gsl_matrix* gsl_matrix_p(void)=0;
    163167
    164168    /**
    165169       @return A const pointer to the internal GSL matrix.
    166170    */
    167     const gsl_matrix* gsl_matrix_p(void) const;
     171    virtual const gsl_matrix* gsl_matrix_p(void) const=0;
    168172
    169173    /**
    … …  
    259263    {
    260264      rhs.move(blas_result_);
    261       detail::move_if_owner(m_, blas_result_);
     265      move_assign(std::move(blas_result_));
    262266      return *this;
    263267    }
    … …  
    294298
    295299      rhs.move(blas_result_);
    296       detail::move_if_owner(m_, blas_result_);
     300      move_assign(std::move(blas_result_));
    297301      return *this;
    298302    }
    … …  
    376380  protected:
    377381    /**
    378        copy data from \c other to \c this
    379      */
    380     void copy_assign(const MatrixBase& other);
    381 
    382     /**
    383        If both this and other own their data, move data from \c other
    384        to \c this. Otherwise copy.
    385      */
    386     void move_assign(MatrixMutable&& other);
    387 
    388     /**
    389        pointer to underlying GSL matrix structure
    390      */
    391     gsl_matrix* m_;
    392 
    393     /**
    394382       blas_result_ is used to temporarily store result in BLAS calls.
    395383       Memory is not allocated for blas_result_ until it is needed.
    … …  
    397385    gsl_matrix* blas_result_;
    398386
     387    /// Behaves like operator=(const MatrixBase&)
     388    void copy_assign(const MatrixBase& other);
     389
     390    /// Behaves like operator=(const MatrixBase&)
     391    virtual void copy_assign(const gsl_matrix* rhs)=0;
     392
     393    /**
     394       Behaves like operator=(MatrixMutable&&)
     395     */
     396    virtual void move_assign(MatrixMutable&& other);
     397
     398    /**
     399       Behaves like operator=(MatrixMutable&&)
     400     */
     401    virtual void move_assign(gsl_matrix*&& rhs)=0;
    399402  private:
    400     /**
    401        \brief Clear all dynamically allocated memory.
    402 
    403        Internal utility function.
    404     */
    405     void delete_allocated_memory(void);
    406 
     403    friend class detail::Mover;
     404    virtual void visit(detail::Mover& mover);
    407405  };
    408406
    … …  
    418416     \throw GSL_error if sizes are not equal.
    419417
    420      \relates Matrix
     418     \relates MatrixMutable
    421419  */
    422420  void swap(MatrixMutable&, MatrixMutable&);
Note: See TracChangeset for help on using the changeset viewer.