vector/Vector.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_VECTOR_CXX

#include "Vector.hxx"

namespace Seldon
{


  // C_VECTOR_BASE //


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



  template <class T, class Allocator>
  inline Vector_Base<T, Allocator>::Vector_Base()
  {
    m_ = 0;
    data_ = NULL;
  }



  template <class T, class Allocator>
  inline Vector_Base<T, Allocator>::Vector_Base(int i)
  {
    m_ = i;
    data_ = NULL;
  }



  template <class T, class Allocator>
  inline Vector_Base<T, Allocator>::
  Vector_Base(const Vector_Base<T, Allocator>& A)
  {
    m_ = A.GetM();
    data_ = NULL;
  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Allocator>
  inline Vector_Base<T, Allocator>::~Vector_Base()
  {

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        if (data_ != NULL)
          {
            vect_allocator_.deallocate(data_, m_);
            m_ = 0;
            data_ = NULL;
          }

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        m_ = 0;
        data_ = NULL;
      }
#endif

  }


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



  template <class T, class Allocator>
  int Vector_Base<T, Allocator>::GetM() const
  {
    return m_;
  }



  template <class T, class Allocator>
  int Vector_Base<T, Allocator>::GetLength() const
  {
    return m_;
  }



  template <class T, class Allocator>
  int Vector_Base<T, Allocator>::GetSize() const
  {
    return m_;
  }



  template <class T, class Allocator>
  typename Vector_Base<T, Allocator>::pointer
  Vector_Base<T, Allocator>::GetData() const
  {
    return data_;
  }



  template <class T, class Allocator>
  typename Vector_Base<T, Allocator>::const_pointer
  Vector_Base<T, Allocator>::GetDataConst() const
  {
    return reinterpret_cast<typename Vector_Base<T,
      Allocator>::const_pointer>(data_);
  }



  template <class T, class Allocator>
  void* Vector_Base<T, Allocator>::GetDataVoid() const
  {
    return reinterpret_cast<void*>(data_);
  }



  template <class T, class Allocator>
  const void* Vector_Base<T, Allocator>::GetDataConstVoid() const
  {
    return reinterpret_cast<const void*>(data_);
  }


  // VECTOR<VECT_FULL> //


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



  template <class T, class Allocator>
  Vector<T, VectFull, Allocator>::Vector():
    Vector_Base<T, Allocator>()
  {
  }



  template <class T, class Allocator>
  Vector<T, VectFull, Allocator>::Vector(int i):
    Vector_Base<T, Allocator>(i)
  {

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        this->data_ = this->vect_allocator_.allocate(i, this);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->data_ = NULL;
      }
    if (this->data_ == NULL)
      this->m_ = 0;
    if (this->data_ == NULL && i != 0)
      throw NoMemory("Vector<VectFull>::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>
  Vector<T, VectFull, Allocator>
  ::Vector(int i, typename Vector<T, VectFull, Allocator>::pointer data):
    Vector_Base<T, Allocator>()
  {
    SetData(i, data);
  }



  template <class T, class Allocator>
  Vector<T, VectFull, Allocator>::
  Vector(const Vector<T, VectFull, Allocator>& V):
    Vector_Base<T, Allocator>(V)
  {

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        this->data_ = this->vect_allocator_.allocate(V.GetM(), this);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        this->m_ = 0;
        this->data_ = NULL;
      }
    if (this->data_ == NULL)
      this->m_ = 0;
    if (this->data_ == NULL && V.GetM() != 0)
      throw NoMemory("Vector<VectFull>::Vector(Vector<VectFull>&)",
                     string("Unable to allocate memory for a vector of size ")
                     + to_str(V.GetM()*sizeof(T)) + " bytes ("
                     + to_str(V.GetM()) + " elements).");
#endif

    this->vect_allocator_.memorycpy(this->data_, V.GetData(), V.GetM());

  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Allocator>
  Vector<T, VectFull, Allocator>::~Vector()
  {
  }


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



  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>::Clear()
  {
    this->~Vector();
  }



  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>::Reallocate(int i)
  {

    if (i != this->m_)
      {

        this->m_ = i;

#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif

            this->data_ =
              reinterpret_cast<pointer>(this->vect_allocator_
                                        .reallocate(this->data_, i, this));

#ifdef SELDON_CHECK_MEMORY
          }
        catch (...)
          {
            this->m_ = 0;
            this->data_ = NULL;
            return;
          }
        if (this->data_ == NULL)
          {
            this->m_ = 0;
            return;
          }
#endif

      }
  }



  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>::Resize(int n)
  {

    if (n == this->m_)
      return;

    Vector<T, VectFull, Allocator> X_new(n);
    for (int i = 0; i < min(this->m_, n); i++)
      X_new(i) = this->data_[i];

    SetData(n, X_new.GetData());
    X_new.Nullify();
  }


  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>
  ::SetData(int i, typename Vector<T, VectFull, Allocator>::pointer data)
  {
    this->Clear();

    this->m_ = i;

    this->data_ = data;
  }



