Changeset 3623


Ignore:
Timestamp:
Feb 9, 2017, 7:11:57 AM (6 years ago)
Author:
Peter
Message:

merge release 0.14.1 into trunk

move functions 'create_gsl_vector' and 'create_gsl_vector_copy' from
Vector private to free functions in namespace detail, i.e., they are
available for other classes/functions, but not part of the public
API. In the same region add functions 'overlap' and 'serial_overlap'
to check if two GSL vectors overlap or can be used safely as LHS and
RHS in a function.

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/NEWS

    r3613 r3623  
    1010
    1111yat 0.14.x series from http://dev.thep.lu.se/yat/svn/branches/0.14-stable
     12
     13version 0.14.1 (released 8 February 2017)
     14  - Fixed bug in element-by-element operations in VectorMutable such
     15    as operator+=, operator-=, ::mul(), and ::div(), when underlying
     16    data structures of lhs and rhs overlap such that elements in rhs
     17    are affected by the operation on lhs before they have been used on
     18    rhs (bug #884).
     19
     20  A complete list of closed tickets can be found here [[br]]
     21  http://dev.thep.lu.se/yat/query?status=closed&milestone=yat+0.14.1
    1222
    1323version 0.14 (released 16 January 2017)
  • trunk/m4/version.m4

    r3566 r3623  
    8282# yat-0.13.1 10:1:0
    8383# yat-0.14   11:0:0
     84# yat-0.14.1 11:1:0
    8485#
    8586# *Accidently, the libtool number was not updated for yat 0.5
  • trunk/test/vector.cc

    r3612 r3623  
    77  Copyright (C) 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    88  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
    9   Copyright (C) 2012, 2016 Peter Johansson
     9  Copyright (C) 2012, 2016, 2017 Peter Johansson
    1010
    1111  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    433433}
    434434
     435
    435436void test_aliasing(test::Suite& suite)
    436437{
     438  suite.out() << "test_aliasing\n";
    437439  utility::Vector vec(10, 1);
    438440  utility::VectorView view(vec, 0, 5, 2);
    439441  utility::VectorConstView const_view(vec, 0, 5, 1);
     442  // vec =     { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
     443  // view =    { 1,    1,    1,    1,    1}
     444  // const_view{ 1, 1, 1, 1, 1}
    440445  view += const_view;
    441   suite.out() << "result: " << view << "\n";
     446  suite.out() << "vec: " << vec << "\n";
     447  suite.out() << "view: " << view << "\n";
     448  suite.out() << "const_view: " << const_view << "\n";
    442449  for (size_t i=0; i<view.size(); ++i) {
    443     if (i==2 || i==4)
    444       suite.xadd(suite.equal(view(i), 2));
    445     else
    446       suite.add(suite.equal(view(i), 2));
    447   }
    448 }
     450    suite.add(suite.equal(view(i), 2));
     451  }
     452
     453  // vec =     { 2, 1, 2, 1, 2, 1, 2, 1, 2, 1}
     454  // view =    { 2,    2,    2,    2,    2}
     455  // const_view{ 2, 1, 2, 1, 2}
     456  suite.out() << "vec: " << vec << "\n";
     457  suite.out() << "view: " << view << "\n";
     458  suite.out() << "const_view: " << const_view << "\n";
     459
     460  view = const_view;
     461  // view =    { 2,    1,    2,    1,    2}
     462  utility::Vector correct(5);
     463  correct(0) = 2;
     464  correct(1) = 1;
     465  correct(2) = 2;
     466  correct(3) = 1;
     467  correct(4) = 2;
     468  if (view != correct) {
     469    suite.add(false);
     470    suite.err() << "error: view: \n" << view << "\nexpected: " << correct
     471                << "\n";
     472  }
     473}
     474
  • trunk/yat/utility/Vector.cc

    r3605 r3623  
    5050
    5151  Vector::Vector(size_t n, double init_value)
    52     : VectorMutable(create_gsl_vector(n))
    53   {
    54     if (n)
    55       all(init_value);
     52    : VectorMutable(detail::create_gsl_vector(n, init_value))
     53  {
    5654    assert(n==0 || vec_);
    5755    assert(n==0 || const_vec_);
     
    122120
    123121    // convert the data to a gsl vector
    124     vec_ = create_gsl_vector(nof_rows*nof_columns);
     122    vec_ = detail::create_gsl_vector(nof_rows*nof_columns);
    125123    size_t n=0;
    126124    // if gsl error handler disabled, out of bounds index will not
     
    148146
    149147    if (size() != other->size) {
    150       gsl_vector* tmp = create_gsl_vector_copy(other);
    151       // delete memory once we know create copy did not fail
     148      gsl_vector* tmp = detail::create_gsl_vector_copy(other);
     149      // delete memory once we know created copy did not fail
    152150      delete_allocated_memory();
    153151      vec_ = tmp;
     
    174172  gsl_vector* Vector::create_gsl_vector_copy(const VectorBase& other) const
    175173  {
    176     return create_gsl_vector_copy(other.gsl_vector_p());
    177   }
    178 
    179 
    180   gsl_vector* Vector::create_gsl_vector_copy(const gsl_vector* other) const
    181   {
    182     if (!other)
    183       return NULL;
    184     gsl_vector* vec = create_gsl_vector(other->size);
    185     if (gsl_vector_memcpy(vec, other))
    186       throw utility::GSL_error("Vector::create_gsl_vector_copy memcpy failed");
    187     return vec;
    188   }
    189 
    190 
    191   gsl_vector* Vector::create_gsl_vector(size_t n) const
    192   {
    193     if (!n)
    194       return NULL;
    195     gsl_vector* vec = gsl_vector_alloc(n);
    196     if (!vec)
    197       throw utility::GSL_error("Vector: failed to allocate memory");
    198     return vec;
     174    return detail::create_gsl_vector_copy(other.gsl_vector_p());
    199175  }
    200176
     
    216192  void Vector::resize(size_t n, double init_value)
    217193  {
     194    if (size() == n) {
     195      all(init_value);
     196      return;
     197    }
     198    if (!n) {
     199      delete_allocated_memory();
     200      return;
     201    }
     202    gsl_vector* tmp = detail::create_gsl_vector(n, init_value);
    218203    delete_allocated_memory();
    219     if (!n)
    220       return;
    221     const_vec_ = vec_ = create_gsl_vector(n);
    222     all(init_value);
     204    const_vec_ = vec_ = tmp;
    223205  }
    224206
  • trunk/yat/utility/Vector.h

    r3605 r3623  
    8686    template<class T>
    8787    Vector(const VectorExpression<T>& other)
    88       : VectorMutable(create_gsl_vector_copy(other.gsl_vector_p()))
     88      : VectorMutable(detail::create_gsl_vector_copy(other.gsl_vector_p()))
    8989    {
    9090    }
     
    176176    {
    177177      // access rhs before deleting vec_
    178       gsl_vector* tmp = create_gsl_vector_copy(rhs.gsl_vector_p());
     178      gsl_vector* tmp = detail::create_gsl_vector_copy(rhs.gsl_vector_p());
    179179      delete_allocated_memory();
    180180      const_vec_ = vec_ = tmp;
     
    239239    gsl_vector* create_gsl_vector_copy(const VectorBase&) const;
    240240
    241     /**
    242        Same as above but takes the internal GSL vector as argument
    243      */
    244     gsl_vector* create_gsl_vector_copy(const gsl_vector*) const;
    245 
    246     /**
    247        \brief Create a new gsl vector
    248 
    249        Necessary memory for the new GSL vector is allocated and the
    250        caller is responsible for freeing the allocated memory.
    251 
    252        \return A pointer to created GSL vector.
    253 
    254        \throw GSL_error if memory cannot be allocated for the new
    255        GSL vector.
    256     */
    257     gsl_vector* create_gsl_vector(size_t n) const;
    258241
    259242    void delete_allocated_memory(void);
  • trunk/yat/utility/VectorBase.cc

    r3550 r3623  
    226226  }
    227227
     228  namespace detail
     229  {
     230    gsl_vector* create_gsl_vector_copy(const gsl_vector* other)
     231    {
     232      if (!other)
     233        return NULL;
     234      gsl_vector* vec = create_gsl_vector(other->size);
     235      if (gsl_vector_memcpy(vec, other))
     236        throw utility::GSL_error("create_gsl_vector_copy memcpy failed");
     237      return vec;
     238    }
     239
     240
     241    gsl_vector* create_gsl_vector(size_t n)
     242    {
     243      if (!n)
     244        return NULL;
     245      gsl_vector* vec = gsl_vector_alloc(n);
     246      if (!vec)
     247        throw utility::GSL_error("gsl_vector_alloc failed");
     248      return vec;
     249    }
     250
     251
     252    gsl_vector* create_gsl_vector(size_t n, double value)
     253    {
     254      if (!n)
     255        return NULL;
     256      if (value) {
     257        gsl_vector* vec = create_gsl_vector(n);
     258        gsl_vector_set_all(vec, value);
     259        return vec;
     260      }
     261      gsl_vector* vec = gsl_vector_calloc(n);
     262      if (!vec)
     263        throw utility::GSL_error("gsl_vector_calloc failed");
     264      return vec;
     265    }
     266
     267
     268    bool overlap(const gsl_vector* a, const gsl_vector* b)
     269    {
     270      assert(a);
     271      assert(b);
     272      assert(a->size);
     273      assert(b->size);
     274      const double* a_first = gsl_vector_const_ptr(a, 0);
     275      const double* a_last = gsl_vector_const_ptr(a, a->size-1);
     276      const double* b_first = gsl_vector_const_ptr(b, 0);
     277      const double* b_last = gsl_vector_const_ptr(b, b->size-1);
     278
     279      return (a_last >= b_first) && (b_last >= a_first);
     280    }
     281
     282
     283    bool serial_overlap(const gsl_vector* a, const gsl_vector* b)
     284    {
     285      assert(a);
     286      assert(b);
     287      assert(a->size);
     288      assert(b->size);
     289      assert(a->size == b->size);
     290      /*
     291        ABCDEFGHIJKLMNOPQRSTUVWXYZ
     292        0ab
     293        1  b a
     294        2   b    a
     295        3    b       a
     296        4     b          a
     297        5      b             a
     298        6       b                a
     299
     300        a_first in mem pos A
     301        a_last in mem pos Y
     302        b_first in mem pos B
     303        b_last in mem pos H
     304
     305        Both a[1] and b[3] point to mem pos E
     306
     307        From this follows that a mem pos is problematic if it
     308        fullfils the following:
     309        1) It lies inside [a_first, a_last]
     310        2) It lies inside [b_first, b_last]
     311
     312        3) The b-line is below the a-line for that position, which is
     313        a linear condition, hence no intermediate extreme point, and
     314        we only need to look at the boundaries.
     315      */
     316
     317      const double* a_first = gsl_vector_const_ptr(a, 0);
     318      const double* a_last = gsl_vector_const_ptr(a, a->size-1);
     319      const double* b_first = gsl_vector_const_ptr(b, 0);
     320      const double* b_last = gsl_vector_const_ptr(b, b->size-1);
     321      if ((a_last >= b_first) && (b_last >= a_first)) {
     322        const double* first = std::max(a_first, b_first);
     323        const double* last = std::min(a_last, b_last);
     324        assert(first <= last);
     325
     326        // check whether b-line is below a-line in first or last,
     327        // which is equivalent index(b) > index(a)
     328
     329        // index in x for a is calculated as:
     330        // (x - a_first) / a->stride
     331        // so
     332        // (x - b_first) * a->stride > (x - a_first) * b->stride)
     333        return ((first - b_first) * a->stride >
     334                (first - a_first) * b->stride) ||
     335          ((last - b_first) * a->stride >
     336           (last - a_first) * b->stride);
     337      }
     338      return false;
     339    }
     340
     341  }
     342
    228343}}} // of namespace utility, yat, and thep
  • trunk/yat/utility/VectorBase.h

    r3605 r3623  
    296296  std::ostream& operator<<(std::ostream& s, const VectorBase& v);
    297297
     298  /// \cond IGNORE_DOXYGEN
     299
     300  // Some convenience functions used in Vector and friends.
     301  namespace detail {
     302
     303    /**
     304       \brief Create a new gsl vector
     305
     306       Necessary memory for the new GSL vector is allocated and the
     307       caller is responsible for freeing the allocated memory.
     308
     309       \return A pointer to created GSL vector.
     310
     311       \throw GSL_error if memory cannot be allocated for the new
     312       GSL vector.
     313    */
     314    gsl_vector* create_gsl_vector(size_t n);
     315
     316    /**
     317       same as above but sets all values to init
     318     */
     319    gsl_vector* create_gsl_vector(size_t n, double init);
     320
     321    /**
     322       \brief Create a new copy of the internal GSL vector.
     323
     324       Necessary memory for the new GSL vector is allocated and the
     325       caller is responsible for freeing the allocated memory.
     326
     327       \return A pointer to a copy of the internal GSL vector.
     328
     329       \throw GSL_error if memory cannot be allocated for the new
     330       copy, or if dimensions mis-match.
     331    */
     332    gsl_vector* create_gsl_vector_copy(const gsl_vector*);
     333
     334
     335    // return true if a and b overlap in memory
     336    bool overlap(const gsl_vector* a, const gsl_vector* b);
     337
     338    /*
     339      Is a stricter condition that overlap(2). Basically it returns
     340      false if it is safe to write code like:
     341      for (size_t i=0; i<a->size; ++i)
     342        a[i] = foo(b[i])
     343
     344      in other words, it returns true if there exists i and j, i<j, such that
     345      gsl_vector_const_ptr(a, j) == gsl_vector_const_ptr(b, i)
     346
     347      For speed reasons function might return true in some cases when
     348      condition above is not true, when both a and b have stride>1 and
     349      elements sit between each other.
     350
     351      Note, a and b must have same size.
     352     */
     353    bool serial_overlap(const gsl_vector* a, const gsl_vector* b);
     354  }
     355
     356  /// \endcond
     357
    298358}}} // of namespace utility, yat, and theplu
    299359
  • trunk/yat/utility/VectorMutable.cc

    r3550 r3623  
    66  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
    8   Copyright (C) 2010, 2012, 2016 Peter Johansson
     8  Copyright (C) 2010, 2012, 2016, 2017 Peter Johansson
    99
    1010  This file is part of the yat library, http://dev.thep.lu.se/trac/yat
     
    2727
    2828#include "VectorMutable.h"
     29#include "Vector.h"
    2930#include "Exception.h"
    3031#include "stl_utility.h"
     
    8485  {
    8586    assert(vec_);
     87    if (detail::serial_overlap(const_vec_, other.gsl_vector_p())) {
     88      Vector tmp(other);
     89      div(tmp);
     90      return;
     91    }
    8692    int status=gsl_vector_div(vec_,other.gsl_vector_p());
    8793    if (status)
     
    108114  {
    109115    assert(vec_);
     116    if (detail::serial_overlap(const_vec_, other.gsl_vector_p())) {
     117      Vector tmp(other);
     118      mul(tmp);
     119      return;
     120    }
    110121    int status=gsl_vector_mul(vec_,other.gsl_vector_p());
    111122    if (status)
     
    128139      throw runtime_error("VectorMutable::operator+= size must be same.");
    129140    assert(vec_);
     141
     142    // if ranges overlap be a bit careful
     143    if (detail::serial_overlap(const_vec_, other.gsl_vector_p())) {
     144      Vector tmp(other);
     145      *this += tmp;
     146      return *this;
     147    }
    130148    int status=gsl_vector_add(vec_, other.gsl_vector_p());
    131149    if (status)
     
    148166      throw runtime_error("VectorMutable::operator-= size must be same.");
    149167    assert(vec_);
     168    if (detail::serial_overlap(const_vec_, other.gsl_vector_p())) {
     169      Vector tmp(other);
     170      *this -= tmp;
     171      return *this;
     172    }
    150173    int status=gsl_vector_sub(vec_, other.gsl_vector_p());
    151174    if (status)
  • trunk/yat/utility/VectorView.cc

    r3605 r3623  
    66  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
    8   Copyright (C) 2012, 2013 Peter Johansson
     8  Copyright (C) 2012, 2013, 2017 Peter Johansson
    99
    1010  This file is part of the yat library, http://dev.thep.lu.se/trac/yat
     
    2828#include "VectorView.h"
    2929#include "VectorMutable.h"
     30#include "Vector.h"
    3031#include "Matrix.h"
    3132#include "utility.h"
     
    101102    if (size() != other->size)
    102103      throw utility::GSL_error("VectorView::assign dimension mis-match");
    103     if (gsl_vector_memcpy(vec_, other))
     104
     105    // indirect assignment
     106    if (detail::serial_overlap(const_vec_, other)) {
     107      gsl_vector* tmp = detail::create_gsl_vector_copy(other);
     108      if (gsl_vector_memcpy(vec_, tmp)) {
     109        gsl_vector_free(tmp);
     110        throw utility::GSL_error("VectorView::assign memcpy failed.");
     111      }
     112      gsl_vector_free(tmp);
     113    }
     114    else if (gsl_vector_memcpy(vec_, other))
    104115      throw utility::GSL_error("VectorView::assign memcpy failed.");
    105116    const_vec_ = vec_;
Note: See TracChangeset for help on using the changeset viewer.