SF_abstract_vector.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#ifndef _SF_ABSTRACT_VECTOR_H
#define _SF_ABSTRACT_VECTOR_H
 
#include <mpi.h>
 
#include <petscvec.h> // TODO: For scattering, VecScatter and friends
#include <petscis.h>
#include <petscsys.h> // TODO: For PETSC_COMM_WORLD
 
#include "SF_container.h"
#include "SF_globals.h"
#include "SF_parallel_layout.h"
#include "SF_parallel_utils.h"
#include "mpi_utils.h"
#include "petsc_utils.h"
 
#include "hashmap.hpp"
 
namespace SF {
 
// Forward declare the scattering class
class scattering;
 
template<class T, class S>
class abstract_vector {
 
  public:
 
  enum ltype {algebraic, nodal, elemwise, unset};
 
  const meshdata<mesh_int_t,mesh_real_t> * mesh   = NULL;    
  int                   dpn    = 0;       
  ltype                 layout = unset;   
 
  virtual ~abstract_vector() = default;
 
  virtual void init(const meshdata<mesh_int_t, mesh_real_t>& imesh,
                    int                                      idpn,
                    ltype                                    inp_layout) = 0;
 
  virtual void init(const abstract_vector<T, S>& vec) = 0;
 
  virtual void init(T igsize, T ilsize, int idpn = 1, ltype ilayout = unset) = 0;
 
  inline std::tuple<T, T> init_common(const meshdata<mesh_int_t, mesh_real_t>& imesh,
                                      int                                      idpn,
                                      ltype                                    inp_layout)
  {
    mesh = &imesh;
    dpn = idpn;
    T N = 0, n = 0;
 
    switch(inp_layout) {
      case algebraic:
        N = static_cast<T>(mesh->pl.num_global_idx()) * static_cast<T>(dpn);
        n = static_cast<T>(mesh->pl.num_algebraic_idx()) * static_cast<T>(dpn);
        layout = algebraic;
        break;
 
      case nodal:
        N = static_cast<T>(mesh->l_numpts) * static_cast<T>(dpn);
        MPI_Allreduce(MPI_IN_PLACE, &N, 1, opencarp::mpi_datatype<T>(), MPI_SUM, mesh->comm);
        n = static_cast<T>(mesh->l_numpts) * static_cast<T>(dpn);
        layout = nodal;
        break;
 
      case elemwise:
        N = static_cast<T>(mesh->g_numelem*dpn);
        n = static_cast<T>(mesh->l_numelem*dpn);
        layout = elemwise;
        break;
 
      default: break;
    }
 
    return std::tuple<T, T>(N, n);
  }
 
  virtual void set(const vector<T>& idx, const vector<S>& vals, const bool additive = false, const bool local = false) = 0;
 
  virtual void set(const vector<T>& idx, const S val, const bool additive = false, const bool local = false) = 0;
 
  virtual void set(const S val) = 0;
 
  virtual void set(const T idx, const S val) = 0;
 
  virtual void get(const vector<T> & idx, S *out) = 0;
 
  virtual S get(const T idx) = 0;
 
  virtual void operator*=(const S sca) = 0;
 
  virtual void operator /= (const S sca) = 0;
 
  virtual void operator*=(const abstract_vector<T, S>& vec) = 0;
 
  virtual void add_scaled(const abstract_vector<T, S>& vec, S k) = 0;
 
  virtual void operator += (const abstract_vector<T,S> & vec) = 0;
 
  virtual void operator-=(const abstract_vector<T, S>& vec) = 0;
 
  virtual void operator+=(S c) = 0;
 
  virtual void operator=(const vector<S>& rhs) = 0;
 
  virtual void operator=(const abstract_vector<T, S>& rhs) = 0;
 
  virtual void shallow_copy(const abstract_vector<T, S>& v) = 0;
 
  virtual void deep_copy(const abstract_vector<T, S>& v) = 0;
 
  virtual void overshadow(const abstract_vector<T, S>& sub, bool member, int offset, int sz, bool share) = 0;
 
  virtual T lsize() const = 0;
 
  virtual T gsize() const = 0;
 
  virtual void get_ownership_range(T& start, T& stop) const = 0;
 
  virtual S* ptr() = 0;
 
  virtual const S* const_ptr() const = 0;
 
  virtual S* device_ptr() = 0;
 
  virtual void release_ptr(S*& p) = 0;
 
  virtual void const_release_ptr(const S*& p) const = 0;
 
  virtual void release_device_ptr(S*& p) = 0;
 
