Quelle AmbiVector.h Sprache: C

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_AMBIVECTOR_H
#define EIGEN_AMBIVECTOR_H

namespace Eigen {

namespace internal {

/** \internal
  * Hybrid sparse/dense vector class designed for intensive read-write operations.
  *
  * See BasicSparseLLT and SparseProduct for usage examples.
  */
template<typename _Scalar, typename _StorageIndex>
class AmbiVector
{
  public:
    typedef _Scalar Scalar;
    typedef _StorageIndex StorageIndex;
    typedef typename NumTraits<Scalar>::Real RealScalar;

    explicit AmbiVector(Index size)
      : m_buffer(0), m_zero(0), m_size(0), m_end(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1)
    {
      resize(size);
    }

    void init(double estimatedDensity);
    void init(int mode);

    Index nonZeros() const;

    /** Specifies a sub-vector to work on */
    void setBounds(Index start, Index end) { m_start = convert_index(start); m_end = convert_index(end); }

    void setZero();

    void restart();
    Scalar& coeffRef(Index i);
    Scalar& coeff(Index i);

    class Iterator;

    ~AmbiVector() { delete[] m_buffer; }

    void resize(Index size)
    {
      if (m_allocatedSize < size)
        reallocate(size);
      m_size = convert_index(size);
    }

    StorageIndex size() const { return m_size; }

  protected:
    StorageIndex convert_index(Index idx)
    {
      return internal::convert_index<StorageIndex>(idx);
    }

    void reallocate(Index size)
    {
      // if the size of the matrix is not too large, let's allocate a bit more than needed such
      // that we can handle dense vector even in sparse mode.
      delete[] m_buffer;
      if (size<1000)
      {
        Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar);
        m_allocatedElements = convert_index((allocSize*sizeof(Scalar))/sizeof(ListEl));
        m_buffer = new Scalar[allocSize];
      }
      else
      {
        m_allocatedElements = convert_index((size*sizeof(Scalar))/sizeof(ListEl));
        m_buffer = new Scalar[size];
      }
      m_size = convert_index(size);
      m_start = 0;
      m_end = m_size;
    }

    void reallocateSparse()
    {
      Index copyElements = m_allocatedElements;
      m_allocatedElements = (std::min)(StorageIndex(m_allocatedElements*1.5),m_size);
      Index allocSize = m_allocatedElements * sizeof(ListEl);
      allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar);
      Scalar* newBuffer = new Scalar[allocSize];
      std::memcpy(newBuffer,  m_buffer,  copyElements * sizeof(ListEl));
      delete[] m_buffer;
      m_buffer = newBuffer;
    }

  protected:
    // element type of the linked list
    struct ListEl
    {
      StorageIndex next;
      StorageIndex index;
      Scalar value;
    };

    // used to store data in both mode
    Scalar* m_buffer;
    Scalar m_zero;
    StorageIndex m_size;
    StorageIndex m_start;
    StorageIndex m_end;
    StorageIndex m_allocatedSize;
    StorageIndex m_allocatedElements;
    StorageIndex m_mode;

    // linked list mode
    StorageIndex m_llStart;
    StorageIndex m_llCurrent;
    StorageIndex m_llSize;
};

/** \returns the number of non zeros in the current sub vector */
template<typename _Scalar,typename _StorageIndex>
Index AmbiVector<_Scalar,_StorageIndex>::nonZeros() const
{
  if (m_mode==IsSparse)
    return m_llSize;
  else
    return m_end - m_start;
}

template<typename _Scalar,typename _StorageIndex>
void AmbiVector<_Scalar,_StorageIndex>::init(double estimatedDensity)
{
  if (estimatedDensity>0.1)
    init(IsDense);
  else
    init(IsSparse);
}

template<typename _Scalar,typename _StorageIndex>
void AmbiVector<_Scalar,_StorageIndex>::init(int mode)
{
  m_mode = mode;
  // This is only necessary in sparse mode, but we set these unconditionally to avoid some maybe-uninitialized warnings
  // if (m_mode==IsSparse)
  {
    m_llSize = 0;
    m_llStart = -1;
  }
}

/** Must be called whenever we might perform a write access
  * with an index smaller than the previous one.
  *
  * Don't worry, this function is extremely cheap.
  */
template<typename _Scalar,typename _StorageIndex>
void AmbiVector<_Scalar,_StorageIndex>::restart()
{
  m_llCurrent = m_llStart;
}

/** Set all coefficients of current subvector to zero */
template<typename _Scalar,typename _StorageIndex>
void AmbiVector<_Scalar,_StorageIndex>::setZero()
{
  if (m_mode==IsDense)
  {
    for (Index i=m_start; i<m_end; ++i)
      m_buffer[i] = Scalar(0);
  }
  else
  {
    eigen_assert(m_mode==IsSparse);
    m_llSize = 0;
    m_llStart = -1;
  }
}

