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 "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(int igsize, int ilsize, int idpn = 1, ltype ilayout = unset) = 0;

  inline std::tuple<int, int> init_common(const meshdata<mesh_int_t, mesh_real_t>& imesh,
                                          int                                      idpn,
                                          ltype                                    inp_layout)
  {
    mesh = &imesh;
    dpn = idpn;
    int N = 0, n = 0;

    switch(inp_layout) {
      case algebraic:
        N = mesh->pl.num_global_idx()*dpn;
        n = mesh->pl.num_algebraic_idx()*dpn;
        layout = algebraic;
        break;

      case nodal:
        N = mesh->l_numpts*dpn;
        MPI_Allreduce(MPI_IN_PLACE, &N, 1, MPI_INT, MPI_SUM, mesh->comm);
        n = mesh->l_numpts*dpn;
        layout = nodal;
        break;

      case elemwise:
        N = mesh->g_numelem;
        n = mesh->l_numelem;
        layout = elemwise;
        break;

      default: break;
    }

    return std::tuple<int, int>(N, n);
  }

  virtual void set(const vector<T>& idx, const vector<S>& vals, const bool additive = false) = 0;

  virtual void set(const vector<T>& idx, const S val) = 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 void release_ptr(S*& p) = 0;

  virtual void const_release_ptr(const S*& p) const = 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);

    int glb_size = this->gsize();
    int 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, "%d\n", glb_size);

      for(int 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++)
      {
        int rsize;
        MPI_Status stat;

        MPI_Recv(&rsize, 1, MPI_INT, pid, SF_MPITAG, comm, &stat);
        wbuff.resize(rsize);
        MPI_Recv(wbuff.data(), rsize*sizeof(S), MPI_BYTE, pid, SF_MPITAG, comm, &stat);

        for(int 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_INT, 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);
    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) VecDestroy(&b_buff);
    if(vec_sc) 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);
  }
};

template<class T>
class scatter_registry
{
  private:
  hashmap::unordered_map<quadruple<int>, scattering*> _registry;

  public:
  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;
    int err = 0;
    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;

    err += ISCreateGeneral(PETSC_COMM_WORLD, idx_a.size(), idx_a.data(), PETSC_COPY_VALUES, &is_a);
    err += ISCreateGeneral(PETSC_COMM_WORLD, idx_b.size(), idx_b.data(), PETSC_COPY_VALUES, &is_b);

    err += ISSetPermutation(is_b);

    Vec a;
    err += VecCreateMPI(PETSC_COMM_WORLD, idx_a.size(), gsize_a*dpn, &a);
    err += VecCreateMPI(PETSC_COMM_WORLD, idx_b.size(), gsize_b*dpn, &sc->b_buff);

    err += VecScatterCreate(a, is_a, sc->b_buff, is_b, &sc->vec_sc);

    err += VecDestroy(&a);
    err += ISDestroy (&is_a);
    err += ISDestroy (&is_b);

    if(err)
      fprintf(stderr, "%s : an error has occurred! Scattering spec: Mesh ID %d, permutation ID %d, numbering %d, dpn %d!\n",
              __func__, int(spec.v1), int(spec.v2), int(spec.v3), int(spec.v4));

    sc->idx_a = idx_a;
    sc->idx_b = idx_b;

    return sc;
  }


  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;
    int err = 0;
    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;

    err += ISCreateGeneral(PETSC_COMM_WORLD, sidx_a.size(), sidx_a.data(), PETSC_COPY_VALUES, &is_a);
    err += ISCreateGeneral(PETSC_COMM_WORLD, sidx_b.size(), sidx_b.data(), PETSC_COPY_VALUES, &is_b);

    err += VecCreateMPI(PETSC_COMM_WORLD, lsize_a*dpn, gsize_a*dpn, &a);
    err += VecCreateMPI(PETSC_COMM_WORLD, lsize_b*dpn, gsize_b*dpn, &sc->b_buff);

    err += VecScatterCreate(a, is_a, sc->b_buff, is_b, &sc->vec_sc);

    err += VecDestroy(&a);
    err += ISDestroy (&is_a);
    err += ISDestroy (&is_b);

    if(err)
      fprintf(stderr, "%s : an error has occurred! Scattering spec: To ID %d, From ID %d, numbering %d, dpn %d!\n",
              __func__, int(spec.v1), int(spec.v2), int(spec.v3), int(spec.v4));

    sc->idx_a = idx_a;
    sc->idx_b = idx_b;

    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