matrix_sparse/Matrix_Sparse.cxx

// Copyright (C) 2001-2009 Vivien Mallet
// Copyright (C) 2003-2009 Marc Duruflé
//
// 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_SPARSE_CXX

#include "Matrix_Sparse.hxx"

#include "../vector/Functions_Arrays.cxx"

#include <set>


namespace Seldon
{


  /****************
   * CONSTRUCTORS *
   ****************/



  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::Matrix_Sparse():
    Matrix_Base<T, Allocator>()
  {
    nz_ = 0;
    ptr_ = NULL;
    ind_ = NULL;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::Matrix_Sparse(int i,
                                                                   int j):
    Matrix_Base<T, Allocator>(i, j)
  {
    nz_ = 0;
    ptr_ = NULL;
    ind_ = NULL;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::
  Matrix_Sparse(int i, int j, int nz):
    Matrix_Base<T, Allocator>(i, j)
  {

    this->nz_ = nz;

#ifdef SELDON_CHECK_DIMENSIONS
    if (nz_ < 0)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                       "Invalid number of non-zero elements: " + to_str(nz)
                       + ".");
      }
    if (nz_ > 0
        && (j == 0
            || static_cast<long int>(nz_-1) / static_cast<long int>(j)
            >= static_cast<long int>(i)))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                       string("There are more values (") + to_str(nz)
                       + " values) than elements in the matrix ("
                       + to_str(i) + " by " + to_str(j) + ").");
      }
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
                                              sizeof(int)) );

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ptr_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ptr_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ")
                     + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1) )
                     + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
                     + " row or column start indices, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ind_ = reinterpret_cast<int*>( calloc(nz_, sizeof(int)) );

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        free(ptr_);
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ind_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        free(ptr_);
        ptr_ = NULL;
        this->data_ = NULL;
      }
    if (ind_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ") + to_str(sizeof(int) * nz)
                     + " bytes to store " + to_str(nz)
                     + " row or column indices, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        this->data_ = this->allocator_.allocate(nz_, this);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        free(ptr_);
        ptr_ = NULL;
        free(ind_);
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (this->data_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        free(ptr_);
        ptr_ = NULL;
        free(ind_);
        ind_ = NULL;
      }
    if (this->data_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ") + to_str(sizeof(int) * nz)
                     + " bytes to store " + to_str(nz) + " values, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

  }



  template <class T, class Prop, class Storage, class Allocator>
  template <class Storage0, class Allocator0,
            class Storage1, class Allocator1,
            class Storage2, class Allocator2>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::
  Matrix_Sparse(int i, int j,
                Vector<T, Storage0, Allocator0>& values,
                Vector<int, Storage1, Allocator1>& ptr,
                Vector<int, Storage2, Allocator2>& ind):
    Matrix_Base<T, Allocator>(i, j)
  {
    nz_ = values.GetLength();

#ifdef SELDON_CHECK_DIMENSIONS
    // Checks whether vector sizes are acceptable.

    if (ind.GetLength() != nz_)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::Matrix_Sparse(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("There are ") + to_str(nz_) + " values but "
                       + to_str(ind.GetLength()) + " row or column indices.");
      }

    if (ptr.GetLength()-1 != Storage::GetFirst(i, j))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::Matrix_Sparse(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("The vector of start indices contains ")
                       + to_str(ptr.GetLength()-1) + string(" row or column")
                       + string(" start indices (plus the number")
                       + " of non-zero entries) but there are "
                       + to_str(Storage::GetFirst(i, j))
                       + " rows or columns (" + to_str(i) + " by "
                       + to_str(j) + " matrix).");
      }

    if (nz_ > 0
        && (j == 0
            || static_cast<long int>(nz_-1) / static_cast<long int>(j)
            >= static_cast<long int>(i)))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::Matrix_Sparse(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("There are more values (")
                       + to_str(values.GetLength())
                       + " values) than elements in the matrix ("
                       + to_str(i) + " by " + to_str(j) + ").");
      }
