power_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/lafem/meta_element.hpp>
#include <kernel/lafem/container.hpp>
 
 
// includes, system
#include <iostream>
 
namespace FEAT
{
  namespace LAFEM
  {
    template<
      typename SubType_,
      int count_>
    class PowerVector
    {
    private:
      // Note: the case = 1 is specialized below
      static_assert(count_ > 1, "invalid block size");
 
      // declare this class template as a friend for recursive inheritance
      template<typename, int>
      friend class PowerVector;
 
      typedef PowerVector<SubType_, count_ - 1> RestClass;
 
      void _read_from_binary(std::istream& file)
      {
        //process first
        _first.read_from(FileMode::fm_binary, file);
 
        // append rest
        _rest._read_from_binary(file);
      }
 
      void _write_out_binary(std::ostream& file) const
      {
        //process first
        _first.write_out(FileMode::fm_binary, file);
 
        // append rest
        _rest._write_out_binary(file);
      }
 
    public:
      typedef SubType_ SubVectorType;
      typedef typename SubVectorType::DataType DataType;
      typedef typename SubVectorType::IndexType IndexType;
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = PowerVector<typename SubType_::template ContainerType<DT2_, IT2_>, count_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      static constexpr int num_blocks = count_;
 
    protected:
      SubVectorType _first;
      RestClass _rest;
 
    public:
      PowerVector()
      {
      }
 
      explicit PowerVector(Index sub_size) :
        _first(sub_size),
        _rest(sub_size)
      {
      }
 
      explicit PowerVector(Index sub_size, DataType value) :
        _first(sub_size, value),
        _rest(sub_size, value)
      {
      }
 
      explicit PowerVector(SubVectorType&& the_first, RestClass&& the_rest) :
        _first(std::move(the_first)),
        _rest(std::move(the_rest))
      {
      }
 
      PowerVector(PowerVector&& other) :
        _first(std::move(other._first)),
        _rest(std::move(other._rest))
      {
      }
 
      PowerVector& operator=(PowerVector&& other)
      {
        if(this != &other)
        {
          _first = std::move(other._first);
          _rest = std::move(other._rest);
        }
        return *this;
      }
 
      PowerVector(const PowerVector&) = delete;
      PowerVector& operator=(const PowerVector&) = delete;
 
      virtual ~PowerVector()
      {
      }
 
      PowerVector clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return PowerVector(_first.clone(mode), _rest.clone(mode));
      }
 
      void clone(const PowerVector& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
        _rest.clone(other._rest, clone_mode);
      }
 
      void clone(const PowerVector& other)
      {
        _first.clone(other._first);
        _rest.clone(other._rest);
      }
 
      template<int i_>
      SubVectorType& at()
      {
        static_assert((0 <= i_) && (i_ < count_), "invalid sub-vector index");
        return PowerElement<i_, SubVectorType>::get(*this);
      }
 
      template<int i_>
      const SubVectorType& at() const
      {
        static_assert((0 <= i_) && (i_ < count_), "invalid sub-vector index");
        return PowerElement<i_, SubVectorType>::get(*this);
      }
 
      SubVectorType& get(int i)
      {
        XASSERTM((0 <= i) && (i < count_), "invalid sub-vector index");
        return (i == 0) ? _first : _rest.get(i-1);
      }
 
      const SubVectorType& get(int i) const
      {
        XASSERTM((0 <= i) && (i < count_), "invalid sub-vector index");
        return (i == 0) ? _first : _rest.get(i-1);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return sizeof(std::uint64_t) + _first.get_checkpoint_size(config) + _rest.get_checkpoint_size(config); //sizeof(std::uint64_t) bits needed to store lenght of checkpointed _first
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        std::uint64_t isize = *(std::uint64_t*) data.data(); //get size of checkpointed _first
        std::vector<char>::iterator start = std::begin(data) + sizeof(std::uint64_t);
        std::vector<char>::iterator last_of_first = std::begin(data) + sizeof(std::uint64_t) + (int) isize;
        std::vector<char> buffer_first(start, last_of_first);
        _first.restore_from_checkpoint_data(buffer_first);
 
        data.erase(std::begin(data), last_of_first);
        _rest.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        std::size_t old_size = data.size();
        data.insert(std::end(data), sizeof(std::uint64_t), 0); //add placeholder
        std::uint64_t ireal_size = _first.set_checkpoint_data(data, config); //add data of _first to the overall checkpoint and save its size
        char* csize = reinterpret_cast<char*>(&ireal_size);
        for(std::size_t i(0) ; i < sizeof(std::uint64_t) ; ++i)  //overwrite the guessed datalength
        {
          data[old_size +i] = csize[i];
        }
 
        return sizeof(std::uint64_t) + ireal_size + _rest.set_checkpoint_data(data, config); //generate and add checkpoint data for the _rest
      }
 