  template <class T, class Allocator>
  template <class Allocator0>
  inline void Vector<T, VectFull, Allocator>
  ::SetData(const Vector<T, VectFull, Allocator0>& V)
  {
    SetData(V.GetLength(), V.GetData());
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::Nullify()
  {
    this->m_ = 0;
    this->data_ = NULL;
  }


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



  template <class T, class Allocator>
  inline typename Vector<T, VectFull, Allocator>::reference
  Vector<T, VectFull, Allocator>::operator() (int i)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongIndex("Vector<VectFull>::operator()",
                       string("Index should be in [0, ") + to_str(this->m_-1)
                       + "], but is equal to " + to_str(i) + ".");
#endif

    return this->data_[i];
  }



  template <class T, class Allocator>
  inline typename Vector<T, VectFull, Allocator>::reference
  Vector<T, VectFull, Allocator>::Get(int i)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongIndex("Vector<VectFull>::Get(int)",
                       string("Index should be in [0, ") + to_str(this->m_-1)
                       + "], but is equal to " + to_str(i) + ".");
#endif

    return this->data_[i];
  }



  template <class T, class Allocator>
  inline typename Vector<T, VectFull, Allocator>::const_reference
  Vector<T, VectFull, Allocator>::operator() (int i) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongIndex("Vector<VectFull>::operator() const",
                       string("Index should be in [0, ") + to_str(this->m_-1)
                       + "], but is equal to " + to_str(i) + ".");
#endif

    return this->data_[i];
  }



  template <class T, class Allocator>
  inline typename Vector<T, VectFull, Allocator>::const_reference
  Vector<T, VectFull, Allocator>::Get(int i) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongIndex("Vector<VectFull>::Get(int) const",
                       string("Index should be in [0, ") + to_str(this->m_-1)
                       + "], but is equal to " + to_str(i) + ".");