#endif

    this->ptr_ = ptr.GetData();
    this->ind_ = ind.GetData();
    this->data_ = values.GetData();

    ptr.Nullify();
    ind.Nullify();
    values.Nullify();
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::
  Matrix_Sparse(const Matrix_Sparse<T, Prop, Storage, Allocator>& A)
  {
    this->m_ = 0;
    this->n_ = 0;
    this->nz_ = 0;
    ptr_ = NULL;
    ind_ = NULL;
    this->Copy(A);
  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>::~Matrix_Sparse()
  {
    this->m_ = 0;
    this->n_ = 0;

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        if (ptr_ != NULL)
          {
            free(ptr_);
            ptr_ = NULL;
          }

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        ptr_ = NULL;
      }
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        if (ind_ != NULL)
          {
            free(ind_);
            ind_ = NULL;
          }

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        ind_ = NULL;
      }
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        if (this->data_ != NULL)
          {
            this->allocator_.deallocate(this->data_, nz_);
            this->data_ = NULL;
          }

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->nz_ = 0;
        this->data_ = NULL;
      }
#endif

    this->nz_ = 0;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_Sparse<T, Prop, Storage, Allocator>::Clear()
  {
    this->~Matrix_Sparse();
  }


  /*********************
   * MEMORY MANAGEMENT *
   *********************/



  template <class T, class Prop, class Storage, class Allocator>
  template <class Storage0, class Allocator0,
            class Storage1, class Allocator1,
            class Storage2, class Allocator2>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::
  SetData(int i, int j,
          Vector<T, Storage0, Allocator0>& values,
          Vector<int, Storage1, Allocator1>& ptr,
          Vector<int, Storage2, Allocator2>& ind)
  {
    this->Clear();
    this->m_ = i;
    this->n_ = j;
    this->nz_ = values.GetLength();

#ifdef SELDON_CHECK_DIMENSIONS
    // Checks whether vector sizes are acceptable.

    if (ind.GetLength() != nz_)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::SetData(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("There are ") + to_str(nz_) + " values but "
                       + to_str(ind.GetLength()) + " row or column indices.");
      }

    if (ptr.GetLength()-1 != Storage::GetFirst(i, j))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::SetData(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("The vector of start indices contains ")
                       + to_str(ptr.GetLength()-1) + string(" row or column")
                       + string(" start indices (plus the number")
                       + " of non-zero entries) but there are "
                       + to_str(Storage::GetFirst(i, j))
                       + " rows or columns (" + to_str(i) + " by "
                       + to_str(j) + " matrix).");
      }

    if (nz_ > 0
        && (j == 0
            || static_cast<long int>(nz_-1) / static_cast<long int>(j)
            >= static_cast<long int>(i)))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_Sparse::SetData(int, int, ")
                       + "const Vector&, const Vector&, const Vector&)",
                       string("There are more values (")
                       + to_str(values.GetLength())
                       + " values) than elements in the matrix ("
                       + to_str(i) + " by " + to_str(j) + ").");
      }
#endif

    this->ptr_ = ptr.GetData();
    this->ind_ = ind.GetData();
    this->data_ = values.GetData();

    ptr.Nullify();
    ind.Nullify();
    values.Nullify();
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_Sparse<T, Prop, Storage, Allocator>
  ::SetData(int i, int j, int nz,
            typename Matrix_Sparse<T, Prop, Storage, Allocator>
            ::pointer values,
            int* ptr, int* ind)
  {
    this->Clear();

    this->m_ = i;
    this->n_ = j;

    this->nz_ = nz;

    this->data_ = values;
    ind_ = ind;
    ptr_ = ptr;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_Sparse<T, Prop, Storage, Allocator>::Nullify()
  {
    this->data_ = NULL;
    this->m_ = 0;
    this->n_ = 0;
    nz_ = 0;
    ptr_ = NULL;
    ind_ = NULL;
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_Sparse<T, Prop, Storage, Allocator>::
  Copy(const Matrix_Sparse<T, Prop, Storage, Allocator>& A)
  {
    this->Clear();
    int nz = A.nz_;
    int i = A.m_;
    int j = A.n_;
    this->nz_ = nz;
    this->m_ = i;
    this->n_ = j;
    if ((i == 0)||(j == 0))
      {
        this->m_ = 0;
        this->n_ = 0;
        this->nz_ = 0;
        return;
      }

#ifdef SELDON_CHECK_DIMENSIONS
    if (nz_ > 0
        && (j == 0
            || static_cast<long int>(nz_-1) / static_cast<long int>(j)
            >= static_cast<long int>(i)))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                       string("There are more values (") + to_str(nz)
                       + " values) than elements in the matrix ("
                       + to_str(i) + " by " + to_str(j) + ").");
      }
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
                                              sizeof(int)) );
        memcpy(this->ptr_, A.ptr_,
               (Storage::GetFirst(i, j) + 1) * sizeof(int));
