array3d/Array3D.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_ARRAY3D_CXX

#include "Array3D.hxx"

namespace Seldon
{


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



  template <class T, class Allocator>
  inline Array3D<T, Allocator>::Array3D()
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;

    length23_ = 0;

    data_ = NULL;
  }



  template <class T, class Allocator>
  inline Array3D<T, Allocator>::Array3D(int i, int j, int k)
  {
    length1_ = i;
    length2_ = j;
    length3_ = k;

    length23_ = length2_ * length3_;

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        data_ = array3D_allocator_.allocate(i*j*k, this);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        length1_ = 0;
        length2_ = 0;
        length3_ = 0;
        length23_ = 0;
        data_ = NULL;
      }
    if (data_ == NULL && i != 0 && j != 0 && k != 0)
      throw NoMemory("Array3D::Array3D(int, int, int)",
                     string("Unable to allocate memory for an array of size ")
                     + to_str(static_cast<long int>(i)
                              * static_cast<long int>(j)
                              * static_cast<long int>(k)
                              * static_cast<long int>(sizeof(T)))
                     + " bytes (" + to_str(i) + " x " + to_str(j)
                     + " x " + to_str(k) + " elements).");
#endif

  }


  template <class T, class Allocator>
  Array3D<T, Allocator>::Array3D(const Array3D<T, Allocator>& A)
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;

    length23_ = 0;

    data_ = NULL;

    Copy(A);
  }


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


  template <class T, class Allocator>
  inline Array3D<T, Allocator>::~Array3D()
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;
    length23_ = 0;

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

        if (data_ != NULL)
          {
            array3D_allocator_.deallocate(data_, length1_ * length23_);
            data_ = NULL;
          }

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

  }


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



  template <class T, class Allocator>
  inline void Array3D<T, Allocator>::Reallocate(int i, int j, int k)
  {
    if (i != length1_ || j != length2_ || k != length3_)
      {
        length1_ = i;
        length2_ = j;
        length3_ = k;

        length23_ = j * k;

#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif

            data_ =
              reinterpret_cast<pointer>(array3D_allocator_.reallocate(data_,
                                                                      i*j*k,
                                                                      this));

#ifdef SELDON_CHECK_MEMORY
          }
        catch (...)
          {
            length1_ = 0;
            length2_ = 0;
            length3_ = 0;
            length23_ = 0;
            data_ = NULL;
          }
        if (data_ == NULL && i != 0 && j != 0 && k != 0)
          throw NoMemory("Array3D::Reallocate(int, int, int)",
                         string("Unable to reallocate memory")
                         + " for an array of size "
                         + to_str(static_cast<long int>(i)
                                  * static_cast<long int>(j)
                                  * static_cast<long int>(k)
                                  * static_cast<long int>(sizeof(T)))
                         + " bytes (" + to_str(i) + " x " + to_str(j)
                         + " x " + to_str(k) + " elements).");
#endif

      }
  }



  template <class T, class Allocator>
  inline void Array3D<T, Allocator>::Clear()
  {
    this->~Array3D();
    this->length1_ = 0;
    this->length2_ = 0;
    this->length3_ = 0;
  }


  /*****************
   * BASIC METHODS *
   *****************/



  template <class T, class Allocator>
  int Array3D<T, Allocator>::GetLength1() const
  {
    return length1_;
  }



  template <class T, class Allocator>
  int Array3D<T, Allocator>::GetLength2() const
  {
    return length2_;
  }



  template <class T, class Allocator>
  int Array3D<T, Allocator>::GetLength3() const
  {
    return length3_;
  }



  template <class T, class Allocator>
  int Array3D<T, Allocator>::GetSize() const
  {
    return length1_ * length23_;
  }



  template <class T, class Allocator>
  int Array3D<T, Allocator>::GetDataSize() const
  {
    return length1_ * length23_;
  }



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


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



  template <class T, class Allocator>
  inline typename Array3D<T, Allocator>::reference
  Array3D<T, Allocator>::operator() (int i, int j, int k)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= length1_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #1 should be in [0, ")
                       + to_str(length1_-1) + "], but is equal to "
                       + to_str(i) + ".");
    if (j < 0 || j >= length2_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #2 should be in [0, ")
                       + to_str(length2_-1) + "], but is equal to "
                       + to_str(j) + ".");
    if (k < 0 || k >= length3_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #3 should be in [0, ")
                       + to_str(length3_-1) + "], but is equal to "
                       + to_str(k) + ".");
