matrix_sparse/Matrix_SymSparse.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_SYMSPARSE_CXX

#include "Matrix_SymSparse.hxx"

namespace Seldon
{


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



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



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

    Reallocate(i, j);
  }



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

    Reallocate(i, j, nz);
  }



  template <class T, class Prop, class Storage, class Allocator>
  template <class Storage0, class Allocator0,
            class Storage1, class Allocator1,
            class Storage2, class Allocator2>
  inline Matrix_SymSparse<T, Prop, Storage, Allocator>::
  Matrix_SymSparse(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_SymSparse::Matrix_SymSparse(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 != i)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_SymSparse::Matrix_SymSparse(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(i)
                       + " rows or columns ("
                       + to_str(i) + " by " + to_str(i) + " matrix).");
      }

    if (static_cast<long int>(2 * nz_ - 2) / static_cast<long int>(i + 1)
        >= static_cast<long int>(i))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_SymSparse::Matrix_SymSparse(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(i) + ").");
      }
#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_SymSparse<T, Prop, Storage, Allocator>::
  Matrix_SymSparse(const Matrix_SymSparse<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_SymSparse<T, Prop, Storage, Allocator>::~Matrix_SymSparse()
  {
    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_SymSparse<T, Prop, Storage, Allocator>::Clear()
  {
    this->~Matrix_SymSparse();
  }


  /*********************
   * 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_SymSparse<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_ = i;
    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_SymSparse::Reallocate(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 != i)
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_SymSparse::Reallocate(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")
                       + string(" non-zero entries) but there are ")
                       + to_str(i) + " rows or columns ("
                       + to_str(i) + " by " + to_str(i) + " matrix).");
      }

    if (static_cast<long int>(2 * nz_ - 2) / static_cast<long int>(i + 1)
        >= static_cast<long int>(i))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim(string("Matrix_SymSparse::Reallocate(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(i) + ").");
      }
#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_SymSparse<T, Prop, Storage, Allocator>
  ::SetData(int i, int j, int nz,
            typename Matrix_SymSparse<T, Prop, Storage, Allocator>
            ::pointer values,
            int* ptr,
            int* ind)
  {
    this->Clear();

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

    this->nz_ = nz;

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



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_SymSparse<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>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>::Reallocate(int i, int j)
  {
    // clearing previous entries
    Clear();

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

    // we try to allocate ptr_
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ptr_ = reinterpret_cast<int*>( calloc(i+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_SymSparse::Reallocate(int, int)",
                     string("Unable to allocate ")
                     + to_str(sizeof(int) * (i+1) )
                     + " bytes to store " + to_str(i+1)
                     + " row or column start indices, for a "
                     + to_str(i) + " by " + to_str(i) + " matrix.");
#endif

    // then filing ptr_ with 0
    for (int k = 0; k <= i; k++)
      ptr_[k] = 0;
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::Reallocate(int i, int j, int nz)
  {
    // clearing previous entries
    Clear();

    this->nz_ = nz;
    this->m_ = i;
    this->n_ = i;

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

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ptr_ = reinterpret_cast<int*>( calloc(i + 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(int, int, int)",
                     string("Unable to allocate ")
                     + to_str(sizeof(int) * (i+1) ) + " bytes to store "
                     + to_str(i+1) + " row or column start indices, for a "
                     + to_str(i) + " by " + to_str(i) + " 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(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(i) + " 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(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(i) + " matrix.");
#endif
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>::Resize(int i, int j)
  {
    if (i < this->m_)
      Resize(i, i, ptr_[i]);
    else
      Resize(i, i, nz_);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::Resize(int i, int j, int nz)
  {

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

    if (nz != nz_)
      {
        // trying to resize ind_ and data_
#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif

            ind_ = reinterpret_cast<int*>( realloc(reinterpret_cast<void*>(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_SymSparse::Resize(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

        Vector<T, VectFull, Allocator> val;
        val.SetData(nz_, this->data_);
        val.Resize(nz);

        this->data_ = val.GetData();
        nz_ = nz;
        val.Nullify();
      }


    if (this->m_ != i)
      {
#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif
            // trying to resize ptr_
            ptr_ = reinterpret_cast<int*>( realloc(ptr_, (i+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_SymSparse::Resize(int, int)",
                         string("Unable to allocate ")
                         + to_str(sizeof(int) * (i+1) )
                         + " bytes to store " + to_str(i+1)
                         + " row or column start indices, for a "
                         + to_str(i) + " by " + to_str(i) + " matrix.");
#endif

        // then filing last values of ptr_ with nz_
        for (int k = this->m_; k <= i; k++)
          ptr_[k] = this->nz_;
      }

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


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_SymSparse<T, Prop, Storage, Allocator>::
  Copy(const Matrix_SymSparse<T, Prop, Storage, Allocator>& A)
  {
    this->Clear();
    int nz = A.GetNonZeros();
    int i = A.GetM();
    int j = A.GetN();
    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 (static_cast<long int>(2 * nz_ - 2) / static_cast<long int>(i + 1)
        >= static_cast<long int>(i))
      {
        this->m_ = 0;
        this->n_ = 0;
        nz_ = 0;
        ptr_ = NULL;
        ind_ = NULL;
        this->data_ = NULL;
        throw WrongDim("Matrix_SymSparse::Matrix_SymSparse(int, int, int)",
                       string("There are more values (") + to_str(nz)
                       + string(" values) than elements in the upper")
                       + " part of the matrix ("
                       + to_str(i) + " by " + to_str(i) + ").");
      }
#endif

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        ptr_ = reinterpret_cast<int*>( calloc(i + 1, sizeof(int)) );
        memcpy(this->ptr_, A.ptr_, (i + 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(int, int, int)",
                     string("Unable to allocate ")
                     + to_str(sizeof(int) * (i+1) ) + " bytes to store "
                     + to_str(i+1) + " row or column start indices, for a "
                     + to_str(i) + " by " + to_str(i) + " 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(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(i) + " 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)
      throw NoMemory("Matrix_SymSparse::Matrix_SymSparse(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(i) + " matrix.");
#endif

  }


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



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



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



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



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



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



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


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



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

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_SymSparse::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_SymSparse::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;

    // Only the upper part is stored.
    if (i > j)
      {
        l = i;
        i = j;
        j = l;
      }

    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_SymSparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_SymSparse<T, Prop, Storage, Allocator>::Val(int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_SymSparse::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_SymSparse::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;

    // Only the upper part is stored.
    if (i > j)
      {
        l = i;
        i = j;
        j = l;
      }

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

    if (a == b)
      throw WrongArgument("Matrix_SymSparse::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_SymSparse::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_SymSparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_SymSparse<T, Prop, Storage, Allocator>::Val(int i, int j) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_SymSparse::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_SymSparse::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;

    // Only the upper part is stored.
    if (i > j)
      {
        l = i;
        i = j;
        j = l;
      }

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

    if (a == b)
      throw WrongArgument("Matrix_SymSparse::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_SymSparse::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 typename Matrix_SymSparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_SymSparse<T, Prop, Storage, Allocator>::Get(int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_SymSparse::Get(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_SymSparse::Get(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;
    // Only the upper part is stored.
    if (i > j)
      {
        l = i;
        i = j;
        j = l;
      }

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

    if (a < b)
      {
        l = Storage::GetSecond(i, j);

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

        if ( (k < b) && (ind_[k] == l))
          return this->data_[k];
      }
    else
      k = a;

    // adding a non-zero entry
    Resize(this->m_, this->n_, nz_+1);

    for (int m = Storage::GetFirst(i, j)+1; m <= this->m_; m++)
      ptr_[m]++;

    for (int m = nz_-1; m >= k+1; m--)
      {
        ind_[m] = ind_[m-1];
        this->data_[m] = this->data_[m-1];
      }

    ind_[k] = Storage::GetSecond(i, j);

    // value of new non-zero entry is set to 0
    SetComplexZero(this->data_[k]);

    return this->data_[k];
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline const typename Matrix_SymSparse<T, Prop, Storage, Allocator>::value_type&
  Matrix_SymSparse<T, Prop, Storage, Allocator>::Get(int i, int j) const
  {
    return Val(i, j);
  }



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



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::Set(int i, int j, const T& val)
  {
    Get(i, j) = val;
  }



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

    return *this;
  }


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



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<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_SymSparse<T, Prop, Storage, Allocator>::SetIdentity()
  {
    int m = this->m_;
    int nz = this->m_;

    if (nz == 0)
      return;

    Clear();

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

    values.Fill(T(1));
    ind.Fill();
    ptr.Fill();

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



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<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_SymSparse<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_SymSparse<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_SymSparse<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_SymSparse<T, Prop, Storage, Allocator>
  ::Write(string FileName) const
  {
    ofstream FileStream;
    FileStream.open(FileName.c_str());

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

    this->Write(FileStream);

    FileStream.close();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::Write(ostream& FileStream) const
  {
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Matrix_SymSparse::Write(ofstream& FileStream)",
                    "Stream is not ready.");
#endif

    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->m_)),
                     sizeof(int));
    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->m_)),
                     sizeof(int));
    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->nz_)),
                     sizeof(int));

    FileStream.write(reinterpret_cast<char*>(this->ptr_),
                     sizeof(int)*(this->m_+1));
    FileStream.write(reinterpret_cast<char*>(this->ind_),
                     sizeof(int)*this->nz_);
    FileStream.write(reinterpret_cast<char*>(this->data_),
                     sizeof(T)*this->nz_);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<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_SymSparse::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 Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::WriteText(ostream& FileStream) const
  {

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

    // Conversion to 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';
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>
  ::Read(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Matrix_SymSparse::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 Matrix_SymSparse<T, Prop, Storage, Allocator>::
  Read(istream& FileStream)
  {

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

    int m, n, nz;
    FileStream.read(reinterpret_cast<char*>(&m), sizeof(int));
    FileStream.read(reinterpret_cast<char*>(&n), sizeof(int));
    FileStream.read(reinterpret_cast<char*>(&nz), sizeof(int));

    Reallocate(m, m, nz);

    FileStream.read(reinterpret_cast<char*>(ptr_),
                    sizeof(int)*(m+1));
    FileStream.read(reinterpret_cast<char*>(ind_), sizeof(int)*nz);
    FileStream.read(reinterpret_cast<char*>(this->data_), sizeof(T)*nz);

#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!FileStream.good())
      throw IOError("Matrix_SymSparse::Read(istream& FileStream)",
                    string("Input operation failed.")
                    + string(" The input file may have been removed")
                    + " or may not contain enough data.");
#endif

  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>::
  ReadText(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());

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

    this->ReadText(FileStream);

    FileStream.close();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_SymSparse<T, Prop, Storage, Allocator>::
  ReadText(istream& FileStream)
  {
    Matrix<T, Prop, Storage, Allocator>& leaf_class =
      static_cast<Matrix<T, Prop, Storage, Allocator>& >(*this);

    T zero; int index = 1;
    ReadCoordinateMatrix(leaf_class, FileStream, zero, index);
  }


  // MATRIX<COLSYMSPARSE> //


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


  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSymSparse, Allocator>::Matrix():
    Matrix_SymSparse<T, Prop, ColSymSparse, Allocator>()
  {
  }



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



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSymSparse, Allocator>::Matrix(int i, int j, int nz):
    Matrix_SymSparse<T, Prop, ColSymSparse, 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, ColSymSparse, Allocator>::
  Matrix(int i, int j,
         Vector<T, Storage0, Allocator0>& values,
         Vector<int, Storage1, Allocator1>& ptr,
         Vector<int, Storage2, Allocator2>& ind):
    Matrix_SymSparse<T, Prop, ColSymSparse, Allocator>(i, j, values, ptr, ind)
  {
  }



  // MATRIX<ROWSYMSPARSE> //


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



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, RowSymSparse, Allocator>::Matrix():
    Matrix_SymSparse<T, Prop, RowSymSparse, Allocator>()
  {
  }



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



  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, RowSymSparse, Allocator>::Matrix(int i, int j, int nz):
    Matrix_SymSparse<T, Prop, RowSymSparse, 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, RowSymSparse, Allocator>::
  Matrix(int i, int j,
         Vector<T, Storage0, Allocator0>& values,
         Vector<int, Storage1, Allocator1>& ptr,
         Vector<int, Storage2, Allocator2>& ind):
    Matrix_SymSparse<T, Prop, RowSymSparse, Allocator>(i, j, values, ptr, ind)
  {
  }


} // namespace Seldon.

#define SELDON_FILE_MATRIX_SYMSPARSE_CXX
#endif