template<typename _Scalar,typename _StorageIndex>
_Scalar& AmbiVector<_Scalar,_StorageIndex>::coeffRef(Index i)
{
  if (m_mode==IsDense)
    return m_buffer[i];
  else
  {
    ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer);
    // TODO factorize the following code to reduce code generation
    eigen_assert(m_mode==IsSparse);
    if (m_llSize==0)
    {
      // this is the first element
      m_llStart = 0;
      m_llCurrent = 0;
      ++m_llSize;
      llElements[0].value = Scalar(0);
      llElements[0].index = convert_index(i);
      llElements[0].next = -1;
      return llElements[0].value;
    }
    else if (i<llElements[m_llStart].index)
    {
      // this is going to be the new first element of the list
      ListEl& el = llElements[m_llSize];
      el.value = Scalar(0);
      el.index = convert_index(i);
      el.next = m_llStart;
      m_llStart = m_llSize;
      ++m_llSize;
      m_llCurrent = m_llStart;
      return el.value;
    }
    else
    {
      StorageIndex nextel = llElements[m_llCurrent].next;
      eigen_assert(i>=llElements[m_llCurrent].index && "you must call restart() before inserting an element with lower or equal index");
      while (nextel >= 0 && llElements[nextel].index<=i)
      {
        m_llCurrent = nextel;
        nextel = llElements[nextel].next;
      }

      if (llElements[m_llCurrent].index==i)
      {
        // the coefficient already exists and we found it !
        return llElements[m_llCurrent].value;
      }
      else
      {
        if (m_llSize>=m_allocatedElements)
        {
          reallocateSparse();
          llElements = reinterpret_cast<ListEl*>(m_buffer);
        }
        eigen_internal_assert(m_llSize<m_allocatedElements && "internal error: overflow in sparse mode");
        // let's insert a new coefficient
        ListEl& el = llElements[m_llSize];
        el.value = Scalar(0);
        el.index = convert_index(i);
        el.next = llElements[m_llCurrent].next;
        llElements[m_llCurrent].next = m_llSize;
        ++m_llSize;
        return el.value;
      }
    }
  }
}

template<typename _Scalar,typename _StorageIndex>
_Scalar& AmbiVector<_Scalar,_StorageIndex>::coeff(Index i)
{
  if (m_mode==IsDense)
    return m_buffer[i];
  else
  {
    ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer);
    eigen_assert(m_mode==IsSparse);
    if ((m_llSize==0) || (i<llElements[m_llStart].index))
    {
      return m_zero;
    }
    else
    {
      Index elid = m_llStart;
      while (elid >= 0 && llElements[elid].index<i)
        elid = llElements[elid].next;

      if (llElements[elid].index==i)
        return llElements[m_llCurrent].value;
      else
        return m_zero;
    }
  }
}

/** Iterator over the nonzero coefficients */
template<typename _Scalar,typename _StorageIndex>
class AmbiVector<_Scalar,_StorageIndex>::Iterator
{
  public:
    typedef _Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;

    /** Default constructor
      * \param vec the vector on which we iterate
      * \param epsilon the minimal value used to prune zero coefficients.
      * In practice, all coefficients having a magnitude smaller than \a epsilon
      * are skipped.
      */
    explicit Iterator(const AmbiVector& vec, const RealScalar& epsilon = 0)
      : m_vector(vec)
    {
      using std::abs;
      m_epsilon = epsilon;
      m_isDense = m_vector.m_mode==IsDense;
      if (m_isDense)
      {
        m_currentEl = 0;   // this is to avoid a compilation warning
        m_cachedValue = 0; // this is to avoid a compilation warning
        m_cachedIndex = m_vector.m_start-1;
        ++(*this);
      }
      else
      {
        ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
        m_currentEl = m_vector.m_llStart;
        while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon)
          m_currentEl = llElements[m_currentEl].next;
        if (m_currentEl<0)
        {
          m_cachedValue = 0; // this is to avoid a compilation warning
          m_cachedIndex = -1;
        }
        else
        {
          m_cachedIndex = llElements[m_currentEl].index;
          m_cachedValue = llElements[m_currentEl].value;
        }
      }
    }

    StorageIndex index() const { return m_cachedIndex; }
    Scalar value() const { return m_cachedValue; }

    operator bool() const { return m_cachedIndex>=0; }

    Iterator& operator++()
    {
      using std::abs;
      if (m_isDense)
      {
        do {
          ++m_cachedIndex;
        } while (m_cachedIndex<m_vector.m_end && abs(m_vector.m_buffer[m_cachedIndex])<=m_epsilon);
        if (m_cachedIndex<m_vector.m_end)
          m_cachedValue = m_vector.m_buffer[m_cachedIndex];
        else
          m_cachedIndex=-1;
      }
      else
      {
        ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
        do {
          m_currentEl = llElements[m_currentEl].next;
        } while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon);
        if (m_currentEl<0)
        {
          m_cachedIndex = -1;
        }
        else
        {
          m_cachedIndex = llElements[m_currentEl].index;
          m_cachedValue = llElements[m_currentEl].value;
        }
      }
      return *this;
    }

  protected:
    const AmbiVector& m_vector; // the target vector
    StorageIndex m_currentEl;   // the current element in sparse/linked-list mode
    RealScalar m_epsilon;       // epsilon used to prune zero coefficients
    StorageIndex m_cachedIndex; // current coordinate
    Scalar m_cachedValue;       // current value
    bool m_isDense;             // mode of the vector
};

} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_AMBIVECTOR_H

quality90%

¤ Dauer der Verarbeitung: 0.14 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.