Matrix_ArraySparseInline.cxx

// Copyright (C) 2003-2011 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_INLINE_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()
    : VirtualMatrix<T>(), val_()
  {
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_ArraySparse<T, Prop, Storage, Allocator>::
  Matrix_ArraySparse(int i, int j) :
    VirtualMatrix<T>(i, j), val_(Storage::GetFirst(i, 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->val_.Clear();
    this->m_ = 0;
    this->n_ = 0;
  }
 
 
  /*******************
   * BASIC FUNCTIONS *
   *******************/
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  inline long 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 *
   **********************************/
  
  
#ifdef SELDON_WITH_MODIFIABLE_PARENTHESIS_OPERATOR
  template <class T, class Prop, class Storage, class Allocator>
  inline T&
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::operator() (int i, int j)  
  {
    return Get(i, j);
  }
#endif
  
  
 
  template <class T, class Prop, class Storage, class Allocator>
  inline const T
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::operator() (int i, int j)
    const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix_ArraySparse");
#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>::Get(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix_ArraySparse");
#endif
 
    return this->val_(Storage::GetFirst(i, j)).Get(Storage::GetSecond(i, j));
  }
 
  
 
  template <class T, class Prop, class Storage, class Allocator>
  inline const T&
  Matrix_ArraySparse<T, Prop, Storage, Allocator>::Get(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix_ArraySparse");
#endif
 
    return this->val_(Storage::GetFirst(i, j)).Get(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
    CheckBounds(i, j, this->m_, this->n_, "Matrix_ArraySparse");
#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
    CheckBounds(i, j, this->m_, this->n_, "Matrix_ArraySparse");
#endif
 
    return
      this->val_(Storage::GetFirst(i, j)).Val(Storage::GetSecond(i, j));
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::Set(int i, int j, const T& x)
  {
    this->Get(i, j) = x;
  }
 
 
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::SetEntry(int i, int j, const T& x)
  {
    this->Get(i, j) = x;
  }
 
 
 
  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
    CheckBounds(i, j, Storage::GetFirst(this->m_, this->n_),
                this->val_(i).GetM(), "Matrix_ArraySparse");
#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
    CheckBounds(i, j, Storage::GetFirst(this->m_, this->n_),
                this->val_(i).GetM(), "Matrix_ArraySparse");
#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
    CheckBounds(i, j, Storage::GetFirst(this->m_, this->n_),
                this->val_(i).GetM(), "Matrix_ArraySparse");
#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
    CheckBounds(i, j, Storage::GetFirst(this->m_, this->n_),
                this->val_(i).GetM(), "Matrix_ArraySparse");
#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)
  {
    this->m_ = m;
    this->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()
  {
    this->m_ = 0;
    this->n_ = 0;
    val_.Nullify();
  }
 
 
#ifdef SELDON_WITH_VIRTUAL
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::ApplySor(const SeldonTranspose& trans, Vector<Treal>& x, const Vector<Treal>& r,
             const typename ClassComplexType<T>::Treal& omega,
             int nb_iter, int stage_ssor) const
  {
    SOR(trans, static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
        x, r, omega, nb_iter, stage_ssor);
  }
  
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::ApplySor(const SeldonTranspose& trans, Vector<Tcplx>& x, const Vector<Tcplx>& r,
             const typename ClassComplexType<T>::Treal& omega,
             int nb_iter, int stage_ssor) const
  {
    SOR(trans,
        static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
        x, r, omega, nb_iter, stage_ssor);
  }
  
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltAddVector(const Treal& alpha, const Vector<Treal>& x,
                 const Treal& beta, Vector<Treal>& y) const
  {
    MltAdd(alpha,
           static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
           x, beta, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltAddVector(const Tcplx& alpha, const Vector<Tcplx>& x,
                 const Tcplx& beta, Vector<Tcplx>& y) const
  {
    MltAdd(alpha,
           static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
           x, beta, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltAddVector(const Treal& alpha, const SeldonTranspose& trans,
                 const Vector<Treal>& x,
                 const Treal& beta, Vector<Treal>& y) const
  {
    MltAdd(alpha, trans,
           static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
           x, beta, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltAddVector(const Tcplx& alpha, const SeldonTranspose& trans,
                 const Vector<Tcplx>& x,
                 const Tcplx& beta, Vector<Tcplx>& y) const
  {
    MltAdd(alpha, trans,
           static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this),
           x, beta, y);
  }
  
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltVector(const Vector<Treal>& x, Vector<Treal>& y) const
  {
    Mlt(static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this), x, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>
  ::MltVector(const Vector<Tcplx>& x, Vector<Tcplx>& y) const
  {
    Mlt(static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this), x, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>  
  ::MltVector(const SeldonTranspose& trans,
              const Vector<Treal>& x, Vector<Treal>& y) const
  {
    Mlt(trans,
        static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this), x, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_ArraySparse<T, Prop, Storage, Allocator>  
  ::MltVector(const SeldonTranspose& trans,
              const Vector<Tcplx>& x, Vector<Tcplx>& y) const
  {
    Mlt(trans,
        static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this), x, y);
  }
 
  template <class T, class Prop, class Storage, class Allocator>
  inline bool Matrix_ArraySparse<T, Prop, Storage, Allocator>  
  ::IsSymmetric() const
  {
    return IsSymmetricMatrix(static_cast<const Matrix<T, Prop, Storage,
                             Allocator>& >(*this));
  }
#endif
 
 
  // MATRIX<ARRAY_COLSPARSE> //
 
 
 
  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSparse, Allocator>::Matrix():
    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);
  }
 
 
  // MATRIX<ARRAY_ROWSPARSE> //
 
 
 
  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSparse, Allocator>::Matrix():
    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);
  }
 
  
  // MATRIX<ARRAY_COLSYMSPARSE> //
 
 
 
  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayColSymSparse, Allocator>::Matrix():
    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 *
   **********************************/
 
 
#ifdef SELDON_WITH_MODIFIABLE_PARENTHESIS_OPERATOR
  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::operator() (int i, int j)  
  {
    return Get(i, j);
  }
#endif
  
  
 
  template <class T, class Prop, class Allocator>
  inline const T
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::operator() (int i, int j)
    const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#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>::Get(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
 
    if (i < j)
      return this->val_(j).Get(i);
    
    return this->val_(i).Get(j);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline const T&
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::Get(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
 
    if (i < j)
      return this->val_(j).Get(i);
    
    return this->val_(i).Get(j);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::Val(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
    
    if (i > j)
      return this->val_(i).Val(j);
    
    return this->val_(j).Val(i);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline const T&
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::Val(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
    
    if (i > j)
      return this->val_(i).Val(j);
    
    return this->val_(j).Val(i);
  }
  
  
 
  template <class T, class Prop, class Allocator>
  inline void
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::Set(int i, int j, const T& x)
  {
    if (i < j)
      this->val_(j).Get(i) = x;
    else
      this->val_(i).Get(j) = x;
  }
 
 
 
  template <class T, class Prop, class Allocator>
  inline void
  Matrix<T, Prop, ArrayColSymSparse, Allocator>::SetEntry(int i, int j, const T& x)
  {
    if (i < j)
      this->val_(j).Get(i) = x;
    else
      this->val_(i).Get(j) = x;
  }
  
  
  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>::
  AddInteraction(int i, int j, const T& val)
  {
    if (i <= j)
      this->val_(j).AddInteraction(i, val);
  }
 
  
  // MATRIX<ARRAY_ROWSYMSPARSE> //
 
 
 
  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Matrix():
    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 *
   **********************************/
 
 
#ifdef SELDON_WITH_MODIFIABLE_PARENTHESIS_OPERATOR
  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::operator() (int i, int j)  
  {
    return Get(i, j);
  }
#endif
 
  
 
  template <class T, class Prop, class Allocator>
  inline const T
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::operator() (int i, int j)
    const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#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>::Get(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
 
    if (i < j)
      return this->val_(i).Get(j);
    
    return this->val_(j).Get(i);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline const T&
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Get(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
 
    if (i < j)
      return this->val_(i).Get(j);
    
    return this->val_(j).Get(i);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline T&
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Val(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
    
    if (i > j)
      return this->val_(j).Val(i);
    
    return this->val_(i).Val(j);
  }
 
  
 
  template <class T, class Prop, class Allocator>
  inline const T&
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Val(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, this->m_, this->n_, "Matrix");
#endif
    
    if (i > j)
      return this->val_(j).Val(i);
    
    return this->val_(i).Val(j);
  }
  
  
 
  template <class T, class Prop, class Allocator>
  inline void
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::Set(int i, int j, const T& x)
  {
    if (i < j)
      this->val_(i).Get(j) = x;
    else
      this->val_(j).Get(i) = x;
  }
 
 
 
  template <class T, class Prop, class Allocator>
  inline void
  Matrix<T, Prop, ArrayRowSymSparse, Allocator>::SetEntry(int i, int j, const T& x)
  {
    if (i < j)
      this->val_(i).Get(j) = x;
    else
      this->val_(j).Get(i) = x;
  }
  
  
  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 ostream& operator <<(ostream& out,
                              const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
  {
    A.WriteText(out);
 
    return out;
  }
 
 
  template <class T, class Prop, class Allocator>
  inline ostream& operator <<(ostream& out,
                              const Matrix<T, Prop, ArrayColSparse, Allocator>& A)
  {
    A.WriteText(out);
 
    return out;
  }
 
 
  template <class T, class Prop, class Allocator>
  inline ostream& operator <<(ostream& out,
                              const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
  {
    A.WriteText(out);
 
    return out;
  }
 
 
  template <class T, class Prop, class Allocator>
  inline ostream& operator <<(ostream& out,
                              const Matrix<T, Prop, ArrayColSymSparse, Allocator>& A)
  {
    A.WriteText(out);
 
    return out;
  }
  
} // namespace Seldon
 
#define SELDON_FILE_MATRIX_ARRAY_SPARSE_INLINE_CXX
#endif