#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ptr_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ptr_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ")
                     + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1) )
                     + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
                     + " row or column start indices, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ind_ = reinterpret_cast<int*>( calloc(nz_, sizeof(int)) );
        memcpy(this->ind_, A.ind_, nz_ * sizeof(int));

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        free(ptr_);
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (ind_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        free(ptr_);
        ptr_ = NULL;
        this->data_ = NULL;
      }
    if (ind_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ") + to_str(sizeof(int) * nz)
                     + " bytes to store " + to_str(nz)
                     + " row or column indices, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        this->data_ = this->allocator_.allocate(nz_, this);
        this->allocator_.memorycpy(this->data_, A.data_, nz_);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->n_ = 0;
        free(ptr_);
        ptr_ = NULL;
        free(ind_);
        ind_ = NULL;
        this->data_ = NULL;
      }
    if (this->data_ == NULL)
      {
        this->m_ = 0;
        this->n_ = 0;
        free(ptr_);
        ptr_ = NULL;
        free(ind_);
        ind_ = NULL;
      }
    if (this->data_ == NULL && i != 0 && j != 0)
      throw NoMemory("Matrix_Sparse::Matrix_Sparse(int, int, int)",
                     string("Unable to allocate ") + to_str(sizeof(int) * nz)
                     + " bytes to store " + to_str(nz) + " values, for a "
                     + to_str(i) + " by " + to_str(j) + " matrix.");
