belosLinProb.hh

Go to the documentation of this file.

 
#ifndef BELOSLINEARPROBLEM_HH_
#define BELOSLINEARPROBLEM_HH_
 
#include "basics.hh"
#include "sparseMatrix.hh"
#include "toolbox/sequence.hh"
#include <BelosLinearProblem.hpp>
#include <Tpetra_CrsMatrix.hpp>
 
namespace concepts {
 
  template<class T, class MV = Tpetra::MultiVector<T, int>,
      class OP = Tpetra::Operator<T> >
  class BelosLinProb: public Belos::LinearProblem<T, MV, OP> {
  public:
 
    BelosLinProb() {
    }
    ;
 
    virtual ~BelosLinProb() {
    }
    ;
 
    BelosLinProb(Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    BelosLinProb(concepts::SparseMatrix<T>& sparse,
        Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    void setConceptsRHS(const Vector<T>& rhs,
        Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    void setConceptsRHS(Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    void writeSolution(Vector<T>& vec,
        Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    void writeSolution(Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    Teuchos::RCP<Tpetra::CrsMatrix<T, int, int> > getCrsMat() {
      return A_;
    }
    ;
    void setCrsMat(Teuchos::RCP<Tpetra::CrsMatrix<T, int, int> > A) {
      A = A_;
      this->setOperator(A_);
    }
 
  private:
 
    Teuchos::RCP<Tpetra::CrsMatrix<T, int, int> > A_;
  };
 
  template<class T, class MV, class OP>
  BelosLinProb<T, MV, OP>::BelosLinProb(SparseMatrix<T> &sparse,
                                        Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
 
    //Thread information
    int MyPID   = Comm->getRank();
    int NumProc = Comm->getSize();
 
    //local matrix data recieved from the assembler thread
    int * nnzPerRow;
    int * col;
    T * val;
 
    //iterators
    int * nnzPerRowIter;
    int * colIter;
    T * valIter;
 
    //Global matrix information
    int globalSize = 0;
    int max_ = 0;
 
    //local matrix information
    int localSize = 0;
    int localNnz = 0;
 
    /*
     ===================================
     Extract data for other processes and send it.
     Just one thread is allowed to know the sparse matrix
     ===================================
     */
    globalSize = sparse.nofRows();
 
    //broadcast globalMatrix Informations
    Comm->broadcast(0, sizeof(int), (char*) (&globalSize));
 
    // compute local size for each processor
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //number of rows of this thread
    int bias = localSizes[MyPID];
    typename SparseMatrix<T>::iterator iter = sparse.begin(0);
 
    for (int i = 1; i < NumProc; ++i) {
      //local Size of the process we collect data
      localSize = localSizes[i];
      //memory for local entries per row array
      nnzPerRow = new int[localSize];
 
      //iterator
      nnzPerRowIter = nnzPerRow;
 
      DEBUGL(0, "localSize = " << localSize);
      DEBUGL(0, "bias      = " << bias);
 
      //count numbers of nonzeros per row
      for (uint j = 0; j < (uint) localSize; ++j) {
        *nnzPerRowIter = 0;
        iter = sparse.begin(j + bias);
        while (iter != sparse.end() && (iter++).row() == (j + bias))
          ++(*nnzPerRowIter);
        ++nnzPerRowIter;
      } //for over rows
 
      //send data to thread number i and delete it
      Comm->send(localSize * sizeof(int), (char*) nnzPerRow, i);
 
      localNnz = 0;
      nnzPerRowIter = nnzPerRow;
      for (int j = 0; j < localSize; ++j)
        localNnz += *nnzPerRowIter++;
 
      delete[] nnzPerRow;
 
      //memory for global values and colums
      val = new T[localNnz];
      col = new int[localNnz];
 
      //extract global data
      valIter = val;
      colIter = col;
      for (uint j = 0; j < (uint) localSize; ++j) {
        iter = sparse.begin(j + bias);
        while (iter.row() == (j + bias) && (iter != sparse.end())) {
          *colIter++ = iter.col();
          *valIter++ = *iter++;
        } //while
      } //inner for
 
      //send data to thread number i and delete it
      Comm->send(localNnz * sizeof(T), (char*) val, i);
      Comm->send(localNnz * sizeof(int), (char*) col, i);
      delete[] val;
      delete[] col;
 
      //increase bias
      bias += localSize;
    }      //for that sends data to other processes
 
    /*
     ===================================
     Extract own data
     ===================================
     */
    //local Size of the process we collect data
    localSize = localSizes[0];
 
    //memory for local entrie per row array
    nnzPerRow = new int[localSize];
 
    //iterator
    nnzPerRowIter = nnzPerRow;
    iter = sparse.begin(0);
    //TODO uint -> int Kontrollieren!
    //count numbers of nonzeros per row
    for (uint i = 0; i < (uint) localSize; ++i) {
      *nnzPerRowIter = 0;
      iter = sparse.begin(i);
      while (iter != sparse.end() && (iter++).row() == i)
        ++(*nnzPerRowIter);
      ++nnzPerRowIter;
    } //for over rows
 
    localNnz = 0;
    nnzPerRowIter = nnzPerRow;
    for (int i = 0; i < localSize; ++i) {
      max_ = std::max(max_, *nnzPerRowIter);
      localNnz += *nnzPerRowIter++;
    }
    //memory for global values and colums
    val = new T[localNnz];
    col = new int[localNnz];
 
    //extract global data
    valIter = val;
    colIter = col;
    for (uint i = 0; i < (uint) localSize; ++i) {
      iter = sparse.begin(i);
      while (iter.row() == i && (iter != sparse.end())) {
        *colIter++ = iter.col();
        *valIter++ = *iter++;
      } //while
    } //inner for
 
    /*
     ===============================
     fill matrix
     ===============================
     */
    //Create continouse and uniform mapping (in terms of rows)
    Teuchos::RCP<const Tpetra::Map<int, int> > mapStiffRow =
        Tpetra::createContigMap<int, int>(globalSize, localSize, Comm);
 
    //declare matrix, rhs and solutionVector
    A_ = Teuchos::RCP<Tpetra::CrsMatrix<T, int, int> >(
        new Tpetra::CrsMatrix<T, int, int>(mapStiffRow, max_));
 
    Teuchos::RCP<Tpetra::MultiVector<T, int> > sol(
        new Tpetra::MultiVector<T, int>(mapStiffRow, 1));
    Teuchos::RCP<Tpetra::MultiVector<T, int> > rhs2(
        new Tpetra::MultiVector<T, int>(mapStiffRow, 1));
 
    //iterators
    nnzPerRowIter = nnzPerRow;
    valIter = val;
    colIter = col;
    int globalIndex = 0;
    int currNnz = 0;
 
    for (int k = 0; k < localSize; ++k) {
      globalIndex = mapStiffRow->getGlobalElement(k);
      currNnz = *nnzPerRowIter;
      if (currNnz > 0) {
        A_->insertGlobalValues(globalIndex,
            Teuchos::ArrayView<int>(colIter, currNnz),
            Teuchos::ArrayView<T>(valIter, currNnz));
        colIter += currNnz;
        valIter += currNnz;
      } //nnz != 0
      nnzPerRowIter++;
    }
 
    A_->fillComplete();
 
    delete[] nnzPerRow;
    delete[] val;
    delete[] col;
 
    this->setOperator(A_);
  }
  ;
 
  template<class T, class MV, class OP>
  BelosLinProb<T, MV, OP>::BelosLinProb(
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
 
    //Thread information
    int MyPID = Comm->getRank();
    int NumProc = Comm->getSize();
 
    //This thread is not the first
    if (MyPID == 0)
      throw std::logic_error("The first thread needs a Matrix");
 
    //local matrix data recieved from the assembler thread
    int * nnzPerRow;
    int * col;
    T * val;
 
    //iterators
    int * nnzPerRowIter;
    int * colIter;
    T * valIter;
 
    //Global matrix information
    int globalSize = 0;
    int max_ = 0;
 
    //local matrix information
    int localSize = 0;
    int localNnz = 0;
    //broadcast globalMatrix Informations
    Comm->broadcast(0, sizeof(int), (char*) (&globalSize));
 
    // compute local data sizes
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //calculate local data size
    localSize = localSizes[MyPID];
 
    //allocate memory for array of number of entries and recieve data from assambler
    nnzPerRow = new int[localSize];
    nnzPerRowIter = nnzPerRow;
    Comm->receive(0, localSize * sizeof(int), (char*) nnzPerRow);
    max_ = *nnzPerRowIter;
    for (int i = 0; i < localSize; ++i) {
      max_ = std::max(max_, *nnzPerRowIter);
      localNnz += *nnzPerRowIter++;
    }
    //allocate memory for locale data
    col = new int[localNnz];
    val = new T[localNnz];
 
    //receive data from assambler thread
    Comm->receive(0, localNnz * sizeof(T), (char*) val);
    Comm->receive(0, localNnz * sizeof(int), (char*) col);
 
    /*
     ===============================
     fill matrix
     ===============================
     */
    //Create continouse and uniform mapping (in terms of rows)
    Teuchos::RCP<const Tpetra::Map<int, int> > mapStiffRow =
        Tpetra::createContigMap<int, int>(globalSize, localSize, Comm);
 
    //declare matrix, rhs and solutionVector
    A_ = Teuchos::RCP<Tpetra::CrsMatrix<T, int, int> >(
        new Tpetra::CrsMatrix<T, int, int>(mapStiffRow, max_));
 
    //iterators
    nnzPerRowIter = nnzPerRow;
    valIter = val;
    colIter = col;
    int globalIndex = 0;
    int currNnz = 0;
 
    for (int k = 0; k < localSize; ++k) {
      globalIndex = mapStiffRow->getGlobalElement(k);
      currNnz = *nnzPerRowIter;
      if (currNnz > 0) {
        A_->insertGlobalValues(globalIndex,
            Teuchos::ArrayView<int>(colIter, currNnz),
            Teuchos::ArrayView<T>(valIter, currNnz));
        colIter += currNnz;
        valIter += currNnz;
      } //nnz != 0
      nnzPerRowIter++;
    }
 
    A_->fillComplete();
 
    delete[] nnzPerRow;
    delete[] val;
    delete[] col;
 
    this->setOperator(A_);
  }
  ;
 
  template<class T, class MV, class OP>
  void BelosLinProb<T, MV, OP>::setConceptsRHS(const Vector<T> &vec,
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
 
    //Thread information
    int MyPID = Comm->getRank();
    int NumProc = Comm->getSize();
 
    T * rhsVal;
    T * rhsValIter;
    T * vecIter;
    int globalSize = vec.size();
 
    //broadcast globalMatrix Informations
    Comm->broadcast(0, sizeof(int), (char*) (&globalSize));
 
    // compute local data sizes
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //number of rows of this thread 
    int bias = localSizes[MyPID];
 
    for (int i = 1; i < NumProc; ++i) {
      int localSize = localSizes[i];;
      rhsVal = new T[localSize];
      rhsValIter = rhsVal;
      vecIter = &(vec[bias]);
      for (int j = 0; j < localSize; ++j)
        *rhsValIter++ = *vecIter++;
      //send and delete data for other processes
      Comm->send(localSize * sizeof(T), (char*) rhsVal, i);
      delete[] rhsVal;
 
      //increase bias
      bias += localSize;
    }
    //my local Size
    int localSize = localSizes[MyPID];
    rhsVal = new T[localSize];
    rhsValIter = rhsVal;
    vecIter = &(vec[0]);
    for (int j = 0; j < localSize; ++j)
      *rhsValIter++ = *vecIter++;
 
    Teuchos::RCP<Tpetra::MultiVector<T, int> > rhs2(
        new Tpetra::MultiVector<T, int>(A_->getRowMap(), 1));
 
    Teuchos::RCP<Tpetra::MultiVector<T, int> > sol(
        new Tpetra::MultiVector<T, int>(A_->getRowMap(), 1));
 
    rhsValIter = rhsVal;
    for (int k = 0; k < localSize; ++k) {
      sol->replaceGlobalValue((A_->getRowMap())->getGlobalElement(k), 0,
          *rhsValIter);
      rhs2->replaceGlobalValue((A_->getRowMap())->getGlobalElement(k), 0,
          *rhsValIter++);
    }
 
    delete[] rhsVal;
 
    this->setRHS(rhs2);
    this->setLHS(sol);
 
  }
 
  template<class T, class MV, class OP>
  void BelosLinProb<T, MV, OP>::setConceptsRHS(
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
 
    //Thread information
    int MyPID = Comm->getRank();
    int NumProc = Comm->getSize();
 
    T * rhsVal;
    T * rhsValIter;
    int globalSize;
 
    //broadcast globalMatrix Informations
    Comm->broadcast(0, sizeof(int), (char*) (&globalSize));
 
    // compute local data sizes
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //calculate local data size
    int localSize = localSizes[MyPID];
 
    rhsVal = new T[localSize];
    Comm->receive(0, localSize * sizeof(T), (char*) rhsVal);
 
    Teuchos::RCP<Tpetra::MultiVector<T, int> > rhs2(
        new Tpetra::MultiVector<T, int>(A_->getRowMap(), 1));
 
    Teuchos::RCP<Tpetra::MultiVector<T, int> > sol(
        new Tpetra::MultiVector<T, int>(A_->getRowMap(), 1));
 
    rhsValIter = rhsVal;
    for (int k = 0; k < localSize; ++k) {
      sol->replaceGlobalValue((A_->getRowMap())->getGlobalElement(k), 0,
          *rhsValIter);
      rhs2->replaceGlobalValue((A_->getRowMap())->getGlobalElement(k), 0,
          *rhsValIter++);
    }
 
    delete[] rhsVal;
 
    this->setRHS(rhs2);
    this->setLHS(sol);
    //solve
  }
 
  template<class T, class MV, class OP>
  void BelosLinProb<T, MV, OP>::writeSolution(Vector<T> & sol,
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
 
    int MyPID = Comm->getRank();
    int NumProc = Comm->getSize();
    int globalSize = sol.size();
 
    // compute local data sizes
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //calculate local data size
    int localSize = localSizes[MyPID];
 
    T *solutionData = new T[globalSize];
    memcpy(solutionData, (this->getLHS()->get1dViewNonConst()).get(),
        sizeof(T) * localSize);
 
    //collects data of all processes
    int ith_size = 0;
    int bias = localSize;
    for (int i = 1; i < NumProc; ++i) {
      ith_size = localSizes[i];
      Comm->receive(i, sizeof(T) * ith_size,
          (char*) (solutionData + bias));
      bias += ith_size;
    }
    T * solutionDataIter = solutionData;
    for (int i = 0; i < sol.size(); ++i) {
      sol[i] = *solutionDataIter++;
    }
    delete[] solutionData;
  }
 
  template<class T, class MV, class OP>
  void BelosLinProb<T, MV, OP>::writeSolution(
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) {
    //sends the part of the solution to the root thread (with PID = 0)
    int MyPID = Comm->getRank();
    int NumProc = Comm->getSize();
    int globalSize = A_->getRowMap()->getGlobalNumElements();
 
    // compute local data sizes
    concepts::Sequence<int> localSizes(NumProc);
    int restsize = globalSize;
    for(int i = 0; i < NumProc; ++i)
      restsize    -= (localSizes[i] = restsize / (NumProc - i));
    DEBUGL(0, "localSizes = " << localSizes);
 
    //calculate local data size
    int localSize = localSizes[MyPID];
 
    Comm->send(sizeof(T) * localSize,
        (const char*) (this->getLHS()->get1dViewNonConst()).get(), 0);
  }
 
} /* namespace concepts */
#endif /* BELOSLINEARPROBLEM_HH_ */