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>
  Array3D<T, Allocator>::Array3D()
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;
 
    length23_ = 0;
 
    data_ = NULL;
  }
 
 
 
  template <class T, class Allocator>
  Array3D<T, Allocator>::Array3D(int i, int j, int k)
  {
    length1_ = i;
    length2_ = j;
    length3_ = k;
 
    length23_ = long(length2_) * long(length3_);
    int nb_elt = long(i)*long(j)*long(k);
 
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif
 
        data_ = Allocator::allocate(nb_elt, 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);
  }
 
 
  /*********************
   * MEMORY MANAGEMENT *
   *********************/
 
 
 
  template <class T, class Allocator>
  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_ = long(j) * long(k);
        int nb_elt = long(i)*long(j)*long(k);
 
#ifdef SELDON_CHECK_MEMORY
        try
          {
#endif
            
            data_ =
              reinterpret_cast<pointer>(Allocator::reallocate(data_,
                                                              nb_elt, 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>
  void Array3D<T, Allocator>::SetData(int i, int j, int k, 
                                      typename Array3D<T, Allocator>
                                             ::pointer data)
  {
    Clear();
    
    length1_ = i;
    length2_ = j;
    length3_ = k;
    length23_ = long(j) * long(k);
    data_ = data;
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>::Nullify()
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;
    length23_ = 0;
    data_ = NULL;
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>::Clear()
  {
    length1_ = 0;
    length2_ = 0;
    length3_ = 0;
    length23_ = 0;
 
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif
 
        if (data_ != NULL)
          {
            Allocator::deallocate(data_, length1_ * length23_);
            data_ = NULL;
          }
 
#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
        data_ = NULL;
      }
#endif
 
    this->length1_ = 0;
    this->length2_ = 0;
    this->length3_ = 0;
  }
 
 
  /************************
   * CONVENIENT FUNCTIONS *
   ************************/
 
 
 
  template <class T, class Allocator>
  size_t Array3D<T, Allocator>::GetMemorySize() const
  {
    return sizeof(*this) + size_t(sizeof(T))*GetDataSize();
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>::Zero()
  {
    Allocator::memoryset(data_, char(0),
                         GetDataSize()*sizeof(value_type));
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>::Fill()
  {
    long taille = this->GetDataSize();
    for (long i = 0; i < taille; i++)
      SetComplexReal(i, data_[i]);
  }
 
 
 
  template <class T, class Allocator>
  template <class T0>
  void Array3D<T, Allocator>::Fill(const T0& x)
  {
    long taille = this->GetDataSize();
    T x_;
    SetComplexReal(x, x_);
    for (long i = 0; i < taille; i++)
      data_[i] = x_;
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>::FillRand()
  {
#ifndef SELDON_WITHOUT_REINIT_RANDOM
    srand(time(NULL));
#endif
    long taille = this->GetDataSize();
    for (long i = 0; i < taille; i++)
      SetComplexReal(rand(), this->data_[i]);
  }
 
 
 
  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, 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("Array3D::Write(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    Write(FileStream, with_size);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator> void Array3D<T, Allocator>
  ::Write(ofstream& FileStream, bool with_size) 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
 
    if (with_size)
      {
        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, 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("Array3D::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
 
    Read(FileStream, with_size);
 
    FileStream.close();
  }
 
 
 
  template <class T, class Allocator>
  void Array3D<T, Allocator>
  ::Read(ifstream& FileStream, bool with_size)
  {
 
#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
 
    if (with_size)
      {
        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("Output operation failed.")
                    + string(" The intput 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 MltScalar(const T0& alpha, Array3D<T, Allocator>& A)
  {
    T* data = A.GetData();
    long taille = A.GetDataSize();
    for (long i = 0; i < taille; i++)
      data[i] *= alpha;
  }
  
} // namespace Seldon.
 
#define SELDON_FILE_ARRAY3D_CXX
#endif