SparseVector.cxx

// Copyright (C) 2003-2011 Marc Duruflé
// Copyright (C) 2001-2011 Vivien Mallet
//
// 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_SPARSE_VECTOR_CXX
 
#include "SparseVector.hxx"
 
namespace Seldon
{
 
  /*********************
   * MEMORY MANAGEMENT *
   *********************/
 
 
 
  template <class T, class Allocator>
  Vector<T, VectSparse, Allocator>::Vector(size_t i):
    Vector<T, VectFull, Allocator>(i)
  {
 
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif
 
        this->index_ = AllocatorInt::allocate(i, this);
 
#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->index_ = NULL;
        this->data_ = NULL;
      }
 
    if (this->index_ == NULL)
      {
        this->m_ = 0;
        this->data_ = NULL;
      }
 
    if (this->data_ == NULL && i != 0)
      throw NoMemory("Vector<VectSparse>::Vector(int)",
                     string("Unable to allocate memory for a vector of size ")
                     + to_str(i * sizeof(T)) + " bytes ("
                     + to_str(i) + " elements).");
#endif
 
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Clear()
  {
    // 'data_' is released.
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif
        if (this->data_ != NULL)
          {
            Allocator::deallocate(this->data_, this->m_);
            this->data_ = NULL;
          }
 
        if (index_ != NULL)
          {
            AllocatorInt::deallocate(index_, this->m_);
            index_ = NULL;
          }
 
        this->m_ = 0;
 
#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->data_ = NULL;
        index_ = NULL;
        this->m_ = 0;
        return;
      }
#endif
 
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::ReallocateVector(size_t i)
  {
    // function implemented in the aim that explicit specialization
    // of Reallocate can call ReallocateVector
    if (i != this->m_)
      {
 
        this->m_ = i;
 
#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif
 
            this->data_ =
              reinterpret_cast<pointer>(Allocator::
                                        reallocate(this->data_, i, this));
 
            index_
              = reinterpret_cast<int*>(AllocatorInt::
                                       reallocate(index_, i, this));
 
#ifdef SELDON_CHECK_MEMORY
          }
        catch (...)
          {
            this->m_ = 0;
            this->data_ = NULL;
            this->index_ = NULL;
            return;
          }
        if (this->data_ == NULL)
          {
            this->m_ = 0;
            this->index_ = NULL;
            return;
          }
#endif
 
      }
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Resize(size_t n)
  {
    ResizeVector(n);
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::ResizeVector(size_t n)
  {
    // function implemented in the aim that explicit specialization
    // of Resize can call ResizeVector
    if (n == this->m_)
      return;
 
    Vector<T, VectFull, Allocator> new_value(n);
    Vector<int> new_index(n);
    size_t Nmin = min(this->m_, n);
    for (size_t i = 0; i < Nmin; i++)
      {
        new_value(i) = this->data_[i];
        new_index(i) = index_[i];
      }
 
    SetData(new_value, new_index);
  }
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>
  ::SetData(size_t i, T* data, int* index)
  {
    this->Clear();
 
    this->m_ = i;
 
    this->data_ = data;
    this->index_ = index;
  }
 
 
  template <class T, class Allocator>
  template<class Allocator2>
  void Vector<T, VectSparse, Allocator>
  ::SetData(Vector<T, VectFull, Allocator2>& data, Vector<int>& index)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    if (data.GetM() != index.GetM())
      throw WrongDim("Vector<VectSparse>::SetData ",
                     string("The data vector and the index vector should")
                     + " have the same size.\n  Size of the data vector: "
                     + to_str(data.GetM()) + "\n  Size of index vector: "
                     + to_str(index.GetM()));
#endif
 
    SetData(data.GetM(), data.GetData(), index.GetData());
    data.Nullify();
    index.Nullify();
  }
 
 
  template <class T, class Allocator>
  template<class Allocator2>
  void Vector<T, VectSparse, Allocator>
  ::SetData(const Vector<T, VectSparse, Allocator2>& V)
  {
    SetData(V.GetM(), V.GetData(), V.GetIndex());
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Nullify()
  {
    this->m_ = 0;
    this->data_ = NULL;
    this->index_ = NULL;
  }
 
 
  /**********************************
   * ELEMENT ACCESS AND AFFECTATION *
   **********************************/
 
 
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::value_type
  Vector<T, VectSparse, Allocator>::operator() (int i) const
  {
    size_t k = 0;
    T zero;
    SetComplexZero(zero);
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      // The entry does not exist, a zero is returned.
      return zero;
 
    return this->data_[k];
  }
 
 
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::value_type
  Vector<T, VectSparse, Allocator>::operator() (int i)
  {
    size_t k = 0;
    T zero;
    SetComplexZero(zero);
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      // The entry does not exist, a zero is returned.
      return zero;
    
    return this->data_[k];
  }
  
  
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::reference
  Vector<T, VectSparse, Allocator>::Get(int i)
  {
    size_t k = 0;
    T zero;
    SetComplexZero(zero);
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      // The entry does not exist yet, so a zero entry is introduced.
      AddInteraction(i, zero);
    
    return this->data_[k];
  }
 
 
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::const_reference
  Vector<T, VectSparse, Allocator>::Get(int i) const
  {
    size_t k = 0;
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      throw WrongArgument("Vector<VectSparse>::Val(int)",
                          "No reference to element " + to_str(i)
                          + " can be returned: it is a zero entry.");
    
    return this->data_[k];
  }
 
 
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::reference
  Vector<T, VectSparse, Allocator>::Val(int i)
  {
    size_t k = 0;
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      throw WrongArgument("Vector<VectSparse>::Val(int)",
                          "the entry " + to_str(i) +
                          " does not belong to the sparsity pattern.");
 
    
    return this->data_[k];
  }
 
 
 
  template <class T, class Allocator>
  typename Vector<T, VectSparse, Allocator>::const_reference
  Vector<T, VectSparse, Allocator>::Val(int i) const
  {
    size_t k = 0;
    // Searching for the entry.
    while (k < this->m_ && index_[k] < i)
      k++;
 
    if (k >= this->m_ || index_[k] != i)
      throw WrongArgument("Vector<VectSparse>::Val(int)",
                          "the entry " + to_str(i) +
                          " does not belong to the sparsity pattern.");
    
    return this->data_[k];
  }
 
  
 
  template <class T, class Allocator> template <class T2, class Alloc2>
  void Vector<T, VectSparse, Allocator>
  ::Copy(const Vector<T2, VectSparse, Alloc2>& X)
  {
    this->Reallocate(X.GetLength());
    AllocatorInt::memorycpy(this->index_, X.GetIndex(), this->m_);
    for (size_t i = 0; i < X.GetLength(); i++)
      this->Value(i) = X.Value(i);
  }
 
 
  /************************
   * CONVENIENT FUNCTIONS *
   ************************/
 
 
 
  template <class T, class Allocator>
  template <class T0>
  Vector<T, VectSparse, Allocator>&
  Vector<T, VectSparse, Allocator>::operator= (const T0& x)
  {
    this->Fill(x);
 
    return *this;
  }
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Print() const
  {
    for (int i = 0; i < this->GetM(); i++)
      cout << (Index(i) + 1) << ' ' << Value(i) << '\n';
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Assemble()
  {
    int new_size = this->m_;
    Vector<T, VectFull, Allocator> values(new_size);
    Vector<int> index(new_size);
    for (int i = 0; i < new_size; i++)
      {
        values(i) = this->data_[i];
        index(i) = index_[i];
      }
 
    Seldon::Assemble(new_size, index, values);
    index.Resize(new_size);
    values.Resize(new_size);
    SetData(values, index);
  }
 
 
 
  template <class T, class Allocator> template<class T0>
  void Vector<T, VectSparse, Allocator>::RemoveSmallEntry(const T0& epsilon)
  {
    size_t new_size = this->m_;
    Vector<T, VectFull, Allocator> values(new_size);
    Vector<int> index(new_size);
    new_size = 0;
    for (size_t i = 0; i < this->m_; i++)
      if (abs(this->data_[i]) > epsilon)
        {
          values(new_size) = this->data_[i];
          index(new_size) = index_[i];
          new_size++;
        }
 
    index.Resize(new_size);
    values.Resize(new_size);
    SetData(values, index);
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::AddInteraction(int i, const T& val)
  {
    // Searching for the position where the entry may be.
    size_t pos = 0;
    while (pos < this->m_ && index_[pos] < i)
      pos++;
 
    // If the entry already exists, adds 'val'.
    if (pos < this->m_ && index_[pos] == i)
      {
        this->data_[pos] += val;
        return;
      }
 
    size_t k;
 
    // If the entry does not exist, the vector is reallocated.
    Resize(this->m_ + 1);
    
    if (this->m_ >= 2)
      for (k = this->m_-1; k > pos; k--)
        {
          this->data_[k] = this->data_[k-1];
          this->index_[k] = this->index_[k-1];
        }
    
    // The new entry.
    this->index_[pos] = i;
    this->data_[pos] = val;
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::
  AddInteractionRow(size_t n, int* index, T* value, bool already_sorted)
  {
    Vector<int> ind;
    Vector<T, VectFull, Allocator> val;
    ind.SetData(n, index);
    val.SetData(n, value);
    AddInteractionRow(n, ind, val, already_sorted);
    ind.Nullify();
    val.Nullify();
  }
 
 
 
  template <class T, class Allocator>
  template<class Allocator0>
  void Vector<T, VectSparse, Allocator>::
  AddInteractionRow(size_t n, const Vector<int>& index2,
                    const Vector<T, VectFull, Allocator0>& value2,
                    bool already_sorted)
  {
    Vector<int> index;
    Vector<T, VectFull, Allocator0> value;
 
    if (!already_sorted)
      {
        // checking if numbers are sorted
        already_sorted = true;
        if (n >= 2)
          for (size_t i = 0; i < n-1; i++)
            if (index2(i+1) <= index2(i))
              already_sorted = false;
      }
    
    if (!already_sorted)
      {
        index.Reallocate(n);
        value.Reallocate(n);
        for (size_t i = 0; i < n; i++)
          {
            index(i) = index2(i);
            value(i) = value2(i);
          }
        
        // Sorts the values to be added according to their indices.
        int ni = n; Seldon::Assemble(ni, index, value); n = ni;
      }
    else
      {
        index.SetData(n, index2.GetData());
        value.SetData(n, value2.GetData());
      }
    
    //  Values that already have an entry 
    
    // Number of values to be added without entry.
    size_t Nnew = 0;
    Vector<bool> new_index(n);
    new_index.Fill(true);
    size_t k = 0;
    for (size_t j = 0; j < n; j++)
      {
        while (k < this->m_ && index_[k] < index(j))
          k++;
 
        if (k < this->m_ && index(j) == index_[k])
          {
            new_index(j) = false;
            this->data_[k] += value(j);
          }
        else
          Nnew++;
      }
 
    if (Nnew > 0)
      {
        // Some values to be added have no entry yet.
        Vector<T> new_val(this->m_ + Nnew);
        Vector<int> new_ind(this->m_ + Nnew);
        int nb = 0;
        k = 0;
        for (size_t j = 0; j < n; j++)
          if (new_index(j))
            {
              while (k < this->m_ && index_[k] < index(j))
                {
                  new_ind(nb) = index_[k];
                  new_val(nb) = this->data_[k];
                  k++;
                  nb++;
                }
 
              // The new entry.
              new_ind(nb) = index(j);
              new_val(nb) = value(j);
              nb++;
            }
 
        // Last entries.
        while (k < this->m_)
          {
            new_ind(nb) = index_[k];
            new_val(nb) = this->data_[k];
            k++;
            nb++;
          }
 
        SetData(new_val, new_ind);
      }
 
    if (already_sorted)
      {
        index.Nullify();
        value.Nullify();
      }
  }
 
  
 
  template <class T, class Allocator>
  typename ClassComplexType<T>::Treal
  Vector<T, VectSparse, Allocator>::GetNormInf() const
  {
    typename ClassComplexType<T>::Treal res(0);
    for (size_t i = 0; i < this->m_; i++)
      res = max(res, abs(this->data_[i]));
    
    return res;
  }
 
 
 
  template <class T, class Allocator>
  int Vector<T, VectSparse, Allocator>::GetNormInfIndex() const
  {
 
#ifdef SELDON_CHECK_DIMENSIONS
    if (this->GetLength() == 0)
      throw WrongDim("Vector<VectSparse>::GetNormInfIndex()",
                     "Vector is null.");
#endif
 
    typename ClassComplexType<T>::Treal res(0), temp;
    size_t j = 0;
    for (size_t i = 0; i < this->GetLength(); i++)
      {
        temp = res;
        res = max(res, abs(this->data_[i]));
        if (temp != res) j = i;
      }
 
    return this->index_[j];
  }
 
  
  /**************************
   * OUTPUT/INPUT FUNCTIONS *
   **************************/
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Write(string FileName) const
  {
    ofstream FileStream;
    FileStream.open(FileName.c_str(), ofstream::binary);
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<VectSparse>::Write(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->Write(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Write(ostream& stream) const
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::Write(ostream& stream)",
                    "Stream is not ready.");
#endif
 
    int n = this->m_;
    stream.write(reinterpret_cast<char*>(&n),
                 sizeof(int));
 
    stream.write(reinterpret_cast<char*>(this->index_),
                 this->m_ * sizeof(int));
 
    stream.write(reinterpret_cast<char*>(this->data_),
                 this->m_ * sizeof(value_type));
 
#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::Write(ostream& stream)",
                    "Output operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::WriteText(string FileName) const
  {
    ofstream FileStream;
    FileStream.precision(cout.precision());
    FileStream.flags(cout.flags());
    FileStream.open(FileName.c_str());
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<VectSparse>::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->WriteText(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::WriteText(ostream& stream) const
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::WriteText(ostream& stream)",
                    "Stream is not ready.");
#endif
 
    // First entries.
    for (size_t i = 0; i < this->m_; i++)
      stream << (Index(i) + 1) << " " << Value(i) << '\n';
    
#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::WriteText(ostream& stream)",
                    "Output operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Read(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str(), ifstream::binary);
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<VectSparse>::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->Read(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::Read(istream& stream)
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::Read(istream& stream)",
                    "Stream is not ready.");
#endif
 
    int m;
    stream.read(reinterpret_cast<char*>(&m), sizeof(int));
    this->Reallocate(m);
 
    stream.read(reinterpret_cast<char*>(this->index_), m * sizeof(int));
 
    stream.read(reinterpret_cast<char*>(this->data_),
                m * sizeof(value_type));
 
#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::Read(istream& stream)",
                    "Input operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, 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("Vector<VectSparse>::ReadText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    this->ReadText(FileStream);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator>
  void Vector<T, VectSparse, Allocator>::ReadText(istream& stream)
  {
    Clear();
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector<VectSparse>::ReadText(istream& stream)",
                    "Stream is not ready.");
#endif
 
    Vector<T, VectFull, Allocator> values;
    Vector<int> index;
    T entry;
    int ind = 0;
    int nb_elt = 0;
    while (!stream.eof())
      {
        // New entry is read.
        stream >> ind >> entry;
 
        if (stream.fail())
          break;
        else
          {
#ifdef SELDON_CHECK_IO
            if (ind < 1)
              throw IOError(string("Vector<VectSparse>::ReadText") +
                            "(istream& stream)",
                            string("Index should be greater ")
                            + "than 0 but is equal to " + to_str(ind) + ".");
#endif
 
            nb_elt++;
            ind--;
 
            // Inserting a new element.
            if (nb_elt > values.GetM())
              {
                values.Resize(2 * nb_elt);
                index.Resize(2 * nb_elt);
              }
 
            values(nb_elt - 1) = entry;
            index(nb_elt - 1) = ind;
          }
      }
 
    if (nb_elt > 0)
      {
        // Allocating to the right size.
        this->Reallocate(nb_elt);
        for (int i = 0; i < nb_elt; i++)
          {
            Index(i) = index(i);
            Value(i) = values(i);
          }
      }
  }
 
 
 
  template <class T, class Allocator>
  ostream& operator << (ostream& out,
                        const Vector<T, VectSparse, Allocator>& V)
  {
    for (size_t i = 0; i < V.GetLength(); i++)
      out << (V.Index(i) + 1) << ' ' << V.Value(i) << '\n';
 
    return out;
  }
  
} // namespace Seldon.
 
#define SELDON_FILE_SPARSE_VECTOR_CXX
#endif