      SubVectorType& first()
      {
        return _first;
      }
 
      const SubVectorType& first() const
      {
        return _first;
      }
 
      RestClass& rest()
      {
        return _rest;
      }
 
      const RestClass& rest() const
      {
        return _rest;
      }
 
      int blocks() const
      {
        return num_blocks;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index size() const
      {
        return first().template size<perspective_>() + rest().template size<perspective_>();
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
        rest().format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
        rest().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
        rest().clear();
      }
 
      static String name()
      {
        return String("PowerVector<") + SubVectorType::name() + "," + stringify(count_) + ">";
      }
 
      template<typename SubType2_>
      void copy(const PowerVector<SubType2_, count_>& x, bool full = false)
      {
        first().copy(x.first(), full);
        rest().copy(x.rest(), full);
      }
 
      void axpy(const PowerVector& x, DataType alpha = DataType(1))
      {
        first().axpy(x.first(), alpha);
        rest().axpy(x.rest(), alpha);
      }
 
      void component_product(const PowerVector & x, const PowerVector & y)
      {
        first().component_product(x.first(), y.first());
        rest().component_product(x.rest(), y.rest());
      }
 
      void component_invert(const PowerVector& x, DataType alpha = DataType(1))
      {
        first().component_invert(x.first(), alpha);
        rest().component_invert(x.rest(), alpha);
      }
 
      void scale(const PowerVector& x, DataType alpha)
      {
        first().scale(x.first(), alpha);
        rest().scale(x.rest(), alpha);
      }
 
      DataType dot(const PowerVector& x) const
      {
        return first().dot(x.first()) + rest().dot(x.rest());
      }
 
      DataType triple_dot(const PowerVector& x, const PowerVector& y) const
      {
        return first().triple_dot(x.first(), y.first())
          + rest().triple_dot(x.rest(), y.rest());
      }
 
      DataType norm2sqr() const
      {
        return first().norm2sqr() + rest().norm2sqr();
      }
 
      DataType norm2() const
      {
        return Math::sqrt(norm2sqr());
      }
 
      DataType max_abs_element() const
      {
        return Math::max(first().max_abs_element(), rest().max_abs_element());
      }
 
      DataType min_abs_element() const
      {
        return Math::min(first().min_abs_element(), rest().min_abs_element());
      }
 
      DataType max_element() const
      {
        return Math::max(first().max_element(), rest().max_element());
      }
 
      DataType min_element() const
      {
        return Math::min(first().min_element(), rest().min_element());
      }
 
      DataType max_rel_diff(const PowerVector & x) const
      {
        return Math::max(first().max_rel_diff(x.first()), rest().max_rel_diff(x.rest()));
      }
 
      bool same_layout(const PowerVector& x) const
      {
        return (first().same_layout(x.first())) && (rest().same_layout(x.rest()));
      }
 
      const DataType operator()(Index index) const
      {
        ASSERT(index < size());
 
        if (index < first().size())
        {
          return first()(index);
        }
        else
        {
          return rest()(index - first().size());
        }
      }
 
      void operator()(Index index, DataType value)
      {
        ASSERT(index < size());
 
        if (index < first().size())
        {
          first()(index, value);
        }
        else
        {
          rest()(index - first().size(), value);
        }
      }
 
      void set_vec(DataType * const pval_set) const
      {
        this->first().set_vec(pval_set);
        this->rest().set_vec(pval_set + this->first().template size<Perspective::pod>());
      }
 
      void set_vec_inv(const DataType * const pval_set)
      {
        this->first().set_vec_inv(pval_set);
        this->rest().set_vec_inv(pval_set + this->first().template size<Perspective::pod>());
      }
 
      template<typename SubType2_>
      void convert(const PowerVector<SubType2_, count_>& other)
      {
        this->first().convert(other.first());
        this->rest().convert(other.rest());
      }
 
      std::size_t bytes() const
      {
        return _first.bytes() + _rest.bytes();
      }
 
      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:
          {
            std::uint64_t magic; // magic_number
            file.read((char *)&magic, (long)(sizeof(std::uint64_t)));
            if (magic != 100)
              XABORTM("Given file or file component is no PowerVector!");
            std::uint64_t count; // subvector count
            file.read((char *)&count, (long)(sizeof(std::uint64_t)));
            if (count != count_)
              XABORTM("PowerVector file read in component count mismatch: class has " + stringify(count_) + "- " + stringify(count) + " read in!");
 
            _read_from_binary(file);
            break;
          }
          default:
            XABORTM("Filemode not supported!");
        }
      }
 