#endif

    return data_[i * length23_ + j * length3_ + k];
  }



  template <class T, class Allocator>
  inline typename Array3D<T, Allocator>::const_reference
  Array3D<T, Allocator>::operator() (int i, int j, int k) const
  {

#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= length1_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #1 should be in [0, ")
                       + to_str(length1_-1) + "], but is equal to "
                       + to_str(i) + ".");
    if (j < 0 || j >= length2_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #2 should be in [0, ")
                       + to_str(length2_-1) + "], but is equal to "
                       + to_str(j) + ".");
    if (k < 0 || k >= length3_)
      throw WrongIndex("Array3D::operator()",
                       string("Index along dimension #3 should be in [0, ")
                       + to_str(length3_-1) + "], but is equal to "
                       + to_str(k) + ".");
#endif

    return data_[i*length23_ + j*length3_ + k];
  }


  template <class T, class Allocator>
  inline Array3D<T, Allocator>& Array3D<T, Allocator>::operator=
  (const Array3D<T, Allocator>& A)
  {
    Copy(A);

    return *this;
  }


  template <class T, class Allocator>
  inline void Array3D<T, Allocator>::Copy(const Array3D<T, Allocator>& A)
  {
    Reallocate(A.GetLength1(), A.GetLength2(), A.GetLength3());

    array3D_allocator_.memorycpy(data_, A.GetData(), GetDataSize());
  }


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



  template <class T, class Allocator>
  void Array3D<T, Allocator>::Zero()
  {
    array3D_allocator_.memoryset(data_, char(0),
                                 GetDataSize()*sizeof(value_type));
  }



  template <class T, class Allocator>
  void Array3D<T, Allocator>::Fill()
  {
    for (int i = 0; i < GetDataSize(); i++)
      data_[i] = i;
  }



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



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



  template <class T, class Allocator>
  void Array3D<T, Allocator>::Print() const
  {
    int i, j, k;

    for (i = 0; i < GetLength1(); i++)
      {
        for (j = 0; j < GetLength2(); j++)
          {
            for (k = 0; k < GetLength3(); k++)
              cout << (*this)(i, j, k) << '\t';
            cout << endl;
          }
        cout << endl;
      }
  }


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



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

    Write(FileStream);

    FileStream.close();
  }



  template <class T, class Allocator> void Array3D<T, Allocator>
  ::Write(ofstream& FileStream) const
  {

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

    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&length1_)),
                     sizeof(int));
    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&length2_)),
                     sizeof(int));
    FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&length3_)),
                     sizeof(int));

    FileStream.write(reinterpret_cast<char*>(data_),
                     length23_ * length1_ * sizeof(value_type));

#ifdef SELDON_CHECK_IO
    // Checks if data was written.
    if (!FileStream.good())
      throw IOError("Array3D::Write(ofstream& FileStream)",
                    string("Output operation failed.")
                    + string("  The output file may have been removed")
                    + " or there is no space left on device.");
#endif

  }



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

    Read(FileStream);

    FileStream.close();
  }



  template <class T, class Allocator>
  void Array3D<T, Allocator>
  ::Read(ifstream& FileStream)
  {

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

    int new_l1, new_l2, new_l3;
    FileStream.read(reinterpret_cast<char*>(&new_l1), sizeof(int));
    FileStream.read(reinterpret_cast<char*>(&new_l2), sizeof(int));
    FileStream.read(reinterpret_cast<char*>(&new_l3), sizeof(int));
    Reallocate(new_l1, new_l2, new_l3);

    FileStream.read(reinterpret_cast<char*>(data_),
                    length23_ * length1_ * sizeof(value_type));

#ifdef SELDON_CHECK_IO
    // Checks if data was read.
    if (!FileStream.good())
      throw IOError("Array3D::Read(ifstream& FileStream)",
                    string("Input operation failed.")
                    + string(" The input file may have been removed")
                    + " or may not contain enough data.");
#endif

  }



  template <class T, class Allocator>
  ostream& operator << (ostream& out,
                        const Array3D<T, Allocator>& A)
  {
    int i, j, k;

    for (i = 0; i < A.GetLength1(); i++)
      {
        for (j = 0; j < A.GetLength2(); j++)
          {
            for (k = 0; k < A.GetLength3(); k++)
              out << A(i, j, k) << '\t';
            out << endl;
          }
        out << endl;
      }

    return out;
  }



  template <class T0, class T, class Allocator>
  void Mlt(const T0& alpha, Array3D<T, Allocator>& A)
  {
    T* data = A.GetData();
    for (int i = 0; i < A.GetDataSize(); i++)
      data[i] *= alpha;
  }


} // namespace Seldon.

#define SELDON_FILE_ARRAY3D_CXX
#endif