  virtual S mag() const = 0;
 
  virtual S sum() const = 0;
 
  virtual S min() const = 0;
 
  virtual S dot(const abstract_vector<T, S>& v) const = 0;
 
  virtual bool is_init() const = 0;
 
  virtual std::string to_string() const = 0;
 
  virtual bool equals(const abstract_vector<T,S> & rhs) const = 0;
 
  virtual void finish_assembly() = 0;
 
  virtual void forward(abstract_vector<T, S>& out, scattering &sc, bool add = false) = 0;
 
  virtual void backward(abstract_vector<T,S> & out, scattering &sc, bool add = false) = 0;
 
  virtual void apply_scattering(scattering& sc, bool fwd) = 0;
 
  inline size_t write_ascii(const char* file, bool write_header)
  {
    int size, rank;
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
    MPI_Comm_rank(comm, &rank);
    MPI_Comm_size(comm, &size);
    const MPI_Datatype mpi_t = opencarp::mpi_datatype<T>();
 
    T glb_size = this->gsize();
    T loc_size = this->lsize();
 
    int err = 0;
    FILE* fd = NULL;
    long int nwr = 0;
 
    if(rank == 0) {
      fd = fopen(file, "w");
      if(!fd) err = 1;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_MAX, comm);
    if(err) {
      treat_file_open_error(file, __func__, errno, false, rank);
      return nwr;
    }
 
    S* p = this->ptr();
 
    if(rank == 0) {
      if(write_header)
        fprintf(fd, "%jd\n", opencarp::printable_int(glb_size));
 
      for(T i=0; i<loc_size/this->dpn; i++) {
        for(int j=0; j<this->dpn; j++)
          nwr += fprintf(fd, "%f ", p[i*this->dpn+j]);
        nwr += fprintf(fd, "\n");
      }
 
      vector<S> wbuff;
      for(int pid=1; pid < size; pid++)
      {
        T rsize;
        MPI_Status stat;
 
        MPI_Recv(&rsize, 1, mpi_t, pid, SF_MPITAG, comm, &stat);
        wbuff.resize(rsize);
        MPI_Recv(wbuff.data(), rsize*sizeof(S), MPI_BYTE, pid, SF_MPITAG, comm, &stat);
 
        for(T i=0; i<rsize/this->dpn; i++) {
          for(int j=0; j<this->dpn; j++)
            nwr += fprintf(fd, "%f ", wbuff[i*this->dpn+j]);
          nwr += fprintf(fd, "\n");
        }
      }
      fclose(fd);
    }
    else {
      MPI_Send(&loc_size, 1, mpi_t, 0, SF_MPITAG, comm);
      MPI_Send(p, loc_size*sizeof(S), MPI_BYTE, 0, SF_MPITAG, comm);
    }
 
    this->release_ptr(p);
    MPI_Bcast(&nwr, 1, MPI_LONG, 0, comm);
 
    return nwr;
  }
 
  template<typename V>
  inline size_t write_binary(FILE* fd)
  {
    int size, rank;
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
    MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &size);
 
    long int loc_size = this->lsize();
    S* p = this->ptr();
 
    vector<V> buff(loc_size);
    for(long int i=0; i<loc_size; i++) buff[i] = p[i];
 
    long int nwr = root_write(fd, buff, comm);
 
    this->release_ptr(p);
    return nwr;
  }
 
  template<typename V>
  inline size_t write_binary(std::string file)
  {
    size_t nwr = 0;
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
    int rank; MPI_Comm_rank(comm, &rank);
 
    FILE* fd = NULL;
    int error = 0;
 
    if(rank == 0) {
      fd = fopen(file.c_str(), "w");
      if(fd == NULL)
        error++;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &error, 1, MPI_INT, MPI_SUM, comm);
 
    if(error == 0) {
      nwr = this->write_binary<V>(fd);
      fclose(fd);
    }
    else {
      treat_file_open_error(file.c_str(), __func__, errno, false, rank);
    }
 
    return nwr;
  }
 
  template<typename V>
  inline size_t read_binary(FILE* fd)
  {
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
 
    size_t loc_size = this->lsize();
    S* p = this->ptr();
    vector<V> buff(loc_size);
 
    size_t nrd = root_read(fd, buff, comm);
 
    for(size_t i=0; i<loc_size; i++)
      p[i] = buff[i];
 
    this->release_ptr(p);
 
    return nrd;
  }
 