      void write_out(FileMode mode, const String& filename) const
      {
        std::ofstream file(filename.c_str(), std::ofstream::out | std::ofstream::binary);
        if (! file.is_open())
          XABORTM("Unable to open Vector file " + filename);
        write_out(mode, file);
        file.close();
      }
 
      void write_out(FileMode mode, std::ostream& file) const
      {
        switch(mode)
        {
          case FileMode::fm_binary:
          {
            size_t gsize(2 * sizeof(std::uint64_t)); // magic_number and subvector count
            std::vector<char> result(gsize);
            char * array(result.data());
            std::uint64_t * uiarray(reinterpret_cast<std::uint64_t *>(array));
            uiarray[0] = 100; 
            uiarray[1] = count_;
 
            file.write(result.data(), long(result.size()));
 
            _write_out_binary(file);
            break;
          }
          default:
            XABORTM("Filemode not supported!");
        }
      }
    }; // class PowerVector<...>
 
 
    template<typename SubType_>
    class PowerVector<SubType_, 1>
    {
    private:
 
      template<typename, int>
      friend class PowerVector;
 
      void _read_from_binary(std::istream& file)
      {
        //process first
        _first.read_from(FileMode::fm_binary, file);
      }
 
      void _write_out_binary(std::ostream& file) const
      {
        //process first
        _first.write_out(FileMode::fm_binary, file);
      }
 
    public:
      typedef SubType_ SubVectorType;
 
      typedef typename SubVectorType::DataType DataType;
      typedef typename SubVectorType::IndexType IndexType;
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = PowerVector<typename SubType_::template ContainerType<DT2_, IT2_>, Index(1)>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      static constexpr int num_blocks = 1;
 
    protected:
      SubVectorType _first;
 
    public:
      PowerVector()
      {
      }
 
      explicit PowerVector(Index sub_size)
        : _first(sub_size)
      {
      }
 
      explicit PowerVector(Index sub_size, DataType value)
        : _first(sub_size, value)
      {
      }
 
      explicit PowerVector(SubVectorType&& the_first) :
        _first(std::move(the_first))
      {
      }
 
      PowerVector(PowerVector&& other) :
        _first(std::move(other._first))
      {
      }
 
      PowerVector& operator=(PowerVector&& other)
      {
        if(this != &other)
        {
          _first = std::move(other._first);
        }
        return *this;
      }
 
      PowerVector(const PowerVector&) = delete;
      PowerVector& operator=(const PowerVector&) = delete;
 
      virtual ~PowerVector()
      {
      }
 
      PowerVector clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return PowerVector(_first.clone(mode));
      }
 
      void clone(const PowerVector& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
      }
 
      void clone(const PowerVector& other)
      {
        _first.clone(other._first);
      }
 
      template<int i>
      SubVectorType& at()
      {
        static_assert(i == 0, "invalid sub-vector index");
        return _first;
      }
 
      template<int i>
      const SubVectorType& at() const
      {
        static_assert(i == 0, "invalid sub-vector index");
        return _first;
      }
 
      SubVectorType& get(int i)
      {
        XASSERTM(i == 0, "invalid sub-vector index");
        return _first;
      }
 
      const SubVectorType& get(int i) const
      {
        XASSERTM(i == 0, "invalid sub-vector index");
        return _first;
      }
 
      SubVectorType& first()
      {
        return _first;
      }
 
      const SubVectorType& first() const
      {
        return _first;
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return _first.get_checkpoint_size(config);
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        _first.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        return _first.set_checkpoint_data(data, config);
      }
 
      int blocks() const
      {
        return 1;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index size() const
      {
        return _first.template size<perspective_>();
      }
 
      void format(DataType value = DataType(0))
      {
        _first.format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
      }
 
      static String name()
      {
        return String("PowerVector<") + SubVectorType::name() + ",1>";
      }
 
      template<typename SubType2_>
      void copy(const PowerVector<SubType2_, 1>& x, bool full = false)
      {
        first().copy(x.first(), full);
      }
 
      void axpy(const PowerVector& x, DataType alpha = DataType(1))
      {
        first().axpy(x.first(), alpha);
      }
 
      void component_product(const PowerVector & x, const PowerVector & y)
      {
        first().component_product(x.first(), y.first());
      }
 
      void component_invert(const PowerVector& x, DataType alpha = DataType(1))
      {
        first().component_invert(x.first(), alpha);
      }
 
      void scale(const PowerVector& x, DataType alpha)
      {
        first().scale(x.first(), alpha);
      }
 
      DataType dot(const PowerVector& x) const
      {
        return first().dot(x.first());
      }
 
      DataType triple_dot(const PowerVector& x, const PowerVector& y) const
      {
        return first().triple_dot(x.first(), y.first());
      }
 
      DataType triple_dot_i(const PowerVector& x, const PowerVector& y) const
      {
        return first().triple_dot_i(x.first(), y.first());
      }
 
      DataType norm2sqr() const
      {
        return first().norm2sqr();
      }
 
      DataType norm2() const
      {
        return Math::sqrt(norm2sqr());
      }
 
      DataType max_abs_element() const
      {
        return first().max_abs_element();
      }
 
      DataType min_abs_element() const
      {
        return first().min_abs_element();
      }
 
      DataType max_element() const
      {
        return first().max_element();
      }
 
      DataType min_element() const
      {
        return first().min_element();
      }
 
      const DataType operator()(Index index) const
      {
        ASSERT(index < size());
 
        return first()(index);
      }
 
      void operator()(Index index, DataType value)
      {
        ASSERT(index < size());
 
        first()(index, value);
      }
 
      void set_vec(DataType * const pval_set) const
      {
        this->first().set_vec(pval_set);
      }
 
      void set_vec_inv(const DataType * const pval_set)
      {
        this->first().set_vec_inv(pval_set);
      }
 
      template<typename SubType2_>
      void convert(const PowerVector<SubType2_, 1>& other)
      {
        this->first().convert(other.first());
      }
 
      std::size_t bytes() const
      {
        return _first.bytes();
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        switch(mode)
        {
        case FileMode::fm_binary:
          read_from_binary(filename);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      void read_from(FileMode mode, std::istream& file)
      {
        switch(mode)
        {
        case FileMode::fm_binary:
          read_from_binary(file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      void read_from_binary(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_binary(file);
        file.close();
      }
 
      void read_from_binary(std::istream& file)
      {
        std::uint64_t magic; // magic_number
        file.read((char *)&magic, (long)(sizeof(std::uint64_t)));
        if (magic != 100)
          XABORTM("Given file or file component is no PowerVector!");
        std::uint64_t count; // subvector count
        file.read((char *)&count, (long)(sizeof(std::uint64_t)));
        if (count != 1)
          XABORTM("PowerVector file read in component count mismatch: class has 1 - " + stringify(count) + " read in!");
 
        _read_from_binary(file);
      }
 
      void write_out(FileMode mode, const String& filename) const
      {
        switch(mode)
        {
        case FileMode::fm_binary:
          write_out_binary(filename);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      void write_out(FileMode mode, std::ostream& file) const
      {
        switch(mode)
        {
        case FileMode::fm_binary:
          write_out_binary(file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      void write_out_binary(const String& filename) const
      {
        std::ofstream file(filename.c_str(), std::ofstream::out | std::ofstream::binary);
        if (! file.is_open())
          XABORTM("Unable to open Vector file " + filename);
        write_out_binary(file);
        file.close();
      }
 
      void write_out_binary(std::ostream& file) const
      {
        size_t gsize(2 * sizeof(std::uint64_t)); // magic_number and subvector count
        std::vector<char> result(gsize);
        char * array(result.data());
        std::uint64_t * uiarray(reinterpret_cast<std::uint64_t *>(array));
        uiarray[0] = 100; 
        uiarray[1] = 1; //fixed count_ for power vector specialization
 
        file.write(result.data(), long(result.size()));
 
        _write_out_binary(file);
      }
    }; // class PowerVector<...,1>
 
    template <typename SubType_, Index count_>
    inline void dump_power_vector(std::ostream & os, const PowerVector<SubType_, count_>& x)
    {
      os << "," << x.first();
      dump_power_vector(os, x.rest());
    }
 
    template <typename SubType_>
    inline void dump_power_vector(std::ostream & os, const PowerVector<SubType_, 1>& x)
    {
      os << x.first();
    }
 
    template <typename SubType_, Index count_>
    inline std::ostream& operator<< (std::ostream & os, const PowerVector<SubType_, count_>& x)
    {
      os << "[";
      dump_power_vector(os, x);
      os << "]";
      return os;
    }
  } // namespace LAFEM
} // namespace FEAT