sparse_vector_blocked.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/type_traits.hpp>
#include <kernel/lafem/container.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/util/tiny_algebra.hpp>
#include <kernel/util/math.hpp>
#include <kernel/lafem/arch/max_abs_index.hpp>
#include <kernel/lafem/arch/min_abs_index.hpp>
#include <kernel/lafem/arch/max_index.hpp>
#include <kernel/lafem/arch/min_index.hpp>
#include <kernel/lafem/arch/max_rel_diff.hpp>
#include <kernel/adjacency/permutation.hpp>
#include <array>
 
namespace FEAT
{
  namespace LAFEM
  {
    namespace Intern
    {
      template<typename DT_, int BlockSize_, Perspective perspective_>
      struct SparseVectorBlockedPerspectiveHelper
      {
        typedef Tiny::Vector<DT_, BlockSize_> Type;
      };
 
      template<typename DT_, int BlockSize_>
      struct SparseVectorBlockedPerspectiveHelper<DT_, BlockSize_, Perspective::pod>
      {
        typedef DT_ Type;
      };
    } // namespace Intern
 
    template <typename DT_, typename IT_, int BlockSize_>
    class SparseVectorBlocked : public Container<DT_, IT_>
    {
    private:
      template <typename T1_, typename T2_>
      static void _insertion_sort(T1_ * key, T2_ * val1, Index size)
      {
        T1_ swap_key;
        T2_ swap1;
        for (Index i(1), j ; i < size ; ++i)
        {
          swap_key = key[i];
          swap1 = val1[i];
          j = i;
          while (j > 0 && key[j-1] > swap_key)
          {
            key[j] = key[j-1];
            val1[j] = val1[j-1];
            --j;
          }
          key[j] = swap_key;
          val1[j] = swap1;
        }
      }
 
      Index & _used_elements()
      {
        return this->_scalar_index.at(1);
      }
 
      Index & _allocated_elements()
      {
        return this->_scalar_index.at(2);
      }
 
      Index & _sorted()
      {
        return this->_scalar_index.at(4);
      }
 
      bool _remove_element(IT_ ind)
      {
        IT_* pindices = this->_indices.at(0);
        IT_* ptr = std::find(pindices, pindices + _used_elements(), ind);
        if(ptr == pindices + _used_elements())
        {
          return false;
        }
        _sorted() = 0;
        *ptr = std::numeric_limits<IT_>::max();
        sort();
        return true;
      }
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      static constexpr int BlockSize = BlockSize_;
      typedef Tiny::Vector<DT_, BlockSize_> ValueType;
 
      explicit SparseVectorBlocked() :
        Container<DT_, IT_> (0)
      {
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(Math::min<Index>(0, 1000));
        this->_scalar_index.push_back(1);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      explicit SparseVectorBlocked(std::vector<char> input) :
         Container<DT_, IT_>(0)
      {
        deserialize<DT2_,IT2_>(input);
      }
 
      explicit SparseVectorBlocked(FileMode mode, const String& filename) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, filename);
      }
 
