MatrixCollection.cxx

// Copyright (C) 2010 INRIA
// Author(s): Marc Fragu
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your option)
// any later version.
//
// Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.
 
 
#ifndef SELDON_FILE_MATRIX_COLLECTION_CXX
 
#include "MatrixCollection.hxx"
 
 
namespace Seldon
{
 
 
  // MATRIXCOLLECTION //
 
 
  /****************
   * CONSTRUCTORS *
   ****************/
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  MatrixCollection<T, Prop, Storage, Allocator>::MatrixCollection():
    Matrix_Base<T, Allocator>(), Mlocal_(), Mlocal_sum_(1),
    Nlocal_(), Nlocal_sum_(1), matrix_()
  {
    nz_ = 0;
    Mmatrix_ = 0;
    Nmatrix_ = 0;
    Mlocal_sum_.Fill(0);
    Nlocal_sum_.Fill(0);
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  MatrixCollection<T, Prop, Storage, Allocator>
  ::MatrixCollection(int i, int j): Matrix_Base<T, Allocator>(i, j),
                                    Mlocal_(i), Mlocal_sum_(i + 1),
                                    Nlocal_(j), Nlocal_sum_(j + 1),
                                    matrix_(i, j)
  {
    nz_ = 0;
    Mmatrix_ = i;
    Nmatrix_ = j;
    Mlocal_.Fill(0);
    Nlocal_.Fill(0);
    Mlocal_sum_.Fill(0);
    Nlocal_sum_.Fill(0);
  }
  
  
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>::Clear()
  {
    for (int i = 0; i < Mmatrix_; i++)
      for (int j = 0; j < Nmatrix_; j++)
        matrix_(i, j).Nullify();
 
    matrix_.Clear();
 
    nz_ = 0;
    Mmatrix_ = 0;
    Nmatrix_ = 0;
    Mlocal_.Clear();
    Nlocal_.Clear();
    Mlocal_sum_.Clear();
    Nlocal_sum_.Clear();
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>::Nullify()
  {
    for (int i = 0; i < Mmatrix_; i++)
      for (int j = 0; j < Nmatrix_; j++)
        matrix_(i, j).Nullify();
 
    nz_ = 0;
    Mmatrix_ = 0;
    Nmatrix_ = 0;
    Mlocal_.Clear();
    Nlocal_.Clear();
    Mlocal_sum_.Clear();
    Nlocal_sum_.Clear();
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Nullify(int i, int j)
  {
    nz_ -= matrix_(i, j).GetDataSize();
    matrix_(i, j).Nullify();
  }
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>::Deallocate()
  {
    for (int i = 0; i < Mmatrix_; i++)
      for (int j = 0; j < Nmatrix_; j++)
        GetMatrix(i, j).Clear();
    this->Clear();
  }
 
 
  /*********************
   * MEMORY MANAGEMENT *
   *********************/
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Reallocate(int i, int j)
  {
    nz_ = 0;
    Mmatrix_ = i;
    Nmatrix_ = j;
    Mlocal_.Reallocate(i);
    Nlocal_.Reallocate(j);
    Mlocal_sum_.Reallocate(i + 1);
    Nlocal_sum_.Reallocate(j + 1);
    Mlocal_.Fill(0.);
    Nlocal_.Fill(0.);
    Mlocal_sum_.Fill(0.);
    Nlocal_sum_.Fill(0.);
    matrix_.Reallocate(i, j);
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  template <class T0, class Prop0, class Storage0, class Allocator0>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::SetMatrix(int i, int j, const Matrix<T0, Prop0, Storage0, Allocator0>& A)
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= Mmatrix_)
      throw WrongRow("MatrixCollection::SetMatrix()",
                     string("Line index should be in [0, ")
                     + to_str(Mmatrix_ - 1) + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= Nmatrix_)
      throw WrongCol("MatrixCollection::SetMatrix()",
                     string("Column index should be in [0, ")
                     + to_str(Nmatrix_ - 1) + "], but is equal to "
                     + to_str(j) + ".");
    if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
      throw WrongDim("MatrixCollection::SetMatrix()",
                     string("The matrix expected should have ")
                     + to_str(this->Mlocal_(i)) + " rows, but has "
                     + to_str(A.GetM()) + " rows.");
    if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
      throw WrongDim("MatrixCollection::SetMatrix()",
                     string("The matrix expected should have ")
                     + to_str(this->Nlocal_(j)) + " columns, but has "
                     + to_str(A.GetN()) + " columns.");
#endif
 
    nz_ = A.GetDataSize() - matrix_(i, j).GetDataSize();
 
    int Mdiff = A.GetM() - Mlocal_(i);
    int Ndiff = A.GetN() - Nlocal_(j);
 
    Mlocal_(i) = A.GetM();
    Nlocal_(j) = A.GetN();
 
    for (int k = i + 1; k < Mmatrix_ + 1; k++)
      Mlocal_sum_(k) += Mdiff;
 
    for (int k = j + 1; k < Nmatrix_ + 1; k++)
      Nlocal_sum_(k) += Ndiff;
 
    this->m_ = Mlocal_sum_(Mmatrix_);
    this->n_ = Nlocal_sum_(Nmatrix_);
 
    matrix_(i, j).SetData(A.GetM(), A.GetN(), A.GetData());
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  template <class T0, class Prop0, class Allocator0>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::SetMatrix(int i, int j, const Matrix<T0, Prop0, RowSparse, Allocator0>& A)
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= Mmatrix_)
      throw WrongRow("MatrixCollection::SetMatrix()",
                     string("Index should be in [0, ")
                     + to_str(Mmatrix_ - 1) + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= Nmatrix_)
      throw WrongCol("MatrixCollection::SetMatrix()",
                     string("Index should be in [0, ")
                     + to_str(Nmatrix_ - 1) + "], but is equal to "
                     + to_str(j) + ".");
    if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
      throw WrongDim("MatrixCollection::SetMatrix()",
                     string("The matrix expected should have ")
                     + to_str(this->Mlocal_(i)) + " rows, but has "
                     + to_str(A.GetM()) + " rows.");
    if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
      throw WrongDim("MatrixCollection::SetMatrix()",
                     string("The matrix expected should have ")
                     + to_str(this->Nlocal_(j)) + " columns, but has "
                     + to_str(A.GetN()) + " columns.");
#endif
 
    nz_ = A.GetDataSize() - matrix_(i, j).GetDataSize();
 
    int Mdiff = A.GetM() - Mlocal_(i);
    int Ndiff = A.GetN() - Nlocal_(j);
 
    Mlocal_(i) = A.GetM();
    Nlocal_(j) = A.GetN();
 
    for (int k = i + 1; k < Mmatrix_ + 1; k++)
      Mlocal_sum_(k) += Mdiff;
 
    for (int k = j + 1; k < Nmatrix_ + 1; k++)
      Nlocal_sum_(k) += Ndiff;
 
    this->m_ = Mlocal_sum_(Mmatrix_);
    this->n_ = Nlocal_sum_(Nmatrix_);
 
    matrix_(i, j).SetData(A.GetM(), A.GetN(), A.GetNonZeros(), A.GetData(),
                          A.GetPtr(), A.GetInd());
  }
 
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  typename MatrixCollection<T, Prop, Storage, Allocator>
  ::value_type
  MatrixCollection<T, Prop, Storage, Allocator>::operator() (int i,
                                                             int j) const
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->Mlocal_sum_(Mmatrix_))
      throw WrongRow("MatrixCollection::operator(int, int)",
                     string("Row index should be in [0, ")
                     + to_str(this->Mlocal_sum_(Mmatrix_) - 1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->Nlocal_sum_(Nmatrix_))
      throw WrongCol("MatrixCollection::operator(int, int)",
                     string("Column index should be in [0, ")
                     + to_str(this->Nlocal_sum_(Nmatrix_) - 1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
 
    int i_global = 0;
    while (i >= Mlocal_sum_(i_global))
      i_global++;
    i_global--;
 
    int j_global = 0;
    while (j >= Nlocal_sum_(j_global))
      j_global++;
    j_global--;
 
    return matrix_(i_global, j_global)(i - Mlocal_sum_(i_global),
                                       j - Nlocal_sum_(j_global));
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  MatrixCollection<T, Prop, Storage, Allocator>&
  MatrixCollection<T, Prop, Storage, Allocator>
  ::operator= (const MatrixCollection<T, Prop, Storage, Allocator>& A)
  {
    this->Copy(A);
    return *this;
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Copy(const MatrixCollection<T, Prop, Storage, Allocator>& A)
  {
    Clear();
 
    this->nz_ = A.nz_;
    this->m_ = A.GetM();
    this->n_ = A.GetN();
    Mmatrix_ = A.Mmatrix_;
    Nmatrix_ = A.Nmatrix_;
    this->Mlocal_ = A.Mlocal_;
    this->Mlocal_sum_ = A.Mlocal_sum_;
    this->Nlocal_ = A.Nlocal_;
    this->Nlocal_sum_ = A.Nlocal_sum_;
 
    matrix_.Reallocate(Mmatrix_, Nmatrix_);
 
    for (int i = 0; i < Mmatrix_; i++)
      for (int j = 0; j < Nmatrix_; j++)
        SetMatrix(i, j, A.GetMatrix(i, j));
  }
 
  
  /************************
   * CONVENIENT FUNCTIONS *
   ************************/
  
  
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>::Print() const
  {
    for (int i = 0; i < Mmatrix_; i++)
      {
        for (int j = 0; j < Nmatrix_; j++)
          cout << GetMatrix(i, j) << endl;
        cout << endl;
      }
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Print(int i, int j) const
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= Mmatrix_)
      throw WrongRow("MatrixCollection::Print(int, int)",
                     string("Row index should be in [0, ")
                     + to_str(Mmatrix_ - 1) + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= Nmatrix_)
      throw WrongCol("MatrixCollection::Print(int, int)",
                     string("Column index should be in [0, ")
                     + to_str(Nmatrix_ - 1) + "], but is equal to "
                     + to_str(j) + ".");
#endif
 
    cout << matrix_(i, j) << endl;
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Write(string FileName, bool with_size) const
  {
    ofstream FileStream;
    FileStream.open(FileName.c_str(), ofstream::binary);
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("MatrixCollection::Write(string FileName, bool)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->Write(FileStream, with_size);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Write(ostream& FileStream, bool with_size = true) const
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("MatrixCollection::Write(ostream& FileStream, bool)",
                    "The stream is not ready.");
#endif
 
    if (with_size)
      {
        FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&Mmatrix_)),
                         sizeof(int));
        FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&Nmatrix_)),
                         sizeof(int));
      }
 
