paralution_dealii.hpp

Go to the documentation of this file.

// Deal II path 
#include <lac/sparse_matrix.h>
#include <lac/sparsity_pattern.h>
#include <lac/vector.h>
#include <lac/block_sparse_matrix.h>


// PARALUTION path
#include <paralution.hpp>

template <typename ValueType>
void import_dealii_matrix(dealii::BlockSparseMatrix<ValueType> &deal_mat,
                          paralution::LocalMatrix<ValueType> *mat) {

  // TODO
  // change for a Base class

  // Matrix size
  int n = deal_mat.n();
  int m = deal_mat.m();
  //  int nnz = deal_mat.n_actually_nonzero_elements();
  int nnz = deal_mat.n_nonzero_elements();

  // COO values
  int *col = NULL;
  int *row = NULL;
  ValueType *val = NULL;

  paralution::allocate_host(nnz, &row);
  paralution::allocate_host(nnz, &col);
  paralution::allocate_host(nnz, &val);

  dealii::BlockSparseMatrix<double>::iterator it = deal_mat.begin(),
    it_end = deal_mat.end();

  int i = 0;
  for (; it!=it_end; ++it) {

    col[i] = it.operator*().column();
    row[i] = it.operator*().row();
    val[i] = it.operator*().value();

    ++i;
  }

  assert(i == nnz);

  mat->SetDataPtrCOO(&row, &col, &val,
                     "Deal II Matrix",
                     nnz, n, m);
  
  mat->ConvertToCSR();

} ;


template <typename ValueType>
void import_dealii_matrix(const dealii::SparsityPattern &sp, 
                          const dealii::SparseMatrix<ValueType> &deal_mat,
                          paralution::LocalMatrix<ValueType> *mat) {

  // Matrix size
  int n = deal_mat.n();
  int m = deal_mat.m();
  //  int nnz = deal_mat.n_actually_nonzero_elements();
  int nnz = deal_mat.n_nonzero_elements();
  
  // COO values
  int *col = NULL;
  int *row = NULL;
  ValueType *val = NULL;

  paralution::allocate_host(nnz, &row);
  paralution::allocate_host(nnz, &col);
  paralution::allocate_host(nnz, &val);

  dealii::SparsityPatternIterators::Iterator it = sp.begin();
  dealii::SparsityPatternIterators::Iterator it_end = sp.end();

  int i = 0;
  for (; it!=it_end; ++it) {

    col[i] = it->column();
    row[i] = it->row();
    val[i] = deal_mat.el(it->row(),it->column());

    ++i;    
  }

  assert(i == nnz);

  mat->SetDataPtrCOO(&row, &col, &val,
                     "Deal II Matrix",
                     nnz, n, m);
  
  // mat->ConvertToCSR();

} ;

template <typename ValueType>
void import_dealii_vector(const dealii::Vector<ValueType> &deal_vec,                             
                             paralution::LocalVector<ValueType> *vec) {

  vec->MoveToHost();

  ValueType *val = NULL;
  paralution::allocate_host(int(deal_vec.size()), &val);

  for (unsigned int i=0; i<deal_vec.size(); ++i)
    val[i] = deal_vec[i];

  vec->SetDataPtr(&val, "DealII vector", int(deal_vec.size()));


  //  assert(int(vec->get_size()) == int(deal_vec.size()));
  //
  //  This is slower
  //
  //  for (unsigned int i=0; i<deal_vec.size(); ++i)
  //    (*vec)[i] = deal_vec[i];

}; 

template <typename ValueType>
void export_dealii_vector(paralution::LocalVector<ValueType> &vec,
                            dealii::Vector<ValueType> *deal_vec) {

  vec.MoveToHost();
  assert(int(vec.get_size()) == int(deal_vec->size()));

  int size = vec.get_size();
  ValueType *val = NULL;

  vec.LeaveDataPtr(&val);

  for (int i=0; i<size; ++i)
    (*deal_vec)[i] = val[i];

  paralution::free_host(&val);
  
  //
  //  This is slower
  //
  //  for (unsigned int i=0; i<deal_vec->size(); ++i)
  //     (*deal_vec)[i] = vec[i];

}; 

template <typename ValueType>
void import_dealii_vector(const dealii::BlockVector<ValueType> &deal_vec,                             
                             paralution::LocalVector<ValueType> *vec) {

  vec->MoveToHost();

  ValueType *val = NULL;
  paralution::allocate_host(int(deal_vec.size()), &val);

  for (unsigned int i=0; i<deal_vec.size(); ++i)
    val[i] = deal_vec[i];

  vec->SetDataPtr(&val, "DealII vector", int(deal_vec.size()));


  //  assert(int(vec->get_size()) == int(deal_vec.size()));
  //
  //  This is slower
  //
  //  for (unsigned int i=0; i<deal_vec.size(); ++i)
  //    (*vec)[i] = deal_vec[i];

}; 

template <typename ValueType>
void export_dealii_vector(paralution::LocalVector<ValueType> &vec,
                            dealii::BlockVector<ValueType> *deal_vec) {

  vec.MoveToHost();
  assert(int(vec.get_size()) == int(deal_vec->size()));

  int size = vec.get_size();
  ValueType *val = NULL;

  vec.LeaveDataPtr(&val);

  for (int i=0; i<size; ++i)
    (*deal_vec)[i] = val[i];

  paralution::free_host(&val);
  
  //
  //  This is slower
  //
  //  for (unsigned int i=0; i<deal_vec->size(); ++i)
  //     (*deal_vec)[i] = vec[i];

};