Vector3.cxx

// Copyright (C) 2010-2012, INRIA
// Author(s): Marc Fragu, 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_VECTOR_VECTOR_3_CXX
 
 
#include "Vector3.hxx"
 
 
namespace Seldon
{
 
 
  // VECTOR3 //
 
 
  /***************
   * CONSTRUCTOR *
   ***************/
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector3<T, Allocator0, Allocator1, Allocator2>::Vector3()
  {
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector3<T, Allocator0, Allocator1, Allocator2>::Vector3(int length)
  {
    data_.Reallocate(length);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Vector3(Vector<int>& length)
  {
    data_.Clear();
    int n, m = length.GetSize();
    data_.Reallocate(m);
    for (int i = 0; i < m; i++)
      {
        n = length(i);
        data_(i).Reallocate(n);
      }
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  template <class Allocator>
  Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Vector3(Vector<Vector<int>, Vect_Full, Allocator>& length)
  {
    data_.Clear();
    int n, m = length.GetSize();
    data_.Reallocate(m);
    for (int i = 0; i < m; i++)
      {
        n = length(i).GetSize();
        data_(i).Reallocate(n);
        for (int j = 0; j < n; j++)
          data_(i)(j).Reallocate(length(i)(j));
      }
  }
 
 
  /**************
   * DESTRUCTOR *
   **************/
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector3<T, Allocator0, Allocator1, Allocator2>::~Vector3()
  {
  }
 
 
  /**********************
   * VECTORS MANAGEMENT *
   **********************/
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>::GetLength() const
  {
    return data_.GetLength();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>::GetSize() const
  {
    return data_.GetSize();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>::GetLength(int i) const
  {
    return data_(i).GetLength();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>::GetSize(int i) const
  {
    return data_(i).GetSize();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetLength(int i, int j) const
  {
    return data_(i)(j).GetLength();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetSize(int i, int j) const
  {
    return data_(i)(j).GetSize();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  size_t Vector3<T, Allocator0, Allocator1, Allocator2>::GetMemorySize() const
  {
    size_t total = sizeof(*this) + sizeof(pointer)*GetLength();
    for (int i = 0; i < GetLength(); i++)
      {
        total += sizeof(int)+sizeof(pointer)+sizeof(pointer)*GetLength(i);
        for (int j = 0; j < GetLength(i); j++)
          total += data_(i)(j).GetMemorySize();
      }
    
    return total;
  }
  
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>::GetNelement() const
  {
    int total = 0;
    for (int i = 0; i < GetLength(); i++)
      for (int j = 0; j < GetLength(i); j++ )
        total += GetLength(i, j);
    return total;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetNelement(int beg, int end) const
  {
    if (beg > end)
      throw WrongArgument("Vector3::GetNelement(int beg, int end)",
                          "The lower bound of the range of inner vectors "
                          "of vectors, [" + to_str(beg) + ", " + to_str(end)
                          + "[, is strictly greater than its upper bound.");
    if (beg < 0 || end > GetLength())
      throw WrongArgument("Vector3::GetNelement(int beg, int end)",
                          "The inner-vector of vectors indexes should be in "
                          "[0," + to_str(GetLength()) + "] but ["
                          + to_str(beg)
                          + ", " + to_str(end) + "[ was provided.");
 
    int total = 0;
    for (int i = beg; i < end; i++)
      for (int j = 0; j < GetLength(i); j++ )
        total += GetLength(i, j);
    return total;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  int Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetNelement(int beg0, int end0, int beg1, int end1) const
  {
    if (beg0 > end0)
      throw WrongArgument("Vector3:::GetNelement(int beg0, int end0, "
                          "int beg1, int end1)",
                          "The lower bound of the range of inner vectors "
                          "of vectors, [" + to_str(beg0) + ", " + to_str(end0)
                          + "[, is strictly greater than its upper bound.");
    if (beg1 > end1)
      throw WrongArgument("Vector3:::GetNelement(int beg0, int end0, "
                          "int beg1, int end1)",
                          "The lower bound of the range of "
                          "of vectors, [" + to_str(beg1) + ", " + to_str(end1)
                          + "[, is strictly greater than its upper bound.");
    if (beg0 < 0 || end0 > GetLength())
      throw WrongArgument("Vector3:::GetNelement(int beg0, int end0, "
                          "int beg1, int end1)",
                          "The inner-vector of vectors indexes should be in "
                          "[0," + to_str(GetLength()) + "] but ["
                          + to_str(beg0)
                          + ", " + to_str(end0) + "[ was provided.");
 
    int total = 0;
    for (int i = beg0; i < end0; i++)
      {
        if (beg1 < 0 || end1 > GetLength(i))
          throw WrongArgument("Vector3:::GetNelement(int beg0, int end0, "
                              "int beg1, int end1)",
                              "For the inner vector of vectors "
                              + to_str(i) + ", the vectors indexes should be "
                              "in [0," + to_str(GetLength(i)) + "] but ["
                              + to_str(beg1)
                              + ", " + to_str(end1) + "[ was provided.");
        for (int j = beg1; j < end1; j++)
          total += GetLength(i, j);
      }
    return total;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<int> Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetShape(int i) const
  {
    Vector<int> shape(GetLength(i));
    for (int j = 0; j < GetLength(i); j++)
      shape(j) = GetLength(i, j);
    return shape;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::GetShape(int i, Vector<int>& shape) const
  {
    shape.Reallocate(GetLength(i));
    for (int j = 0; j < GetLength(i); j++)
      shape(j) = GetLength(i, j);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::Reallocate(int N)
  {
    data_.Reallocate(N);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Reallocate(int i, int N)
  {
    data_(i).Reallocate(N);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Reallocate(int i, int j, int N)
  {
    data_(i)(j).Reallocate(N);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  template <class Td, class Allocatord>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Flatten(Vector<Td, VectFull, Allocatord>& data) const
  {
    data.Reallocate(GetNelement());
    int i, j, k, n(0);
    for (i = 0; i < GetLength(); i++)
      for (j = 0; j < GetLength(i); j++)
        for (k = 0; k < GetLength(i, j); k++)
          data(n++) = data_(i)(j)(k);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  template <class Td, class Allocatord>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Flatten(int beg, int end, Vector<Td, VectFull, Allocatord>& data) const
  {
    if (beg > end)
      throw WrongArgument("Vector3:::Flatten(int beg, int end, Vector& data)",
                          "The lower bound of the range of inner vectors "
                          "of vectors, [" + to_str(beg) + ", " + to_str(end)
                          + "[, is strictly greater than its upper bound.");
    if (beg < 0 || end > GetLength())
      throw WrongArgument("Vector3:::Flatten(int beg, int end, Vector& data)",
                          "The inner-vector of vectors indexes should be in "
                          "[0," + to_str(GetLength()) + "] but ["
                          + to_str(beg)
                          + ", " + to_str(end) + "[ was provided.");
 
    data.Reallocate(GetNelement(beg, end));
    int i, j, k, n(0);
    for (i = beg; i < end; i++)
      for (j = 0; j < GetLength(i); j++)
        for (k = 0; k < GetLength(i, j); k++)
          data(n++) = data_(i)(j)(k);
  }
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  template <class Td, class Allocatord>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Flatten(int beg0, int end0, int beg1, int end1,
            Vector<Td, VectFull, Allocatord>& data) const
  {
    if (beg0 > end0)
      throw WrongArgument("Vector3:::Flatten(int beg0, int end0, int beg1, "
                          "int end1, Vector& data)",
                          "The lower bound of the range of inner vectors "
                          "of vectors, [" + to_str(beg0) + ", " + to_str(end0)
                          + "[, is strictly greater than its upper bound.");
    if (beg1 > end1)
      throw WrongArgument("Vector3:::Flatten(int beg0, int end0, int beg1, "
                          "int end1, Vector& data)",
                          "The lower bound of the range of "
                          "of vectors, [" + to_str(beg1) + ", " + to_str(end1)
                          + "[, is strictly greater than its upper bound.");
    if (beg0 < 0 || end0 > GetLength())
      throw WrongArgument("Vector3:::Flatten(int beg0, int end0, int beg1, "
                          "int end1, Vector& data)",
                          "The inner-vector of vectors indexes should be in "
                          "[0," + to_str(GetLength()) + "] but ["
                          + to_str(beg0)
                          + ", " + to_str(end0) + "[ was provided.");
 
    data.Reallocate(GetNelement(beg0, end0, beg1, end1));
    int i, j, k, n(0);
    for (i = beg0; i < end0; i++)
      {
        if (beg1 < 0 || end1 > GetLength(i))
          throw WrongArgument("Vector3:::Flatten(int beg0, int end0, int beg1, "
                              "int end1, Vector& data)",
                              "For the inner vector of vectors "
                              + to_str(i) + ", the vectors indexes should be "
                              "in [0," + to_str(GetLength(i)) + "] but ["
                              + to_str(beg1)
                              + ", " + to_str(end1) + "[ was provided.");
        for (j = beg1; j < end1; j++)
          for (k = 0; k < GetLength(i, j); k++)
            data(n++) = data_(i)(j)(k);
      }
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0,
               Allocator1, Allocator2>::PushBack(int i, int j, const T& x)
  {
    data_(i)(j).PushBack(x);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::
  PushBack(int i, const Vector<T, Vect_Full, Allocator0>& X)
  {
    data_(i).PushBack(X);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::PushBack(const Vector<Vector<T, Vect_Full, Allocator0>,
             Vect_Full, Allocator1>& X)
  {
    data_.PushBack(X);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::PushBack(const Vector<Vector<Vector<T, Vect_Full, Allocator0>,
             Vect_Full, Allocator1>, Vect_Full, Allocator2>& X)
  {
    for (int i = 0; i < X.GetM(); i++)
      data_.PushBack(X(i));
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::PushBack(const Vector3<T, Allocator0, Allocator1, Allocator2>& X)
  {
    for (int i = 0; i < X.GetLength(); i++)
      data_.PushBack(X.GetVector());
  }
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::Clear()
  {
    data_.Clear();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0,
               Allocator1, Allocator2>::Clear(int i)
  {
    data_(i).Clear();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::Clear(int i, int j)
  {
    data_(i)(j).Clear();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::Fill(const T& x)
  {
    for (int i = 0; i < data_.GetM(); i++)
      for (int j = 0; j < data_(i).GetM(); j++)
        data_(i)(j).Fill(x);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<Vector<Vector<T, Vect_Full, Allocator0>, Vect_Full, Allocator1>,
         Vect_Full, Allocator2>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector()
  {
    return data_;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  const Vector<Vector<Vector<T, Vect_Full, Allocator0>,
                      Vect_Full, Allocator1>, Vect_Full, Allocator2>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector() const
  {
    return data_;
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<Vector<T, Vect_Full, Allocator0>, VectFull, Allocator1>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector(int i)
  {
    return data_(i);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  const Vector<Vector<T, Vect_Full, Allocator0>, VectFull, Allocator1>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector(int i) const
  {
    return data_(i);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<T, Vect_Full, Allocator0>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector(int i, int j)
  {
    return data_(i)(j);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  const Vector<T, Vect_Full, Allocator0>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::GetVector(int i, int j)
    const
  {
    return data_(i)(j);
  }
 
 
  /*********************************
   * ELEMENT ACCESS AND ASSIGNMENT *
   *********************************/
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  const Vector<Vector<T, Vect_Full, Allocator0>, VectFull, Allocator1>&
  Vector3<T, Allocator0,
          Allocator1, Allocator2>::operator() (int i) const
  {
    return data_(i);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<Vector<T, Vect_Full, Allocator0>, VectFull, Allocator1>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::operator() (int i)
  {
    return data_(i);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  const Vector<T, Vect_Full, Allocator0>&
  Vector3<T, Allocator0,
          Allocator1, Allocator2>::operator() (int i, int j) const
  {
    return data_(i)(j);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  Vector<T, Vect_Full, Allocator0>&
  Vector3<T, Allocator0, Allocator1, Allocator2>::operator() (int i, int j)
  {
    return data_(i)(j);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  typename Vector3<T, Allocator0, Allocator1, Allocator2>::const_reference
  Vector3<T, Allocator0,
          Allocator1, Allocator2>::operator() (int i, int j, int k) const
  {
    return data_(i)(j)(k);
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  typename Vector3<T, Allocator0, Allocator1, Allocator2>::reference
  Vector3<T, Allocator0, Allocator1, Allocator2>::operator()
    (int i, int j, int k)
  {
    return data_(i)(j)(k);
  }
 
 
  /**********************
   * CONVENIENT METHODS *
   *********************/
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>::Print() const
  {
    for (int i = 0; i < data_.GetM(); i++)
      for(int j = 0; j < data_(i).GetM(); j++)
        {
          cout << "Vector " << i << ", " << j << ": ";
          data_(i)(j).Print();
        }
  }
 
 
  /**************************
   * INPUT/OUTPUT FUNCTIONS *
   **************************/
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Write(string file_name, bool with_size) const
  {
    ofstream file_stream;
    file_stream.open(file_name.c_str());
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!file_stream.is_open())
      throw IOError("Vector3::Write(string file_name, bool with_size)",
                    string("Unable to open file \"") + file_name + "\".");
#endif
 
    this->Write(file_stream, with_size);
 
    file_stream.close();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Write(ostream& stream, bool with_size) const
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector3::Write(ostream& stream, bool with_size)",
                    "The stream is not ready.");
#endif
 
    if (with_size)
      {
        int m = GetLength();
        stream.write(reinterpret_cast<char*>(const_cast<int*>(&m)),
                     sizeof(int));
      }
 
    for (int i = 0; i < GetLength(); i++)
      {
        if (with_size)
          {
            int m = GetLength(i);
            stream.write(reinterpret_cast<char*>(const_cast<int*>(&m)),
                         sizeof(int));
          }
        for (int j = 0; j < GetLength(i); j++)
          data_(i)(j).Write(stream, with_size);
      }
 
#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!stream.good())
      throw IOError("Vector3::Write(ostream& stream, bool with_size)",
                    "Output operation failed.");
#endif
 
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Read(string file_name, bool with_size)
  {
    ifstream file_stream;
    file_stream.open(file_name.c_str());
 
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!file_stream.is_open())
      throw IOError("Vector3::Read(string file_name, bool with_size)",
                    string("Unable to open file \"") + file_name + "\".");
#endif
 
    this->Read(file_stream, with_size);
 
    file_stream.close();
  }
 
 
 
  template <class T, class Allocator0, class Allocator1, class Allocator2>
  void Vector3<T, Allocator0, Allocator1, Allocator2>
  ::Read(istream& stream, bool with_size)
  {
 
#ifdef SELDON_CHECK_IO
    // Checks if the stream is ready.
    if (!stream.good())
      throw IOError("Vector3::Read(istream& stream, bool with_size)",
                    "The stream is not ready.");
#endif
 
    if (with_size)
      {
        int new_size;
        stream.read(reinterpret_cast<char*>(&new_size), sizeof(int));
        this->Reallocate(new_size);
      }
 
    for (int i = 0; i < GetLength(); i++)
      {
        if (with_size)
          {
            int new_size;
            stream.read(reinterpret_cast<char*>(&new_size), sizeof(int));
            this->Reallocate(i, new_size);
          }
        for (int j = 0; j < GetLength(i); j++)
          data_(i)(j).Read(stream, with_size);
      }
 
#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!stream.good())
      throw IOError("Vector3::Read(istream& stream, bool with_size)",
                    "Output operation failed.");
#endif
 
  }
 
 
} // namespace Seldon.
 
 
#define SELDON_FILE_VECTOR_VECTOR_3_CXX
#endif