  template<typename V>
  inline size_t read_binary(std::string file)
  {
    MPI_Comm comm = mesh != NULL ? mesh->comm : PETSC_COMM_WORLD;
 
    size_t nrd = 0;
    int rank; MPI_Comm_rank(comm, &rank);
 
    FILE* fd = NULL;
    int error = 0;
 
    if(rank == 0) {
      fd = fopen(file.c_str(), "r");
      if(fd == NULL)
        error++;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &error, 1, MPI_INT, MPI_SUM, comm);
 
    if(error == 0) {
      nrd = read_binary<V>(fd);
      fclose(fd);
    }
    else {
      treat_file_open_error(file.c_str(), __func__, errno, false, rank);
    }
 
    return nrd;
  }
 
  inline size_t read_ascii(FILE* fd)
  {
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
 
    size_t loc_size = this->lsize();
    S* p = this->ptr();
 
    size_t nrd = root_read_ascii(fd, p, loc_size, comm, false);
    this->release_ptr(p);
    return nrd;
  }
 
  inline size_t read_ascii(std::string file)
  {
    MPI_Comm comm = this->mesh != NULL ? this->mesh->comm : PETSC_COMM_WORLD;
    int rank; MPI_Comm_rank(comm, &rank);
 
    size_t nrd = 0;
 
    FILE* fd = NULL;
    int error = 0;
 
    if(rank == 0) {
      fd = fopen(file.c_str(), "r");
      if(fd == NULL)
        error++;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &error, 1, MPI_INT, MPI_SUM, comm);
 
    if(error == 0) {
      nrd = read_ascii(fd);
      if(fd) fclose(fd);
    }
    else {
      treat_file_open_error(file.c_str(), __func__, errno, false, rank);
    }
 
    return nrd;
  }
};
 
 
template<class T, class S>
inline bool is_init(const abstract_vector<T,S>* v)
{
  return (v && v->is_init());
}
 
 
// TODO: This is a type of Vector and can most likely be integrated into the
// `abstract_vector` type in some way. Everything here and also in
// `SF_parallel_layout.h` probably need to be abstracted away, or made less
// PETSc specific?
class scattering
{
  public:
  Vec        b_buff;  
  VecScatter vec_sc;  
 
  vector<SF_int> idx_a, idx_b;
 
  scattering() : b_buff(NULL), vec_sc(NULL)
  {}
 
  ~scattering()
  {
    if(b_buff) PETSC_CHECK(VecDestroy(&b_buff));
    if(vec_sc) PETSC_CHECK(VecScatterDestroy(&vec_sc));
  }
  template<class T, class S>
  inline void forward(abstract_vector<T,S> & in, abstract_vector<T,S> & out, bool add = false)
  {
    in.forward(out, *this, add);
  }
 
  template<class T, class S>
  inline void backward(abstract_vector<T,S> & in, abstract_vector<T,S> & out, bool add = false)
  {
      in.backward(out, *this, add);
  }
 
  template<class T, class S>
  inline void operator()(abstract_vector<T,S> & v, bool fwd)
  {
    v.apply_scattering(*this, fwd);
  }
};
 
class scatter_registry
{
  private:
  hashmap::unordered_map<quadruple<int>, scattering*> _registry;
 
