SF_abstract_matrix.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_MATRIX_H
#define _SF_ABSTRACT_MATRIX_H

#include <mpi.h>

#include "SF_abstract_vector.h"
#include "SF_container.h"
#include "SF_globals.h"
#include "SF_parallel_layout.h"

namespace SF {


template <class T, class S>
class abstract_matrix
{
 protected:
  abstract_matrix() : mesh(NULL),
                      NRows(0),
                      NCols(0),
                      row_dpn(0),
                      col_dpn(0),
                      lsize(0),
                      start(0),
                      stop(0)
  {
  }

  const meshdata<mesh_int_t, mesh_real_t> * mesh;

  int  NRows;         
  int  NCols;         
  int  row_dpn;       
  int  col_dpn;       

  int  lsize;         
  int  start;         
  int  stop;          

  public:

  virtual ~abstract_matrix() = default;

  virtual inline void init(T iNRows, T iNCols, T ilrows, T ilcols,
                           T loc_offset, T mxent)
  {
    this->NRows = iNRows;
    this->NCols = iNCols;
    this->lsize = ilrows;
    this->start = loc_offset;
    this->stop  = this->start + this->lsize;
  }

  inline void init(const meshdata<mesh_int_t, mesh_real_t> & imesh,
                   const T irow_dpn,
                   const T icol_dpn,
                   const T max_edges)
  {
    mesh = &imesh;
    this->row_dpn = irow_dpn;
    this->col_dpn = icol_dpn;

    int rank;
    MPI_Comm_rank(mesh->comm, &rank);

    T M = mesh->pl.num_global_idx() * this->row_dpn;
    T N = mesh->pl.num_global_idx() * this->col_dpn;
    T m = mesh->pl.num_algebraic_idx() * this->row_dpn;
    T n = mesh->pl.num_algebraic_idx() * this->col_dpn;

    T nmax = max_edges;
    if(max_edges < 0) nmax = max_nodal_edgecount(*mesh);

    T loc_offset = mesh->pl.algebraic_layout()[rank]*this->row_dpn;

    init(M, N, m, n, loc_offset, nmax*this->col_dpn);
    zero();
  }

  inline const meshdata<mesh_int_t, mesh_real_t>* mesh_ptr() const
  {
    return this->mesh;
  }

  inline int dpn_row() const { return this->row_dpn; }

  inline int dpn_col() const { return this->col_dpn; }

  virtual void zero() = 0;

  virtual void mult(const abstract_vector<T,S> & x, abstract_vector<T,S> & b) const = 0;

  virtual void mult_LR(const abstract_vector<T, S>& L, const abstract_vector<T, S>& R) = 0;

  virtual void diag_add(const abstract_vector<T, S>& diag) = 0;

  virtual void get_diagonal(abstract_vector<T, S>& vec) const = 0;

  virtual void finish_assembly() = 0;

  virtual void scale(S s) = 0;

  virtual void add_scaled_matrix(const abstract_matrix<T, S>& A, const S s, const bool same_nnz) = 0;

  virtual void duplicate(const abstract_matrix<T,S>& M) = 0;

  virtual void set_values(const vector<T>& row_idx, const vector<T>& col_idx, const vector<S>& vals, bool add) = 0;

  virtual void set_values(const vector<T> & row_idx, const vector<T> & col_idx,
                         const S* vals, bool add) = 0;

  virtual void set_value(T row_idx, T col_idx, S val, bool add) = 0;

  inline void get_values(const vector<T> & row_idx, const vector<T> & col_idx,
                         vector<S> & values) const {
    assert(row_idx.size() == col_idx.size() == values.size());

    for (int i = 0; i < row_idx.size(); i++)
      values[i] = get_value(row_idx[i], col_idx[i]);
  }

  virtual S get_value(SF_int row_idx, SF_int col_idx) const = 0;

  virtual void write(const char* filename) const = 0;


  inline bool equals(const abstract_matrix<T,S> & rhs) const
  {
    if (NRows != rhs.NRows or NCols != rhs.NCols)
      return false;

    bool equal = true;
    for (size_t row = 0; row < NRows; row++)
      for (size_t col = 0; col < NCols; col++)
        if (fabs(get_value(row, col) - rhs.get_value(row, col)) > 0.0001)
          equal = false;

    return equal;
  }

  inline std::string to_string() const
  {
    std::ostringstream stream;
    stream << "matrix (" << this->NRows << " x " << this->NCols << ") [\n";

    for (auto row = 0; row < this->NRows; row++) {
      stream << "   [";
      for (auto col = 0; col < this->NCols; col++) {
        if (col > 0) stream << ", ";
        stream << get_value(row, col);
      }
      stream << "]\n";
    }
    stream << "]";

    return stream.str();
  }
};

} // namespace SF



#endif // _SF_ABSTRACT_MATRIX_H