#endif

    return this->data_[i];
  }



  template <class T, class Allocator>
  inline Vector<T, VectFull, Allocator>& Vector<T, VectFull, Allocator>
  ::operator= (const Vector<T, VectFull, Allocator>& X)
  {
    this->Copy(X);

    return *this;
  }



  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>
  ::Copy(const Vector<T, VectFull, Allocator>& X)
  {
    this->Reallocate(X.GetLength());

    this->vect_allocator_.memorycpy(this->data_, X.GetData(), this->m_);
  }



  template <class T, class Allocator>
  inline Vector<T, VectFull, Allocator>
  Vector<T, VectFull, Allocator>::Copy() const
  {
    return Vector<T, VectFull, Allocator>(*this);
  }



  template <class T, class Allocator> template<class T0>
  inline Vector<T, VectFull, Allocator>& Vector<T, VectFull, Allocator>
  ::operator*= (const T0& alpha)
  {
    for (int i = 0; i < this->m_; i++)
      this->data_[i] *= alpha;

    return *this;
  }



  template <class T, class Allocator>
  inline void Vector<T, VectFull, Allocator>::Append(const T& x)
  {
    int i = this->GetLength();
    this->Reallocate(i + 1);
    this->data_[i] = x;
  }



  template <class T, class Allocator> template<class T0>
  inline void Vector<T, VectFull, Allocator>::PushBack(const T0& x)
  {
    Resize(this->m_+1);
    this->data_[this->m_-1] = x;
  }



  template <class T, class Allocator> template<class Allocator0>
  inline void Vector<T, VectFull, Allocator>
  ::PushBack(const Vector<T, VectFull, Allocator0>& X)
  {
    int Nold = this->m_;
    Resize(this->m_ + X.GetM());
    for (int i = 0; i < X.GetM(); i++)
      this->data_[Nold+i] = X(i);
  }


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



  template <class T, class Allocator>
  int Vector<T, VectFull, Allocator>::GetDataSize()
  {
    return this->m_;
  }


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



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::Zero()
  {
    this->vect_allocator_.memoryset(this->data_, char(0),
                                    this->GetDataSize() * sizeof(value_type));
  }


  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::Fill()
  {
    for (int i = 0; i < this->m_; i++)
      this->data_[i] = i;
  }



  template <class T, class Allocator>
  template <class T0>
  void Vector<T, VectFull, Allocator>::Fill(const T0& x)
  {
    for (int i = 0; i < this->m_; i++)
      this->data_[i] = x;
  }



  template <class T, class Allocator>
  template <class T0>
  Vector<T, VectFull, Allocator>&
  Vector<T, VectFull, Allocator>::operator= (const T0& x)
  {
    this->Fill(x);

    return *this;
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::FillRand()
  {
    srand(time(NULL));
    for (int i = 0; i < this->m_; i++)
      this->data_[i] = rand();
  }


  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::Print() const
  {
    for (int i = 0; i < this->GetLength(); i++)
      cout << (*this)(i) << "\t";
    cout << endl;
  }


  /*********
   * NORMS *
   *********/



  template <class T, class Allocator>
  typename Vector<T, VectFull, Allocator>::value_type
  Vector<T, VectFull, Allocator>::GetNormInf() const
  {
    value_type res = value_type(0);
    for (int i = 0; i < this->GetLength(); i++)
      {
        res = max(res, this->data_[i]);
        res = max(res, T(-this->data_[i]));
      }

    return res;
  }



  template <class T, class Allocator>
  int Vector<T, VectFull, Allocator>::GetNormInfIndex() const
  {

#ifdef SELDON_CHECK_DIMENSIONS
    if (this->GetLength() == 0)
      throw WrongDim("Vector<VectFull>::GetNormInfIndex()",
                     "Vector is null.");
#endif

    value_type res = value_type(0), temp;
    int j = 0;
    for (int i = 0; i < this->GetLength(); i++)
      {
        temp = res;
        res = max(res, this->data_[i]);
        res = max(res, T(-this->data_[i]));
        if (temp != res) j = i;
      }

    return j;
  }


  /**************************
   * OUTPUT/INPUT FUNCTIONS *
   **************************/



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>
  ::Write(string FileName, bool with_size) 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<VectFull>::Write(string FileName, "
                    "bool with_size)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->Write(FileStream, with_size);

    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>
  ::Write(ostream& FileStream, bool with_size) const
  {

#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::Write(ostream& FileStream, "
                    "bool with_size)",
                    "The stream is not ready.");
#endif

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

    FileStream.write(reinterpret_cast<char*>(this->data_),
                     this->m_ * sizeof(value_type));

#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::Write(ostream& FileStream, "
                    "bool with_size)",
                    "Output operation failed.");
#endif

  }



  template <class T, class Allocator>
  void Vector<T, VectFull, 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<VectFull>::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->WriteText(FileStream);

    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::WriteText(ostream& FileStream) const
  {

#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::WriteText(ostream& FileStream)",
                    "The stream is not ready.");
#endif

    if (this->GetLength() != 0)
      FileStream << (*this)(0);

    for (int i = 1; i < this->GetLength(); i++)
      FileStream << "\t" << (*this)(i);

#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::WriteText(ostream& FileStream)",
                    "Output operation failed.");
#endif

  }


