dense_vector.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/lafem/forward.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/type_traits.hpp>
#include <kernel/util/math.hpp>
#include <kernel/util/random.hpp>
#include <kernel/lafem/container.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/lafem/arch/component_invert.hpp>
#include <kernel/lafem/arch/dot_product.hpp>
#include <kernel/lafem/arch/norm.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/lafem/arch/scale.hpp>
#include <kernel/lafem/arch/axpy.hpp>
#include <kernel/lafem/arch/component_product.hpp>
#include <kernel/adjacency/permutation.hpp>
#include <kernel/util/statistics.hpp>
#include <kernel/util/time_stamp.hpp>
 
#include <iostream>
#include <fstream>
#include <string>
#include <stdint.h>
 
namespace FEAT
{
  namespace LAFEM
  {
    template <typename DT_, typename IT_ = Index>
    class DenseVector : public Container<DT_, IT_>
    {
    public:
      class ScatterAxpy
      {
      public:
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
        typedef DT_ DataType;
        typedef IT_ IndexType;
 
      private:
        Index _num_entries;
        DT_ * _data;
 
      public:
        explicit ScatterAxpy(VectorType & vector) :
          _num_entries(vector.size()),
          _data(vector.elements())
        {
        }
 
        template <typename LocalVector_, typename Mapping_>
        void operator()(const LocalVector_ & loc_vec, const Mapping_ & mapping, DT_ alpha = DT_(1))
        {
          // loop over all local entries
          for (int i(0); i < mapping.get_num_local_dofs(); ++i)
          {
            // get dof index
            Index dof_idx = mapping.get_index(i);
            ASSERT(dof_idx < _num_entries);
 
            // update vector entry
            _data[dof_idx] += alpha * loc_vec[i];
          }
        }
      }; // class ScatterAxpy
 
      class GatherAxpy
      {
      public:
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
        typedef DT_ DataType;
        typedef IT_ IndexType;
 
      private:
        Index _num_entries;
        const DT_ * _data;
 
      public:
        explicit GatherAxpy(const VectorType & vector) :
          _num_entries(vector.size()),
          _data(vector.elements())
        {
        }
 
        template <typename LocalVector_, typename Mapping_>
        void operator()(LocalVector_ & loc_vec, const Mapping_ & mapping, DT_ alpha = DT_(1))
        {
          // loop over all local entries
          for (int i(0); i < mapping.get_num_local_dofs(); ++i)
          {
            // get dof index
            Index dof_idx = mapping.get_index(i);
            ASSERT(dof_idx < _num_entries);
 
            // update local vector data
            loc_vec[i] += alpha * _data[dof_idx];
          }
        }
      }; // class GatherAxpy
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      typedef DT_ ValueType;
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using ContainerType = DenseVector<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      explicit DenseVector() :
        Container<DT_, IT_> (0)
      {
      }
 
      explicit DenseVector(Index size_in) :
        Container<DT_, IT_>(size_in)
      {
        if (size_in == Index(0))
          return;
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(size_in));
 
        this->_elements_size.push_back(size_in);
      }
 
      explicit DenseVector(Index size_in, DT_ value) :
        Container<DT_, IT_>(size_in)
      {
        if (size_in == Index(0))
          return;
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(size_in));
        this->_elements_size.push_back(size_in);
 
        MemoryPool::set_memory(this->_elements.at(0), value, size_in);
      }
 
