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>(i, i)
  {
    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, int nz):
    Matrix_Base<T, Allocator>(i, i)
  {
    this->nz_ = 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

#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>
  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>
  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 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>::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>
  ::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';
  }


  // MATRIX<COLSYMSPARSE> //


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


  template <class T, class Prop, class Allocator>
  Matrix<T, Prop, ColSymSparse, Allocator>::Matrix()  throw():
    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()  throw():
    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