matrix_sparse/Matrix_ArraySparse.cxx

// 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_ARRAY_SPARSE_CXX

#include "Matrix_ArraySparse.hxx"

namespace Seldon
{


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



  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  Matrix_ArraySparse(int i, int j) :
    val_(Storage::GetFirst(i, j))
  {
    this->m_ = i;
    this->n_ = j;
  }


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



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


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



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  Reallocate(int i, int j)
  {
    // Clears previous entries.
    Clear();

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

    int n = Storage::GetFirst(i, j);
    val_.Reallocate(n);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Resize(int i, int j)
  {
    int n = Storage::GetFirst(this->m_, n_);
    int new_n = Storage::GetFirst(i, j);
    if (n != new_n)
      {
        Vector<Vector<T, VectSparse, Allocator>, VectFull,
          NewAlloc<Vector<T, VectSparse, Allocator> > > new_val;

        new_val.Reallocate(new_n);

        for (int k = 0 ; k < min(n, new_n) ; k++)
          Swap(new_val(k), this->val_(k));

        val_.SetData(new_n, new_val.GetData());
        new_val.Nullify();

      }

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


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



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetM() const
  {
    return m_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetN() const
  {
    return n_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::GetM(const SeldonTranspose& status) const
  {
    if (status.NoTrans())
      return m_;
    else
      return n_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::GetN(const SeldonTranspose& status) const
  {
    if (status.NoTrans())
      return n_;
    else
      return m_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetNonZeros()
    const
  {
    int nnz = 0;
    for (int i = 0; i < this->val_.GetM(); i++)
      nnz += this->val_(i).GetM();

    return nnz;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetDataSize()
    const
  {
    return GetNonZeros();
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline int* Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetIndex(int i)
    const
  {
    return val_(i).GetIndex();
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline T*
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::GetData(int i) const
  {
    return val_(i).GetData();
  }



  template <class T, class Prop, class Storage, class Allocat>
  inline Vector<T, VectSparse, Allocat>*
  Matrix_ArraySparse<T, Prop, Storage, Allocat>::GetData() const
  {
    return val_.GetData();
  }


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



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

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

    return this->val_(Storage::GetFirst(i, j))(Storage::GetSecond(i, j));
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline T&
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::operator() (int i, int j)
  {

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

    return this->val_(Storage::GetFirst(i, j))(Storage::GetSecond(i, j));
  }



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

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

    return
      this->val_(Storage::GetFirst(i, j)).Val(Storage::GetSecond(i, j));
  }



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

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

    return
      this->val_(Storage::GetFirst(i, j)).Val(Storage::GetSecond(i, j));
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline const T& Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  Value (int i, int j) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_ArraySparse::value", "Index should be in [0, "
                     + to_str(this->m_-1) + "], but is equal to "
                     + to_str(i) + ".");

    if ((j < 0)||(j >= this->val_(i).GetM()))
      throw WrongCol("Matrix_ArraySparse::value", "Index should be in [0, " +
                     to_str(this->val_(i).GetM()-1) + "], but is equal to "
                     + to_str(j) + ".");
#endif

    return val_(i).Value(j);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline T&
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::Value (int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_ArraySparse::value", "Index should be in [0, "
                     + to_str(this->m_-1) + "], but is equal to "
                     + to_str(i) + ".");

    if ((j < 0)||(j >= this->val_(i).GetM()))
      throw WrongCol("Matrix_ArraySparse::value", "Index should be in [0, " +
                     to_str(this->val_(i).GetM()-1) + "], but is equal to "
                     + to_str(j) + ".");
#endif

    return val_(i).Value(j);
  }



  template <class T, class Prop, class Storage, class Allocator> inline
  int Matrix_ArraySparse<T, Prop, Storage, Allocator>::Index(int i, int j)
    const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_ArraySparse::index", "Index should be in [0, "
                     + to_str(this->m_-1) + "], but is equal to "
                     + to_str(i) + ".");

    if ((j < 0)||(j >= this->val_(i).GetM()))
      throw WrongCol("Matrix_ArraySparse::index", "Index should be in [0, " +
                     to_str(this->val_(i).GetM()-1) + "], but is equal to "
                     + to_str(j) + ".");
#endif

    return val_(i).Index(j);
  }



  template <class T, class Prop, class Storage, class Allocator> inline
  int& Matrix_ArraySparse<T, Prop, Storage, Allocator>::Index(int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_ArraySparse::index", "Index should be in [0, "
                     + to_str(this->m_-1) + "], but is equal to "
                     + to_str(i) + ".");

    if (j < 0 || j >= this->val_(i).GetM())
      throw WrongCol("Matrix_ArraySparse::index", "Index should be in [0, " +
                     to_str(this->val_(i).GetM()-1) + "], but is equal to "
                     + to_str(j) + ".");
#endif

    return val_(i).Index(j);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  SetData(int i, int n, T* val, int* ind)
  {
    val_(i).SetData(n, val, ind);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Nullify(int i)
  {
    val_(i).Nullify();
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  SetData(int m, int n, Vector<T, VectSparse, Allocator>* val)
  {
    m_ = m;
    n_ = n;
    val_.SetData(Storage::GetFirst(m, n), val);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Nullify()
  {
    m_ = 0;
    n_ = 0;
    val_.Nullify();
  }


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



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Print() const
  {
    if (Storage::GetFirst(1, 0) == 1)
      for (int i = 0; i < this->m_; i++)
        {
          for (int j = 0; j < this->val_(i).GetM(); j++)
            cout << (i+1) << " " << this->val_(i).Index(j)+1
                 << " " << this->val_(i).Value(j) << endl;
        }
    else
      for (int i = 0; i < this->n_; i++)
        {
          for (int j = 0; j < this->val_(i).GetM(); j++)
            cout << this->val_(i).Index(j)+1 << " " << i+1
                 << " " << this->val_(i).Value(j) << endl;
        }
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Assemble()
  {
    for (int i = 0; i < val_.GetM(); i++)
      val_(i).Assemble();
  }



  template <class T, class Prop, class Storage, class Allocator>
  template<class T0>
  void Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  RemoveSmallEntry(const T0& epsilon)
  {
    for (int i = 0; i < val_.GetM(); i++)
      val_(i).RemoveSmallEntry(epsilon);
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::SetIdentity()
  {
    this->n_ = this->m_;
    for (int i = 0; i < this->m_; i++)
      {
        val_(i).Reallocate(1);
        val_(i).Index(0) = i;
        val_(i).Value(0) = T(1);
      }
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Zero()
  {
    for (int i = 0; i < val_.GetM(); i++)
      val_(i).Zero();
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::Fill()
  {
    int value = 0;
    for (int i = 0; i < val_.GetM(); i++)
      for (int j = 0; j < val_(i).GetM(); j++)
        val_(i).Value(j) = value++;
  }


  template <class T, class Prop, class Storage, class Allo> template<class T0>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allo>::Fill(const T0& x)
  {
    for (int i = 0; i < val_.GetM(); i++)
      val_(i).Fill(x);
  }


  template <class T, class Prop, class Storage, class Allocator>
  template <class T0>
  inline Matrix_ArraySparse<T, Prop, Storage, Allocator>&
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::operator= (const T0& x)
  {
    this->Fill(x);
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>::FillRand()
  {
    for (int i = 0; i < val_.GetM(); i++)
      val_(i).FillRand();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<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_ArraySparse::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_ArraySparse<T, Prop, Storage, Allocator>::
  Write(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

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

    for (int i = 0; i < this->m_; i++)
      val_(i).Write(FileStream);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<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_ArraySparse<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';

  }



  template <class T, class Prop, class Storage, class Allocator>
  void Matrix_ArraySparse<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_ArraySparse::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_ArraySparse<T, Prop, Storage, Allocator>::
  Read(istream& FileStream)
  {

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

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

    val_.Reallocate(this->m_);
    for (int i = 0; i < this->m_; i++)
      val_(i).Read(FileStream);

#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!FileStream.good())
      throw IOError("Matrix_ArraySparse::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_ArraySparse<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_ArraySparse::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_ArraySparse<T, Prop, Storage, Allocator>::
  ReadText(istream& FileStream)
  {
    // previous elements are removed
    Clear();

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

    Vector<T, VectFull, Allocator> values;
    Vector<int> row_numbers, col_numbers;
    T entry; int row = 0, col;
    int nb_elt = 0;
    while (!FileStream.eof())
      {
        // new entry is read (1-index)
        FileStream >> row >> col >> entry;

        if (FileStream.fail())
          break;
        else
          {
#ifdef SELDON_CHECK_IO
            if (row < 1)
              throw IOError(string("Matrix_ArraySparse::ReadText")+
                            "(istream& FileStream)",
                            string("Error : Row number should be greater ")
                            + "than 0 but is equal to " + to_str(row));

            if (col < 1)
              throw IOError(string("Matrix_ArraySparse::ReadText")+
                            "(istream& FileStream)",
                            string("Error : Column number should be greater")
                            + " than 0 but is equal to " + to_str(col));
#endif

            nb_elt++;

            // inserting new element
            if (nb_elt > values.GetM())
              {
                values.Resize(2*nb_elt);
                row_numbers.Resize(2*nb_elt);
                col_numbers.Resize(2*nb_elt);
              }

            values(nb_elt-1) = entry;
            row_numbers(nb_elt-1) = row;
            col_numbers(nb_elt-1) = col;
          }
      }

    if (nb_elt > 0)
      {
        Matrix<T, Prop, Storage, Allocator>& leaf_class =
          static_cast<Matrix<T, Prop, Storage, Allocator>& >(*this);
        row_numbers.Resize(nb_elt);
        col_numbers.Resize(nb_elt);
        values.Resize(nb_elt);
        ConvertMatrix_from_Coordinates(row_numbers, col_numbers, values,
                                       leaf_class, 1);
      }
  }


  // MATRIX<ARRAY_COLSPARSE> //



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSparse, Allocator>::Matrix()  throw():
    Matrix_ArraySparse<T, Prop, ArrayColSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSparse, Allocator>::Matrix(int i, int j):
    Matrix_ArraySparse<T, Prop, ArrayColSparse, Allocator>(i, j)
  {
  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::ClearColumn(int i)
  {
    this->val_(i).Clear();
  }



  template <class T, class Prop, class Alloc> inline
  void Matrix<T, Prop, ArrayColSparse, Alloc>::ReallocateColumn(int i,int j)
  {
    this->val_(i).Reallocate(j);
  }



  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayColSparse, Allocator>::ResizeColumn(int i,int j)
  {
    this->val_(i).Resize(j);
  }



  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayColSparse, Allocator>::SwapColumn(int i,int j)
  {
    Swap(this->val_(i), this->val_(j));
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  ReplaceIndexColumn(int i, IVect& new_index)
  {
    for (int j = 0; j < this->val_(i).GetM(); j++)
      this->val_(i).Index(j) = new_index(j);
  }



  template <class T, class Prop, class Allocator>
  inline int Matrix<T, Prop, ArrayColSparse, Allocator>::
  GetColumnSize(int i) const
  {
    return this->val_(i).GetSize();
  }


  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayColSparse, Allocator>::PrintColumn(int i) const
  {
    this->val_(i).Print();
  }



  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayColSparse, Allocator>::AssembleColumn(int i)
  {
    this->val_(i).Assemble();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  AddInteraction(int i, int j, const T& val)
  {
    this->val_(j).AddInteraction(i, val);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  AddInteractionRow(int i, int nb, int* col_, T* value_)
  {
    IVect col;
    col.SetData(nb, col_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionRow(i, nb, col, val);
    col.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  AddInteractionColumn(int i, int nb, int* row_, T* value_)
  {
    IVect row;
    row.SetData(nb, row_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionColumn(i, nb, row, val);
    row.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  AddInteractionRow(int i, int nb, const IVect& col,
                    const Vector<T, VectFull, Alloc1>& val)
  {
    for (int j = 0; j < nb; j++)
      this->val_(col(j)).AddInteraction(i, val(j));
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayColSparse, Allocator>::
  AddInteractionColumn(int i, int nb, const IVect& row,
                       const Vector<T, VectFull, Alloc1>& val)
  {
    this->val_(i).AddInteractionRow(nb, row, val);
  }


  // MATRIX<ARRAY_ROWSPARSE> //



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSparse, Allocator>::Matrix()  throw():
    Matrix_ArraySparse<T, Prop, ArrayRowSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSparse, Allocator>::Matrix(int i, int j):
    Matrix_ArraySparse<T, Prop, ArrayRowSparse, Allocator>(i, j)
  {
  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::ClearRow(int i)
  {
    this->val_(i).Clear();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  ReallocateRow(int i, int j)
  {
    this->val_(i).Reallocate(j);
  }



  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayRowSparse, Allocator>::ResizeRow(int i, int j)
  {
    this->val_(i).Resize(j);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::SwapRow(int i,int j)
  {
    Swap(this->val_(i), this->val_(j));
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  ReplaceIndexRow(int i, IVect& new_index)
  {
    for (int j = 0; j < this->val_(i).GetM(); j++)
      this->val_(i).Index(j) = new_index(j);
  }



  template <class T, class Prop, class Allocator> inline
  int Matrix<T, Prop, ArrayRowSparse, Allocator>::GetRowSize(int i) const
  {
    return this->val_(i).GetSize();
  }


  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayRowSparse, Allocator>::PrintRow(int i) const
  {
    this->val_(i).Print();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::AssembleRow(int i)
  {
    this->val_(i).Assemble();
  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  AddInteraction(int i, int j, const T& val)
  {
    this->val_(i).AddInteraction(j, val);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  AddInteractionRow(int i, int nb, int* col_, T* value_)
  {
    IVect col;
    col.SetData(nb, col_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionRow(i, nb, col, val);
    col.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  AddInteractionColumn(int i, int nb, int* row_, T* value_)
  {
    IVect row;
    row.SetData(nb, row_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionColumn(i, nb, row, val);
    row.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  AddInteractionRow(int i, int nb, const IVect& col,
                    const Vector<T, VectFull, Alloc1>& val)
  {
    this->val_(i).AddInteractionRow(nb, col, val);
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayRowSparse, Allocator>::
  AddInteractionColumn(int i, int nb, const IVect& row,
                       const Vector<T, VectFull, Alloc1>& val)
  {
    for (int j = 0; j < nb; j++)
      this->val_(row(j)).AddInteraction(i, val(j));
  }


  // MATRIX<ARRAY_COLSYMSPARSE> //



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSymSparse, Allocator>::Matrix()  throw():
    Matrix_ArraySparse<T, Prop, ArrayColSymSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSymSparse, Allocator>::Matrix(int i, int j):
    Matrix_ArraySparse<T, Prop, ArrayColSymSparse, Allocator>(i, j)
  {
  }


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



  template <class T, class Prop, class Allocator>
  inline T
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::operator() (int i, int j)
    const
  {

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

    if (i <= j)
      return this->val_(j)(i);

    return this->val_(i)(j);
  }



  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::operator() (int i, int j)
  {

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

    if (i < j)
      return this->val_(j)(i);

    return this->val_(i)(j);
  }


  template <class T, class Prop, class Allocator> inline
  void Matrix<T, Prop, ArrayColSymSparse, Allocator>::ClearColumn(int i)
  {
    this->val_(i).Clear();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  ReallocateColumn(int i, int j)
  {
    this->val_(i).Reallocate(j);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  ResizeColumn(int i, int j)
  {
    this->val_(i).Resize(j);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  SwapColumn(int i, int j)
  {
    Swap(this->val_(i), this->val_(j));
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  ReplaceIndexColumn(int i, IVect& new_index)
  {
    for (int j = 0; j < this->val_(i).GetM(); j++)
      this->val_(i).Index(j) = new_index(j);
  }



  template <class T, class Prop, class Allocator>
  inline int Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  GetColumnSize(int i) const
  {
    return this->val_(i).GetSize();
  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  PrintColumn(int i) const
  {
    this->val_(i).Print();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AssembleColumn(int i)
  {
    this->val_(i).Assemble();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AddInteractionRow(int i, int nb, int* col_, T* value_)
  {
    IVect col;
    col.SetData(nb, col_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionRow(i, nb, col, val);
    col.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AddInteractionColumn(int i, int nb, int* row_, T* value_)
  {
    IVect row;
    row.SetData(nb, row_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionColumn(i, nb, row, val);
    row.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AddInteraction(int i, int j, const T& val)
  {
    if (i <= j)
      this->val_(j).AddInteraction(i, val);
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AddInteractionRow(int i, int nb, const IVect& col,
                    const Vector<T, VectFull, Alloc1>& val)
  {
    AddInteractionColumn(i, nb, col, val);
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayColSymSparse, Allocator>::
  AddInteractionColumn(int i, int nb, const IVect& row,
                       const Vector<T, VectFull, Alloc1>& val)
  {
    IVect new_row(nb);
    Vector<T, VectFull, Alloc1> new_val(nb);
    nb = 0;
    for (int j = 0; j < new_row.GetM(); j++)
      if (row(j) <= i)
        {
          new_row(nb) = row(j);
          new_val(nb) = val(j); nb++;
        }

    this->val_(i).AddInteractionRow(nb, new_row, new_val);
  }


  // MATRIX<ARRAY_ROWSYMSPARSE> //



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Matrix()  throw():
    Matrix_ArraySparse<T, Prop, ArrayRowSymSparse, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Matrix(int i, int j):
    Matrix_ArraySparse<T, Prop, ArrayRowSymSparse, Allocator>(i, j)
  {
  }


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



  template <class T, class Prop, class Allocator>
  inline T
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::operator() (int i, int j)
    const
  {

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

    if (i <= j)
      return this->val_(i)(j);

    return this->val_(j)(i);
  }



  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::operator() (int i, int j)
  {

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

    if (i <= j)
      return this->val_(i)(j);

    return this->val_(j)(i);

  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::ClearRow(int i)
  {
    this->val_(i).Clear();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  ReallocateRow(int i,int j)
  {
    this->val_(i).Reallocate(j);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  ResizeRow(int i,int j)
  {
    this->val_(i).Resize(j);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  SwapRow(int i,int j)
  {
    Swap(this->val_(i), this->val_(j));
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  ReplaceIndexRow(int i,IVect& new_index)
  {
    for (int j = 0; j < this->val_(i).GetM(); j++)
      this->val_(i).Index(j) = new_index(j);
  }



  template <class T, class Prop, class Allocator>
  inline int Matrix<T, Prop, ArrayRowSymSparse, Allocator>::GetRowSize(int i)
    const
  {
    return this->val_(i).GetSize();
  }


  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::PrintRow(int i)
    const
  {
    this->val_(i).Print();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>
  ::AssembleRow(int i)
  {
    this->val_(i).Assemble();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  AddInteraction(int i, int j, const T& val)
  {
    if (i <= j)
      this->val_(i).AddInteraction(j, val);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  AddInteractionRow(int i, int nb, int* col_, T* value_)
  {
    IVect col;
    col.SetData(nb, col_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionRow(i, nb, col, val);
    col.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  AddInteractionColumn(int i, int nb, int* row_, T* value_)
  {
    IVect row;
    row.SetData(nb, row_);
    Vector<T> val;
    val.SetData(nb, value_);
    AddInteractionColumn(i, nb, row, val);
    row.Nullify();
    val.Nullify();
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  AddInteractionRow(int i, int nb, const IVect& col,
                    const Vector<T, VectFull, Alloc1>& val)
  {
    IVect new_col(nb);
    Vector<T, VectFull, Alloc1> new_val(nb);
    nb = 0;
    for (int j = 0; j < new_col.GetM(); j++)
      if (i <= col(j))
        {
          new_col(nb) = col(j);
          new_val(nb) = val(j); nb++;
        }

    this->val_(i).AddInteractionRow(nb, new_col, new_val);
  }



  template <class T, class Prop, class Allocator> template <class Alloc1>
  inline void Matrix<T, Prop, ArrayRowSymSparse, Allocator>::
  AddInteractionColumn(int i, int nb, const IVect& row,
                       const Vector<T,VectFull,Alloc1>& val)
  {
    // Symmetric matrix, row = column.
    AddInteractionRow(i, nb, row, val);
  }


  template <class T, class Prop, class Allocator>
  ostream& operator <<(ostream& out,
                       const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
  {
    A.WriteText(out);

    return out;
  }


  template <class T, class Prop, class Allocator>
  ostream& operator <<(ostream& out,
                       const Matrix<T, Prop, ArrayColSparse, Allocator>& A)
  {
    A.WriteText(out);

    return out;
  }


  template <class T, class Prop, class Allocator>
  ostream& operator <<(ostream& out,
                       const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
  {
    A.WriteText(out);

    return out;
  }


  template <class T, class Prop, class Allocator>
  ostream& operator <<(ostream& out,
                       const Matrix<T, Prop, ArrayColSymSparse, Allocator>& A)
  {
    A.WriteText(out);

    return out;
  }


} // namespace Seldon

#define SELDON_FILE_MATRIX_ARRAY_SPARSE_CXX
#endif