  public:
  template<class T> inline scattering*
  register_permutation(const quadruple<int> spec,
                       const vector<T> & nbr_a,
                       const vector<T> & nbr_b,
                       const size_t gsize_a,
                       const size_t gsize_b,
                       const short dpn)
  {
    assert(_registry.count(spec) == 0);
 
    scattering* sc = new scattering();
    _registry[spec] = sc;
 
    IS is_a, is_b;
    vector<SF_int> idx_a(nbr_a.size() * dpn);
    vector<SF_int> idx_b(nbr_b.size() * dpn);
 
    // we copy indexing into new containers because of dpn and also
    // because of the implicit typecast between T and SF_int
    for(size_t i=0; i<nbr_a.size(); i++)
      for(short j=0; j<dpn; j++) idx_a[i*dpn+j] = nbr_a[i]*dpn+j;
 
    for(size_t i=0; i<nbr_b.size(); i++)
      for(short j=0; j<dpn; j++) idx_b[i*dpn+j] = nbr_b[i]*dpn+j;
 
    PETSC_CHECK(ISCreateGeneral(PETSC_COMM_WORLD, idx_a.size(), idx_a.data(), PETSC_COPY_VALUES, &is_a));
    PETSC_CHECK(ISCreateGeneral(PETSC_COMM_WORLD, idx_b.size(), idx_b.data(), PETSC_COPY_VALUES, &is_b));
 
    PETSC_CHECK(ISSetPermutation(is_b));
 
    Vec a;
    PETSC_CHECK(VecCreateMPI(PETSC_COMM_WORLD, idx_a.size(), gsize_a*dpn, &a));
    PETSC_CHECK(VecCreateMPI(PETSC_COMM_WORLD, idx_b.size(), gsize_b*dpn, &sc->b_buff));
    PETSC_CHECK(VecSetFromOptions(a));
    PETSC_CHECK(VecSetFromOptions(sc->b_buff));
 
    PETSC_CHECK(VecScatterCreate(a, is_a, sc->b_buff, is_b, &sc->vec_sc));
 
    PETSC_CHECK(VecDestroy(&a));
    PETSC_CHECK(ISDestroy (&is_a));
    PETSC_CHECK(ISDestroy (&is_b));
 
    sc->idx_a.assign(idx_a.begin(), idx_a.end());
    sc->idx_b.assign(idx_b.begin(), idx_b.end());
 
    return sc;
  }
 
 
  template<class T> inline scattering*
  register_scattering(const quadruple<int> spec,
                      const vector<T> & layout_a,
                      const vector<T> & layout_b,
                      const vector<T> & idx_a,
                      const vector<T> & idx_b,
                      const int rank,
                      const int dpn)
  {
    // scattering spec must not exist yet
    assert(_registry.count(spec) == 0);
 
    // the local index sets have to be of equal size.
    assert(idx_a.size() == idx_b.size());
 
    scattering* sc = new scattering();
    _registry[spec] = sc;
 
    IS is_a, is_b;
    Vec a;
    vector<SF_int> sidx_a(idx_a.size() * dpn);
    vector<SF_int> sidx_b(idx_b.size() * dpn);
 
    T lsize_a = layout_a[rank+1] - layout_a[rank];
    T lsize_b = layout_b[rank+1] - layout_b[rank];
    T gsize_a = layout_a[layout_a.size()-1];
    T gsize_b = layout_b[layout_a.size()-1];
 
    // we copy indexing into new containers because of dpn and also
    // because of the implicit typecast between T and SF_int
    for(size_t i=0; i<idx_a.size(); i++)
      for(short j=0; j<dpn; j++) sidx_a[i*dpn+j] = idx_a[i]*dpn+j;
 
    for(size_t i=0; i<idx_b.size(); i++)
      for(short j=0; j<dpn; j++) sidx_b[i*dpn+j] = idx_b[i]*dpn+j;
 
    PETSC_CHECK(ISCreateGeneral(PETSC_COMM_WORLD, sidx_a.size(), sidx_a.data(), PETSC_COPY_VALUES, &is_a));
    PETSC_CHECK(ISCreateGeneral(PETSC_COMM_WORLD, sidx_b.size(), sidx_b.data(), PETSC_COPY_VALUES, &is_b));
 
    PETSC_CHECK(VecCreateMPI(PETSC_COMM_WORLD, lsize_a*dpn, gsize_a*dpn, &a));
    PETSC_CHECK(VecCreateMPI(PETSC_COMM_WORLD, lsize_b*dpn, gsize_b*dpn, &sc->b_buff));
    PETSC_CHECK(VecSetFromOptions(a));
    PETSC_CHECK(VecSetFromOptions(sc->b_buff));
 
    PETSC_CHECK(VecScatterCreate(a, is_a, sc->b_buff, is_b, &sc->vec_sc));
 
    PETSC_CHECK(VecDestroy(&a));
    PETSC_CHECK(ISDestroy (&is_a));
    PETSC_CHECK(ISDestroy (&is_b));
 
    sc->idx_a.assign(idx_a.begin(), idx_a.end());
    sc->idx_b.assign(idx_b.begin(), idx_b.end());
 
    return sc;
  }
 
  inline scattering* get_scattering(const quadruple<int> spec)
  {
    scattering* ret = NULL;
 
    // scattering must be present
    if(_registry.count(spec))
      ret = _registry[spec];
 
    return ret;
  }
 
  inline void free_scatterings()
  {
    typename hashmap::unordered_map<quadruple<int>,scattering*>::iterator it;
    for(it = _registry.begin(); it != _registry.end(); ++it) {
      if(it->second) {
        delete it->second;
        it->second = NULL;
      }
    }
  }
};
 
 
} // namespace SF
 
 
#endif // _SF_ABSTRACT_VECTOR_H