      explicit DenseVector(Index size_in, DT_ * data) :
        Container<DT_, IT_>(size_in)
      {
        if (size_in == Index(0))
          return;
 
        this->_elements.push_back(data);
        this->_elements_size.push_back(size_in);
 
        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));
      }
 
      explicit DenseVector(const DenseVector & dv_in, Index size_in, Index offset_in) :
        Container<DT_, IT_>(size_in)
      {
        XASSERT(size_in > Index(0));
        XASSERTM(size_in + offset_in <= dv_in.size(), "Ranged vector part exceeds original vector size!");
 
        this->_foreign_memory = true;
 
        DT_ * te(const_cast<DT_ *>(dv_in.elements()));
        this->_elements.push_back(te + offset_in);
        this->_elements_size.push_back(size_in);
      }
 
      template <int BS_>
      explicit DenseVector(const DenseVectorBlocked<DT_, IT_, BS_> & other) :
        Container<DT_, IT_>(other.template size<Perspective::pod>())
      {
        convert(other);
      }
 
      explicit DenseVector(FileMode mode, String filename) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, filename);
      }
 
      explicit DenseVector(FileMode mode, std::istream& file) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, file);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      explicit DenseVector(std::vector<char> input) :
        Container<DT_, IT_>(0)
      {
        deserialize<DT2_, IT2_>(input);
      }
 
      explicit DenseVector(Random & rng, Index size_in, DataType min, DataType max) :
        Container<DT_, IT_>(size_in)
      {
        if (size_in == Index(0))
          return;
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(size_in));
        this->_elements_size.push_back(size_in);
 
        this->format(rng, min, max);
      }
 
      DenseVector(DenseVector && other) :
        Container<DT_, IT_>(std::forward<DenseVector>(other))
      {
      }
 
      virtual ~DenseVector()
      {
      }
 
      DenseVector & operator=(DenseVector && other)
      {
        this->move(std::forward<DenseVector>(other));
 
        return *this;
      }
 
      DenseVector clone(CloneMode clone_mode = CloneMode::Deep) const
      {
        DenseVector t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template <typename DT2_, typename IT2_>
      void clone(const DenseVector<DT2_, IT2_> & other, CloneMode clone_mode = CloneMode::Deep)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const DenseVector<DT2_, IT2_> & other)
      {
        this->assign(other);
      }
 
      template <typename DT2_, typename IT2_, int BS2_>
      void convert(const DenseVectorBlocked< DT2_, IT2_, BS2_> & other)
      {
        this->clear();
 
        this->_scalar_index.push_back(other.template size<Perspective::pod>());
        this->_elements.push_back(other.get_elements().at(0));
        this->_elements_size.push_back(this->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));
      }
 
      template <typename VT_>
      void convert(const VT_ & a)
      {
        DenseVector vec(a.template size<Perspective::pod>());
        a.set_vec(vec.elements());
        this->move(std::move(vec));
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      void deserialize(std::vector<char> input)
      {
        this->template _deserialize<DT2_, IT2_>(FileMode::fm_dv, input);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      std::vector<char> serialize(const LAFEM::SerialConfig& config = LAFEM::SerialConfig()) const
      {
        return this->template _serialize<DT2_, IT2_>(FileMode::fm_dv, config);
      }
 
      template <int BlockSize_>
      DenseVectorBlocked<DT_, IT_, BlockSize_> inflate_to_blocks()
      {
        DenseVectorBlocked<DT_, IT_, BlockSize_> result(this->size());
 
        const auto* p = elements();
        auto* rp = result.template elements<Perspective::native>();
        for (Index i(0) ; i < this->size() ; ++i)
        {
          rp[i] = typename decltype(result)::ValueType(p[i]);
        }
        return result;
      }
 
      template <typename VT_>
      void copy(const VT_ & a)
      {
        XASSERTM(this->template size<Perspective::pod>() == a.template size<Perspective::pod>(), "Vectors have not the same size!");
 
        a.set_vec(this->elements());
      }
 
      template <typename VT_>
      void copy_inv(VT_ & a) const
      {
        XASSERTM(this->template size<Perspective::pod>() == a.template size<Perspective::pod>(), "Vectors have not the same size!");
 
        a.set_vec_inv(this->elements());
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        std::ios_base::openmode bin = std::ifstream::in | std::ifstream::binary;
        if (mode == FileMode::fm_mtx)
          bin = std::ifstream::in;
        std::ifstream file(filename.c_str(), bin);
        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_mtx:
        {
          this->clear();
          this->_scalar_index.push_back(0);
 
          Index rows;
          String line;
          std::getline(file, line);
          if (line.find("%%MatrixMarket matrix array real general") == String::npos)
          {
            XABORTM("Input-file is not a compatible mtx-vector-file");
          }
          while (! file.eof())
          {
            std::getline(file, line);
            if (file.eof())
              XABORTM("Input-file is empty");
 
            String::size_type begin(line.find_first_not_of(" "));
            if (line.at(begin) != '%')
              break;
          }
          {
            String::size_type begin(line.find_first_not_of(" "));
            line.erase(0, begin);
            String::size_type end(line.find_first_of(" "));
            String srows(line, 0, end);
            rows = (Index)atol(srows.c_str());
            line.erase(0, end);
 
            begin = line.find_first_not_of(" ");
            line.erase(0, begin);
            end = line.find_first_of(" ");
            String scols(line, 0, end);
            Index cols((Index)atol(scols.c_str()));
            line.erase(0, end);
            XASSERTM(cols == 1, "Input-file is no dense-vector-file");
          }
 
          DenseVector<DT_, IT_> tmp(rows);
          DT_ * pval(tmp.elements());
 
          while (! file.eof())
          {
            std::getline(file, line);
            if (file.eof())
              break;
 
            String::size_type begin(line.find_first_not_of(" "));
            line.erase(0, begin);
            String::size_type end(line.find_first_of(" "));
            String sval(line, 0, end);
            DT_ tval((DT_)atof(sval.c_str()));
 
            *pval = tval;
            ++pval;
          }
          this->move(std::move(tmp));
          break;
        }
        case FileMode::fm_exp:
        {
          this->clear();
          this->_scalar_index.push_back(0);
 
          std::vector<DT_> data;
 
          while (! file.eof())
          {
            std::string line;
            std::getline(file, line);
            if (line.find("#", 0) < line.npos)
              continue;
            if (file.eof())
              break;
 
            std::string n_z_s;
 
            std::string::size_type first_digit(line.find_first_not_of(" "));
            line.erase(0, first_digit);
            std::string::size_type eol(line.length());
            for (std::string::size_type i(0); i < eol; ++i)
            {
              n_z_s.append(1, line[i]);
            }
 
            DT_ n_z((DT_)atof(n_z_s.c_str()));
 
            data.push_back(n_z);
          }
 
          this->_scalar_index.at(0) = Index(data.size());
          this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(Index(data.size())));
          this->_elements_size.push_back(Index(data.size()));
          MemoryPool::copy(this->_elements.at(0), &data[0], Index(data.size()));
          break;
        }
        case FileMode::fm_dv:
        case FileMode::fm_binary:
          this->template _deserialize<double, std::uint64_t>(FileMode::fm_dv, 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 || mode == FileMode::fm_exp)
          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:
        {
          DenseVector<DT_, IT_> temp;
          temp.convert(*this);
 
          const Index tsize(temp.size());
          file << "%%MatrixMarket matrix array real general\n";
          file << tsize << " " << 1 << "\n";
 
          const DT_ * pval(temp.elements());
          for (Index i(0); i < tsize; ++i, ++pval)
          {
            file << stringify_fp_sci(*pval) << "\n";
          }
          break;
        }
        case FileMode::fm_exp:
        {
          DT_ * temp = MemoryPool::template allocate_memory<DT_>((this->size()));
          MemoryPool::copy(temp, this->_elements.at(0), this->size());
 
          for (Index i(0); i < this->size(); ++i)
          {
            file << stringify_fp_sci(temp[i]) << "\n";
          }
          MemoryPool::release_memory(temp);
          break;
        }
        case FileMode::fm_dv:
        case FileMode::fm_binary:
          this->template _serialize<double, std::uint64_t>(FileMode::fm_dv, file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      template <Perspective = Perspective::native>
      DT_ * elements()
      {
        if (this->_elements.size() == 0)
          return nullptr;
 
        return this->_elements.at(0);
      }
 
      template <Perspective = Perspective::native>
      DT_ const * elements() const
      {
        if (this->_elements.size() == 0)
          return nullptr;
 
        return this->_elements.at(0);
      }
 
      const DT_ operator()(Index index) const
      {
        ASSERT(index < this->size());
 
        MemoryPool::synchronize();
 
        return elements()[index];
      }
 
      void operator()(Index index, DT_ value)
      {
        ASSERT(index < this->size());
        elements()[index] = value;
        MemoryPool::synchronize();
      }
 
      static String name()
      {
        return "DenseVector";
      }
 
      void copy(const DenseVector & x, bool full = false)
      {
        this->_copy_content(x, full);
      }
 
 
      void axpy(
        const DenseVector & x,
        const DT_ alpha = DT_(1))
      {
        XASSERTM(x.size() == this->size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size() * 2);
        Arch::Axpy::value(this->elements(), alpha, x.elements(), this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void component_product(const DenseVector & x, const DenseVector & y)
      {
        XASSERTM(this->size() == x.size(), "Vector size does not match!");
        XASSERTM(this->size() == y.size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size());
        Arch::ComponentProduct::value(this->elements(), x.elements(), y.elements(), this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void component_invert(const DenseVector & x, const DT_ alpha = DT_(1))
      {
        XASSERTM(this->size() == x.size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size());
        Arch::ComponentInvert::value(this->elements(), x.elements(), alpha, this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void scale(const DenseVector & x, const DT_ alpha)
      {
        XASSERTM(x.size() == this->size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size());
        Arch::Scale::value(this->elements(), x.elements(), alpha, this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      DataType triple_dot(const DenseVector & x, const DenseVector & y) const
      {
        XASSERTM(x.size() == this->size(), "Vector size does not match!");
        XASSERTM(y.size() == this->size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size() * 3);
        DataType result = Arch::TripleDotProduct::value(this->elements(), x.elements(), y.elements(), this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DataType dot(const DenseVector & x) const
      {
        XASSERTM(x.size() == this->size(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->size() * 2);
        DataType result = Arch::DotProduct::value(this->elements(), x.elements(), this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ norm2() const
      {
        TimeStamp ts_start;
        Statistics::add_flops(this->size() * 2);
 
        DT_ result = Arch::Norm2::value(this->elements(), this->size());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ norm2sqr() const
      {
        // fallback
        return Math::sqr(this->norm2());
      }
 
      DT_ max_abs_element() const
      {
        TimeStamp ts_start;
 
        Index max_abs_index = Arch::MaxAbsIndex::value(this->template elements<Perspective::pod>(), this->template size<Perspective::pod>());
        ASSERT(max_abs_index < this->template size<Perspective::pod>());
        DT_ result;
        MemoryPool::copy(&result, this->template elements<Perspective::pod>() + max_abs_index, 1);
        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 size<Perspective::pod>());
        ASSERT(min_abs_index < this->template size<Perspective::pod>());
        DT_ result;
        MemoryPool::copy(&result, this->template elements<Perspective::pod>() + min_abs_index, 1);
        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 size<Perspective::pod>());
        ASSERT(max_index < this->template size<Perspective::pod>());
        DT_ result;
        MemoryPool::copy(&result, this->template elements<Perspective::pod>() + max_index, 1);
 
        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 size<Perspective::pod>());
        ASSERT(min_index < this->template size<Perspective::pod>());
        DT_ result;
        MemoryPool::copy(&result, this->template elements<Perspective::pod>() + min_index, 1);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_rel_diff(const DenseVector & 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 size<Perspective::pod>());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return max_rel_diff;
      }
 
      bool same_layout(const DenseVector& x) const
      {
        if (this->size() == 0 && x.size() == 0 && this->get_elements().size() == 0 && x.get_elements().size() == 0)
          return true;
        if (this->size() != x.size())
          return false;
        if (this->get_elements().size() != x.get_elements().size())
          return false;
        if (this->get_indices().size() != x.get_indices().size())
          return false;
 
        return true;
      }
 
 
      void permute(Adjacency::Permutation & perm)
      {
        if (perm.size() == 0)
          return;
 
        XASSERTM(perm.size() == this->size(), "Container size does not match permutation size");
 
        perm.apply(elements());
      }
 
      void set_vec(DT_ * const pval_set) const
      {
        MemoryPool::copy(pval_set, this->elements(), this->size());
      }
 
      void set_vec_inv(const DT_ * const pval_set)
      {
        MemoryPool::copy(this->elements(), pval_set, this->size());
      }
 
 
      friend std::ostream & operator<<(std::ostream & lhs, const DenseVector & b)
      {
        lhs << "[";
        for (Index i(0); i < b.size(); ++i)
        {
          lhs << "  " << stringify(b(i));
        }
        lhs << "]";
 
        return lhs;
      }
    }; // class DenseVector<...>
 
  #ifdef FEAT_EICKT
    extern template class DenseVector<float, std::uint32_t>;
    extern template class DenseVector<double, std::uint32_t>;
    extern template class DenseVector<float, std::uint64_t>;
    extern template class DenseVector<double, std::uint64_t>;
  #endif
 
  } // namespace LAFEM
} // namespace FEAT