    int i, j;
    for (i = 0; i < this->GetMmatrix(); i++)
      {
        for (j = 0; j < this->GetNmatrix(); j++)
          GetMatrix(i, j).Write(FileStream, with_size);
      }
 
#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!FileStream.good())
      throw IOError("MatrixCollection::Write(ostream& FileStream, bool)",
                    "Output operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::WriteText(string FileName) const
  {
    ofstream FileStream;
    FileStream.precision(cout.precision());
    FileStream.flags(cout.flags());
    FileStream.open(FileName.c_str());
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("MatrixCollection::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->WriteText(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::WriteText(ostream& FileStream) const
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the file is ready.
    if (!FileStream.good())
      throw IOError("MatrixCollection::WriteText(ostream& FileStream)",
                    "The stream is not ready.");
#endif
 
    int i, j;
    for (i = 0; i < this->GetMmatrix(); i++)
      {
        for (j = 0; j < this->GetNmatrix(); j++)
          GetMatrix(i, j).WriteText(FileStream);
        FileStream << endl;
      }
 
#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!FileStream.good())
      throw IOError("MatrixCollection::WriteText(ostream& FileStream)",
                    "Output operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>::Read(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str(), ifstream::binary);
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("MatrixCollection::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->Read(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  void MatrixCollection<T, Prop, Storage, Allocator>
  ::Read(istream& FileStream)
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("MatrixCollection::Read(istream& FileStream)",
                    "The stream is not ready.");
#endif
 
    int *new_m, *new_n;
    new_m = new int;
    new_n = new int;
 
    FileStream.read(reinterpret_cast<char*>(new_m), sizeof(int));
    FileStream.read(reinterpret_cast<char*>(new_n), sizeof(int));
 
    this->Reallocate(*new_m, *new_n);
 
    T working_matrix;
    int i, j;
    for (i = 0; i < *new_m; i++)
      for (j = 0; j < *new_n; j++)
        {
          working_matrix.Read(FileStream);
          SetMatrix(i, j, working_matrix);
          working_matrix.Nullify();
        }
 
 
    delete new_n;
    delete new_m;
 
#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!FileStream.good())
      throw IOError("MatrixCollection::Read(istream& FileStream)",
                    "Input operation failed.");
#endif
 
  }
  
  
} // namespace Seldon.
 
 
#define SELDON_FILE_MATRIX_COLLECTION_CXX
#endif