      explicit SparseVectorBlocked(FileMode mode, std::istream& file) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, file);
      }
 
 
      explicit SparseVectorBlocked(Index size_in) :
        Container<DT_, IT_>(size_in)
      {
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(Math::min<Index>(size_in, 1000));
        this->_scalar_index.push_back(1);
      }
 
      explicit SparseVectorBlocked(Index size_in, DenseVectorBlocked<DT_, IT_, BlockSize_> & elements_in,
                            DenseVector<IT_, IT_> & indices_in, bool is_sorted = true) :
        Container<DT_, IT_>(size_in)
      {
        XASSERT(size_in != Index(0));
        XASSERTM(indices_in.size() == elements_in.size(), "Vector size mismatch!");
 
        this->_scalar_index.push_back(elements_in.size());
        this->_scalar_index.push_back(elements_in.size());
        this->_scalar_index.push_back(Math::min<Index>(size_in, 1000));
        this->_scalar_index.push_back(Index(is_sorted));
 
        this->_elements.push_back(elements_in.template elements<Perspective::pod>());
        this->_elements_size.push_back(elements_in.template size<Perspective::pod>());
        this->_indices.push_back(indices_in.elements());
        this->_indices_size.push_back(indices_in.size());
 
        for (Index i(0) ; i < this->_elements.size() ; ++i)
          MemoryPool::increase_memory(this->_elements.at(i));
        for (Index i(0) ; i < this->_indices.size() ; ++i)
          MemoryPool::increase_memory(this->_indices.at(i));
 
        this->sort();
      }
 
      SparseVectorBlocked(SparseVectorBlocked && other) :
        Container<DT_, IT_>(std::forward<SparseVectorBlocked>(other))
      {
      }
 
      SparseVectorBlocked & operator= (SparseVectorBlocked && other)
      {
        this->move(std::forward<SparseVectorBlocked>(other));
 
        return *this;
      }
 
      SparseVectorBlocked clone(CloneMode clone_mode = CloneMode::Deep) const
      {
        SparseVectorBlocked t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template<typename DT2_, typename IT2_>
      void clone(const SparseVectorBlocked<DT2_, IT2_, BlockSize_> & other, CloneMode clone_mode = CloneMode::Deep)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseVectorBlocked<DT2_, IT2_, BlockSize_> & other)
      {
        this->sort();
        this->clone(other);
      }
 
      template <Perspective perspective_ = Perspective::native>
      auto elements() const -> const typename Intern::SparseVectorBlockedPerspectiveHelper<DT_, BlockSize_, perspective_>::Type *
      {
        if (this->empty())
          return nullptr;
 
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
        return (const typename Intern::SparseVectorBlockedPerspectiveHelper<DT_, BlockSize_, perspective_>::Type *)(this->_elements.at(0));
      }
 
      template <Perspective perspective_ = Perspective::native>
      auto elements() -> typename Intern::SparseVectorBlockedPerspectiveHelper<DT_, BlockSize_, perspective_>::Type *
      {
        if (this->empty())
          return nullptr;
 
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
        return (typename Intern::SparseVectorBlockedPerspectiveHelper<DT_, BlockSize_, perspective_>::Type *)(this->_elements.at(0));
      }
 
      IT_ * indices()
      {
        if(this->_indices.empty())
          return nullptr;
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
        return this->_indices.at(0);
      }
 
      IT_ const * indices() const
      {
        if(this->_indices.empty())
          return nullptr;
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
        return this->_indices.at(0);
      }
 
    template <Perspective perspective_ = Perspective::native>
      Index size() const
      {
        if (perspective_ == Perspective::pod)
          return static_cast<const Container<DT_, IT_> *>(this)->size() * Index(BlockSize_);
        else
          return static_cast<const Container<DT_, IT_> *>(this)->size();
      }
 
      const ValueType operator()(Index index) const
      {
        ASSERTM(index < this->_scalar_index.at(0), "index exceeds sparse vector size");
 
        MemoryPool::synchronize();
 
        if (this->_elements.empty())
          return ValueType(0.);
 
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
 
        Index i(0);
        while (i < used_elements())
        {
          if (indices()[i] >= index)
            break;
          ++i;
        }
 
        if (i < used_elements() && indices()[i] == index)
        {
          return this->elements()[i];
        }
        else
          return ValueType(0.);
      }
 
 
      void operator()(Index index, const ValueType& val)
      {
        ASSERTM(index < this->_scalar_index.at(0), "index exceeds sparse vector size");
 
        // flag container as not sorted anymore
        // CAUTION: do not use any method triggering resorting until we are finished
        _sorted() = 0;
 
        // vector is empty, no arrays allocated
        if (this->_elements.size() == 0)
        {
          this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(alloc_increment() * Index(BlockSize_)));
          this->_elements_size.push_back(alloc_increment() * Index(BlockSize_));
          MemoryPool::template set_memory<DT_>(this->_elements.back(), DT_(4711), alloc_increment() * Index(BlockSize_));
          this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(alloc_increment()));
          this->_indices_size.push_back(alloc_increment());
          MemoryPool::template set_memory<IT_>(this->_indices.back(), IT_(4711), alloc_increment());
          _allocated_elements() = alloc_increment();
          MemoryPool::copy(this->_elements.at(0), val.v, Index(BlockSize_));
          MemoryPool::set_memory(this->_indices.at(0), IT_(index));
          _used_elements() = 1;
        }
 
        // append element in already allocated arrays
        else if(_used_elements() < allocated_elements())
        {
          MemoryPool::copy(this->_elements.at(0) + _used_elements() * Index(BlockSize_), val.v,
          Index(BlockSize_));
          MemoryPool::set_memory(this->_indices.at(0) + _used_elements(), IT_(index));
          ++_used_elements();
        }
 
        // reallocate arrays, append element
        else
        {
          _allocated_elements() += alloc_increment();
 
          DT_ * elements_new(MemoryPool::template allocate_memory<DT_>(
            allocated_elements() * Index(BlockSize_)));
          MemoryPool::template set_memory<DT_>(elements_new, DT_(4711),
          allocated_elements() * Index(BlockSize_));
          IT_ * indices_new(MemoryPool::template allocate_memory<IT_>(allocated_elements()));
          MemoryPool::template set_memory<IT_>(indices_new, IT_(4711), allocated_elements());
 
          MemoryPool::copy(elements_new, this->_elements.at(0), _used_elements() * Index(BlockSize_));
          MemoryPool::copy(indices_new, this->_indices.at(0), _used_elements());
 
          MemoryPool::release_memory(this->_elements.at(0));
          MemoryPool::release_memory(this->_indices.at(0));
 
          this->_elements.at(0) = elements_new;
          this->_indices.at(0) = indices_new;
 
          MemoryPool::copy(this->_elements.at(0) + used_elements() * Index(BlockSize_), val.v,
          Index(BlockSize_));
          MemoryPool::set_memory(this->_indices.at(0) + _used_elements(), IT_(index));
 
          ++_used_elements();
          this->_elements_size.at(0) = allocated_elements() * Index(BlockSize_);
          this->_indices_size.at(0) = allocated_elements();
        }
      }
 
      void sort()
      {
        if (sorted() == 0)
        {
          //first of all, mark vector as sorted, because otherwise we would call ourselves inifite times
          // CAUTION: do not use any method triggering resorting until we are finished
          _sorted() = 1;
 
          // check if there is anything to be sorted
          if(_used_elements() == Index(0))
            return;
 
          IT_ * pindices;
          ValueType * pelements;
          pindices = this->_indices.at(0);
          pelements = this->elements();
 
          _insertion_sort(pindices, pelements, _used_elements());
 
          // find and mark duplicate entries
          for (Index i(1) ; i < _used_elements() ; ++i)
          {
            if (pindices[i-1] == pindices[i])
            {
              pindices[i-1] = std::numeric_limits<IT_>::max();
            }
          }
 
          // sort out marked duplicated elements
          _insertion_sort(pindices, pelements, _used_elements());
          Index junk(0);
          while (pindices[_used_elements() - 1 - junk] == std::numeric_limits<IT_>::max() && junk < _used_elements())
            ++junk;
          _used_elements() -= junk;
        }
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      void deserialize(std::vector<char> input)
      {
        this->template _deserialize<DT2_, IT2_>(FileMode::fm_svb, input);
      }
 
 
      template<typename IndexContainer>
      bool remove_elements(const IndexContainer& inds)
      {
        bool success = true;
        for(auto i : inds)
        {
          success &= _remove_element(IndexType(i));
        }
        return success;
      }
 
      bool remove_element(IndexType ind)
      {
        return  _remove_element(ind);
      }
 
      bool element_exists(IndexType ind) const
      {
        const IT_* pindices = this->_indices.at(0);
        return std::find(pindices, pindices + this->_scalar_index.at(1), ind) != pindices + this->_scalar_index.at(1);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      std::vector<char> serialize(const LAFEM::SerialConfig& config = SerialConfig()) const
      {
        return this->template _serialize<DT2_, IT2_>(FileMode::fm_svb, config);
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        std::ifstream file(filename.c_str(), std::ifstream::in | std::ifstream::binary);
        if (! file.is_open())
          XABORTM("Unable to open Vector file " + filename);
        read_from(mode, file);
        file.close();
      }
 
      void read_from(FileMode mode, std::istream& file)
      {
        switch(mode)
        {
          case FileMode::fm_binary:
          case FileMode::fm_svb:
            this->template _deserialize<double, std::uint64_t>(FileMode::fm_svb, file);
            break;
          default:
            XABORTM("Filemode not supported!");
        }
      }
 
      void write_out(FileMode mode, const String& filename) const
      {
        std::ios_base::openmode bin = std::ofstream::out | std::ofstream::binary;
        if(mode == FileMode::fm_mtx)
          bin = std::ofstream::out;
        std::ofstream file;
        char* buff = nullptr;
        if(mode == FileMode::fm_mtx)
        {
          buff = new char[LAFEM::FileOutStreamBufferSize];
          file.rdbuf()->pubsetbuf(buff, LAFEM::FileOutStreamBufferSize);
        }
        file.open(filename.c_str(), bin);
        if(! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
        write_out(mode, file);
        file.close();
        delete[] buff;
      }
 
      void write_out(FileMode mode, std::ostream& file) const
      {
        switch(mode)
        {
          case FileMode::fm_mtx:
          {
 
            file << "%%MatrixMarket matrix coordinate real general\n";
            file << this->size()*BlockSize << " " << 1 << " " << this->used_elements()*BlockSize << "\n";
 
            const Index u_elem(this->used_elements());
            const IT_ * pind(this->indices());
            const DT_ * pval(this->elements<Perspective::pod>());
            for (Index i(0) ; i < u_elem ; ++i, ++pind, pval+=BlockSize)
            {
              for(Index k(0); k < BlockSize; ++k)
                file << (*pind)*BlockSize+k+1 << " " << 1 << " " << stringify_fp_sci(*(pval+k)) << "\n";
            }
            break;
          }
 
          case FileMode::fm_binary:
          case FileMode::fm_svb:
            this->template _serialize<double, std::uint64_t>(FileMode::fm_svb, file);
            break;
          default:
            XABORTM("Filemode not supported!");
        }
      }
 
 
      DT_ max_abs_element() const
      {
        TimeStamp ts_start;
 
        Index max_abs_index = Arch::MaxAbsIndex::value(this->template elements<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(max_abs_index < this->template used_elements<Perspective::pod>());
        DT_ result(this->template elements<Perspective::pod>()[max_abs_index]);
        result = Math::abs(result);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ min_abs_element() const
      {
        TimeStamp ts_start;
 
        Index min_abs_index = Arch::MinAbsIndex::value(this->template elements<Perspective::pod>(), this->template  used_elements<Perspective::pod>());
        ASSERT(min_abs_index < this->template used_elements<Perspective::pod>());
        DT_ result(this->template elements<Perspective::pod>()[min_abs_index]);
        result = Math::abs(result);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_element() const
      {
        TimeStamp ts_start;
 
        Index max_index = Arch::MaxIndex::value(this->template elements<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(max_index < this->template used_elements<Perspective::pod>());
        DT_ result(this->template elements<Perspective::pod>()[max_index]);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ min_element() const
      {
        TimeStamp ts_start;
 
        Index min_index = Arch::MinIndex::value(this->template elements<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(min_index < this->template used_elements<Perspective::pod>());
        DT_ result(this->template elements<Perspective::pod>()[min_index]);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_rel_diff(const SparseVectorBlocked & x) const
      {
        XASSERTM(x.used_elements() == this->used_elements(), "Nonzero count does not match!");
        TimeStamp ts_start;
 
        DataType max_rel_diff = Arch::MaxRelDiff::value(this->template elements<Perspective::pod>(),
          x.template elements<Perspective::pod>(), this->template used_elements<Perspective::pod>());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return max_rel_diff;
      }
 
 
      void permute(Adjacency::Permutation & perm)
      {
        if (perm.size() == 0)
          return;
 
        XASSERTM(perm.size() == this->size(), "Container size does not match permutation size");
 
        SparseVectorBlocked<DT_, IT_, BlockSize_> target(this->size());
 
        auto inv = perm.inverse();
        const Index * const inv_pos(inv.get_perm_pos());
        for (Index i(0) ; i < this->used_elements() ; ++i)
        {
          const Index col = this->indices()[i];
          target(inv_pos[col], (*this)(col));
        }
 
        target.sort();
        this->assign(target);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        if (sorted() == 0)
          const_cast<SparseVectorBlocked *>(this)->sort();
        if (perspective_ == Perspective::pod)
          return this->_scalar_index.at(1) * Index(BlockSize_);
        else
          return this->_scalar_index.at(1);
      }
 
      Index allocated_elements() const
      {
        if(this->_scalar_index.empty())
          return Index(0);
        return this->_scalar_index.at(2);
      }
 
      Index alloc_increment() const
      {
        if(this->_scalar_index.empty())
          return Index(0);
        return this->_scalar_index.at(3);
      }
 
      Index sorted() const
      {
        if(this->_scalar_index.empty())
          return Index(0);
        return this->_scalar_index.at(4);
      }
 
      static String name()
      {
        return "SparseVectorBlocked";
      }
 
      bool same_layout(const SparseVectorBlocked& other) const
      {
        if(this->size() != other.size())
          return false;
        if(this->used_elements() != other.used_elements())
          return false;
 
        Index n = this->used_elements();
        if(n == Index(0))
          return true; // both vectors are empty
 
        const IT_* a = this->indices();
        const IT_* b = other.indices();
 
        // shallow copy? (aka same array)
        if(a == b)
          return true;
 
        // check all array entries
        for(Index i(0); i < n; ++i)
        {
          if(a[i] != b[i])
            return false;
        }
 
        // okay, arrays are identical
        return true;
      }
 
      friend std::ostream & operator<< (std::ostream & lhs, const SparseVectorBlocked & b)
      {
        lhs << "[";
        for (Index i(0) ; i < b.size() ; ++i)
        {
          ValueType t = b(i);
          for (int j(0) ; j < BlockSize_ ; ++j)
            lhs << "  " << stringify(t[j]);
        }
        lhs << "]";
 
        return lhs;
      }
    }; // class SparseVectorBlocked<...>
  } // namespace LAFEM
} // namespace FEAT