#ifdef SELDON_WITH_HDF5


  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::WriteHDF5(string FileName,
                                                 string group_name,
                                                 string dataset_name) const
  {
    hid_t file_id = H5Fopen(FileName.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!file_id)
      throw IOError("Vector<VectFull>::WriteHDF5(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    hid_t dataset_id, dataspace_id, group_id;
    herr_t status;

    T x(0);
    hid_t datatype = GetH5Type(x);

    if (!H5Lexists(file_id, group_name.c_str(), H5P_DEFAULT))
      group_id = H5Gcreate(file_id, group_name.c_str(), 0);
    group_id = H5Gopen(file_id, group_name.c_str());

    if (!H5Lexists(group_id, dataset_name.c_str(), H5P_DEFAULT))
      {
        // Create the initial dataspace.
        hsize_t dim[1] = {0};
        hsize_t maxdims[1] = {H5S_UNLIMITED};
        dataspace_id = H5Screate_simple(1, dim, maxdims);

        // Define chunking parameters to allow extension of the dataspace.
        hid_t cparms = H5Pcreate(H5P_DATASET_CREATE);
        hsize_t chunk_dims[1] = {this->GetLength()};
        status = H5Pset_chunk(cparms, 1, chunk_dims);

        // Create the dataset.
        hid_t filetype = H5Tvlen_create(datatype);
        dataset_id = H5Dcreate(group_id, dataset_name.c_str(),
                               filetype, dataspace_id, cparms);
      }

    // Opens the dataset, and extend it to store a new vector.
    dataset_id = H5Dopen(group_id, dataset_name.c_str());
    dataspace_id = H5Dget_space(dataset_id);
    hsize_t dims_out[1];
    status  = H5Sget_simple_extent_dims(dataspace_id, dims_out, NULL);
    hsize_t new_dim[1]= {dims_out[0] + 1};
    status = H5Dextend(dataset_id, new_dim);

    // Selects the last memory part of the dataset, to store the vector.
    hsize_t offset[1] = {dims_out[0]};
    hsize_t dim2[1] = {1};
    dataspace_id = H5Dget_space(dataset_id);
    status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
                                 offset, NULL,
                                 dim2, NULL);

    // Defines memory space.
    hid_t dataspace = H5Screate_simple(1, dim2, NULL);

    // Writes the data to the memory hyperslab selected.
    hid_t memtype = H5Tvlen_create(datatype);
    hvl_t wdata[1];
    wdata[0].len = this->GetLength();
    wdata[0].p = this->GetDataVoid();
    status = H5Dwrite(dataset_id, memtype, dataspace,
                      dataspace_id, H5P_DEFAULT, wdata);

    // Closes the dataset, group and file.
    status = H5Dclose(dataset_id);
    status = H5Gclose(group_id);
    status = H5Fclose(file_id);
  }
#endif



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>
  ::Read(string FileName, bool with_size)
  {
    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<VectFull>::Read(string FileName, bool with_size)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->Read(FileStream, with_size);

    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>
  ::Read(istream& FileStream, bool with_size)
  {

#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::Read(istream& FileStream, "
                    "bool with_size)",
                    "The stream is not ready.");
#endif

    if (with_size)
      {
        int new_size;
        FileStream.read(reinterpret_cast<char*>(&new_size), sizeof(int));
        this->Reallocate(new_size);
      }

    FileStream.read(reinterpret_cast<char*>(this->data_),
                    this->GetLength() * sizeof(value_type));

#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::Read(istream& FileStream, "
                    "bool with_size)",
                    "Output operation failed.");
#endif

  }



  template <class T, class Allocator>
  void Vector<T, VectFull, 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<VectFull>::ReadText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->ReadText(FileStream);

    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, VectFull, Allocator>::ReadText(istream& FileStream)
  {
    // Previous values of the vector are cleared.
    Clear();

#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!FileStream.good())
      throw IOError("Vector<VectFull>::ReadText(istream& FileStream)",
                    "The stream is not ready.");
#endif

    T entry;
    int number_element = 0;
    while (!FileStream.eof())
      {
        // Reads a new entry.
        FileStream >> entry;

        if (FileStream.fail())
          break;
        else
          {
            number_element++;

            // If needed, resizes the vector. Its size is already doubled so
            // that the vector should be resized a limited number of times.
            if (number_element > this->m_)
              this->Resize(2 * number_element);

            this->data_[number_element - 1] = entry;
          }
      }

    // Resizes to the actual size.
    if (number_element > 0)
      this->Resize(number_element);
    else
      this->Clear();
  }



  template <class T, class Storage, class Allocator>
  ostream& operator << (ostream& out,
                        const Vector<T, Storage, Allocator>& V)
  {
    for (int i = 0; i < V.GetLength() - 1; i++)
      out << V(i) << '\t';
    if (V.GetLength() != 0)
      out << V(V.GetLength() - 1);

    return out;
  }


} // namespace Seldon.

#define SELDON_FILE_VECTOR_CXX
#endif