#endif

  }


  /*******************
   * BASIC FUNCTIONS *
   *******************/



  template <class T, class Prop, class Storage, class Allocator>
  int Matrix_Sparse<T, Prop, Storage, Allocator>::GetNonZeros() const
  {
    return nz_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  int Matrix_Sparse<T, Prop, Storage, Allocator>::GetDataSize() const
  {
    return nz_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  int* Matrix_Sparse<T, Prop, Storage, Allocator>::GetPtr() const
  {
    return ptr_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  int* Matrix_Sparse<T, Prop, Storage, Allocator>::GetInd() const
  {
    return ind_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  int Matrix_Sparse<T, Prop, Storage, Allocator>::GetPtrSize() const
  {
    return (Storage::GetFirst(this->m_, this->n_) + 1);
  }



  template <class T, class Prop, class Storage, class Allocator>
  int Matrix_Sparse<T, Prop, Storage, Allocator>::GetIndSize() const
  {
    return nz_;
  }


  /**********************************
   * ELEMENT ACCESS AND AFFECTATION *
   **********************************/



  template <class T, class Prop, class Storage, class Allocator>
  inline typename Matrix_Sparse<T, Prop, Storage, Allocator>::value_type
  Matrix_Sparse<T, Prop, Storage, Allocator>::operator() (int i, int j) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_Sparse::operator()",
                     string("Index should be in [0, ") + to_str(this->m_-1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_Sparse::operator()",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif

    int k,l;
    int a,b;

    a = ptr_[Storage::GetFirst(i, j)];
    b = ptr_[Storage::GetFirst(i, j) + 1];

    if (a == b)
      return T(0);

    l = Storage::GetSecond(i, j);

    for (k = a; (k < b-1) && (ind_[k] < l); k++);

    if (ind_[k] == l)
      return this->data_[k];
    else
      return T(0);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline typename Matrix_Sparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_Sparse<T, Prop, Storage, Allocator>::Val(int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_Sparse::Val(int, int)",
                     string("Index should be in [0, ") + to_str(this->m_-1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_Sparse::Val(int, int)",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif

    int k, l;
    int a, b;

    a = ptr_[Storage::GetFirst(i, j)];
    b = ptr_[Storage::GetFirst(i, j) + 1];

    if (a == b)
      throw WrongArgument("Matrix_Sparse::Val(int, int)",
                          "No reference to element (" + to_str(i) + ", "
                          + to_str(j)
                          + ") can be returned: it is a zero entry.");

    l = Storage::GetSecond(i, j);

    for (k = a; (k < b-1) && (ind_[k] < l); k++);

    if (ind_[k] == l)
      return this->data_[k];
    else
      throw WrongArgument("Matrix_Sparse::Val(int, int)",
                          "No reference to element (" + to_str(i) + ", "
                          + to_str(j)
                          + ") can be returned: it is a zero entry.");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline
  const typename Matrix_Sparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_Sparse<T, Prop, Storage, Allocator>::Val(int i, int j) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_Sparse::Val(int, int)",
                     string("Index should be in [0, ") + to_str(this->m_-1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_Sparse::Val(int, int)",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif

    int k, l;
    int a, b;

    a = ptr_[Storage::GetFirst(i, j)];
    b = ptr_[Storage::GetFirst(i, j) + 1];

    if (a == b)
      throw WrongArgument("Matrix_Sparse::Val(int, int)",
                          "No reference to element (" + to_str(i) + ", "
                          + to_str(j)
                          + ") can be returned: it is a zero entry.");

    l = Storage::GetSecond(i, j);

    for (k = a; (k < b-1) && (ind_[k] < l); k++);

    if (ind_[k] == l)
      return this->data_[k];
    else
      throw WrongArgument("Matrix_Sparse::Val(int, int)",
                          "No reference to element (" + to_str(i) + ", "
                          + to_str(j)
                          + ") can be returned: it is a zero entry.");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_Sparse<T, Prop, Storage, Allocator>
  ::AddInteraction(int i, int j, const T& val)
  {
    Val(i, j) += val;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_Sparse<T, Prop, Storage, Allocator>&
  Matrix_Sparse<T, Prop, Storage, Allocator>
  ::operator= (const Matrix_Sparse<T, Prop, Storage, Allocator>& A)
  {
    this->Copy(A);

    return *this;
  }


  /************************
   * CONVENIENT FUNCTIONS *
   ************************/



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::Zero()
  {
    this->allocator_.memoryset(this->data_, char(0),
                               this->nz_ * sizeof(value_type));
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::SetIdentity()
  {
    int m = this->m_;
    int n = this->n_;
    int nz = min(m, n);

    if (nz == 0)
      return;

    Clear();

    Vector<T, VectFull, Allocator> values(nz);
    Vector<int, VectFull, CallocAlloc<int> > ptr(Storage::GetFirst(m, n) + 1);
    Vector<int, VectFull, CallocAlloc<int> > ind(nz);

    values.Fill(T(1));
    ind.Fill();
    int i;
    for (i = 0; i < nz + 1; i++)
      ptr(i) = i;
    for (i = nz + 1; i < ptr.GetLength(); i++)
      ptr(i) = nz;

    SetData(m, n, values, ptr, ind);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::Fill()
  {
    for (int i = 0; i < this->GetDataSize(); i++)
      this->data_[i] = i;
  }



  template <class T, class Prop, class Storage, class Allocator>
  template <class T0>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::Fill(const T0& x)
  {
    for (int i = 0; i < this->GetDataSize(); i++)
      this->data_[i] = x;
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::FillRand()
  {
    srand(time(NULL));
    for (int i = 0; i < this->GetDataSize(); i++)
      this->data_[i] = rand();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::Print() const
  {
    for (int i = 0; i < this->m_; i++)
      {
        for (int j = 0; j < this->n_; j++)
          cout << (*this)(i, j) << "\t";
        cout << endl;
      }
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::
  WriteText(string FileName) const
  {
    ofstream FileStream; FileStream.precision(14);
    FileStream.open(FileName.c_str());

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Matrix_ArraySparse::Write(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->WriteText(FileStream);

    FileStream.close();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_Sparse<T, Prop, Storage, Allocator>::
  WriteText(ostream& FileStream) const
  {

#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Matrix_ArraySparse::Write(ofstream& FileStream)",
                    "Stream is not ready.");
#endif

    // conversion in coordinate format (1-index convention)
    IVect IndRow, IndCol; Vector<T> Value;
    const Matrix<T, Prop, Storage, Allocator>& leaf_class =
      static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this);

    ConvertMatrix_to_Coordinates(leaf_class, IndRow, IndCol,
                                 Value, 1, true);

    for (int i = 0; i < IndRow.GetM(); i++)
      FileStream << IndRow(i) << " " << IndCol(i) << " " << Value(i) << '\n';

  }


  // MATRIX<COLSPARSE> //


  /****************
   * CONSTRUCTORS *
   ****************/


  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSparse, Allocator>::Matrix()  throw():
    Matrix_Sparse<T, Prop, ColSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSparse, Allocator>::Matrix(int i, int j):
    Matrix_Sparse<T, Prop, ColSparse, Allocator>(i, j, 0)
  {
  }



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSparse, Allocator>::Matrix(int i, int j, int nz):
    Matrix_Sparse<T, Prop, ColSparse, Allocator>(i, j, nz)
  {
  }



  template <class T, class Prop, class Allocator>
  template <class Storage0, class Allocator0,
            class Storage1, class Allocator1,
            class Storage2, class Allocator2>
  Matrix<T, Prop, ColSparse, Allocator>::
  Matrix(int i, int j,
         Vector<T, Storage0, Allocator0>& values,
         Vector<int, Storage1, Allocator1>& ptr,
         Vector<int, Storage2, Allocator2>& ind):
    Matrix_Sparse<T, Prop, ColSparse, Allocator>(i, j, values, ptr, ind)
  {
  }



  // MATRIX<ROWSPARSE> //


  /****************
   * CONSTRUCTORS *
   ****************/


  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, RowSparse, Allocator>::Matrix()  throw():
    Matrix_Sparse<T, Prop, RowSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, RowSparse, Allocator>::Matrix(int i, int j):
    Matrix_Sparse<T, Prop, RowSparse, Allocator>(i, j, 0)
  {
  }



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, RowSparse, Allocator>::Matrix(int i, int j, int nz):
    Matrix_Sparse<T, Prop, RowSparse, Allocator>(i, j, nz)
  {
  }



  template <class T, class Prop, class Allocator>
  template <class Storage0, class Allocator0,
            class Storage1, class Allocator1,
            class Storage2, class Allocator2>
  Matrix<T, Prop, RowSparse, Allocator>::
  Matrix(int i, int j,
         Vector<T, Storage0, Allocator0>& values,
         Vector<int, Storage1, Allocator1>& ptr,
         Vector<int, Storage2, Allocator2>& ind):
    Matrix_Sparse<T, Prop, RowSparse, Allocator>(i, j, values, ptr, ind)
  {
  }



  template <class T, class Prop, class Allocator>
  void Matrix<T, Prop, RowSparse, Allocator>::FillRand(int Nelement)
  {
    if (this->m_ == 0 || this->n_ == 0)
      return;

    Vector<int> i(Nelement), j(Nelement);
    Vector<T> value(Nelement);

    set<pair<int, int> > skeleton;
    set<pair<int, int> >::iterator it;

    srand(time(NULL));

    while (static_cast<int>(skeleton.size()) != Nelement)
      skeleton.insert(make_pair(rand() % this->m_, rand() % this->n_));

    int l = 0;
    for (it = skeleton.begin(); it != skeleton.end(); it++)
      {
        i(l) = it->first;
        j(l) = it->second;
        value(l) = double(rand());
        l++;
      }

    ConvertMatrix_from_Coordinates(i, j, value, *this);
  }



  template <class T, class Prop, class Allocator>
  void Matrix<T, Prop, RowSparse, Allocator>
  ::FillRand(int Nelement, const T& x)
  {
    if (this->m_ == 0 || this->n_ == 0)
      return;

    Vector<int> i(Nelement), j(Nelement);
    Vector<T> value(Nelement);
    value.Fill(x);

    srand(time(NULL));
    for (int l = 0; l < Nelement; l++)
      {
        i(l) = rand() % this->m_;
        j(l) = rand() % this->n_;
      }

    ConvertMatrix_from_Coordinates(i, j, value, *this);
  }


} // namespace Seldon.

#define SELDON_FILE_MATRIX_SPARSE_CXX
#endif