// -----------------------------------------------------------
// SPL (Solitonic POOMA-based Library).
// Copyright (C) 2000 by Serge Mingaleev <address@hidden>
// Last modified: Sep 17, 2000.
// -----------------------------------------------------------

#ifndef SPL_ENGINE_SPARSEENGINE_H
#error "Should be included from SparseEngine.h only!"
#endif

#ifndef SPL_ENGINE_SPARSEENGINE4_H
#define SPL_ENGINE_SPARSEENGINE4_H

//=============================================================
//   Realization of the Engine<4,T,Sparse>
//=============================================================

#define Dim 4

template<class T>
class Engine<4,T,Sparse>
{
public:

  //============================================================
  // Exported typedefs and constants
  //============================================================

  typedef Engine<Dim,T,Sparse>             This_t;
  typedef Engine<Dim,T,Sparse>             Engine_t;
//  typedef Pooma::SparseBase<Dim>          Base_t;
//  typedef typename Base_t::Domain_t        Domain_t;
  typedef Interval<Dim>                    Domain_t;
  typedef DomainLayout<Dim>                Layout_t;
  typedef T                                Element_t;
  typedef T&                               ElementRef_t;
  typedef Sparse                           Tag_t;

  enum { sparse        = true  };
  enum { brick         = false };
  enum { dimensions    = Dim   };
  enum { hasDataObject = false };
  enum { dynamic       = false };
  enum { zeroBased     = false };
  enum { multiPatch    = false };

  //==========================================================
  // Constructors and Factory Methods
  //==========================================================

  // Default constructor. Creates a Brick-Engine with no data
  // and an "empty" domain.  This is not really usable until it
  // has been initialized (via operator=) to a new engine with an
  // actual domain.

  Engine() { }

  // These constructors take an Interval<Dim> and create a new
  // Brick-Engine with data of type T on this Domain. This is where
  // storage gets allocated. The second version uses a model data
  // element to initialize storage.
  
  Engine(const Domain_t& dom)
  {
    data_m.initialize(dom[0], dom[1]);
    for (int i=dom[0].first(); i<=dom[0].last(); i++) 
      for (int j=dom[1].first(); j<=dom[1].last(); j++) 
         data_m(i,j).initialize(dom[2], dom[3]);

    for (int d = 0; d < Dim; d++) {
      domain_m[d] = dom[d];
      firsts_m[d] = dom[d].first();
    }

    tol_m=data_m(firsts_m[0], firsts_m[1]).engine().tolerance();

    size_m=dom[0].size()*dom[1].size()
          *data_m(firsts_m[0], firsts_m[1]).engine().size();

#ifdef SPL_DEBUG_SPARSE
    cerr << "Creating 4D Sparse Array ("
         << domain_m[0] << "x" 
         << domain_m[1] << "x" 
         << domain_m[2] << "x" 
         << domain_m[3] 
         << ") with size = " << size_m 
         << endl;
#endif
  }

  // Subsetting Constructors.
  //   There are none - you cannot create a Brick-Engine by taking
  //   a "view" of another Engine.  Brick-Engines, by definition,
  //   view all of the data. See Engine<D,T,BrickView> below.

  //============================================================
  // Destructor
  //============================================================

  ~Engine(){ 
#ifdef SPL_DEBUG_SPARSE
    cerr << "Deleting 4D Sparse Array ("
         << domain_m[0] << "x" 
         << domain_m[1] << "x" 
         << domain_m[2] << "x" 
         << domain_m[3] 
         << ") with size = " << size_m 
         << endl;
#endif
  };

  //============================================================
  // Assignment operators
  //============================================================

  // Assigment is SHALLOW, to be consistent with copy.

  Engine_t &operator=(const Engine_t &model)
  {
    if (this == &model) return *this;

    tol_m = model.tol_m;
    size_m = model.size_m;
    for (int d = 0; d < Dim; ++d) {
      firsts_m[d] = model.firsts_m[d];
    };
    domain_m = model.domain_m;
    data_m.initialize(model.data_m);
    for (int i=domain_m[0].first(); i<=domain_m[0].last(); i++) 
      for (int j=domain_m[1].first(); j<=domain_m[1].last(); j++) {
         data_m(i,j).initialize(model.data_m(i,j));
//	 data_m(i,j)=model.data_m(i,j);
      };

    return *this;
  }

  //============================================================
  // Accessor and Mutator functions:
  //============================================================

  // Element access via Loc.

  Element_t read(const Loc<Dim> &loc) const
  {
    Loc<2> L01(loc[0], loc[1]), 
           L23(loc[2], loc[3]);
    return data_m(L01)(L23);
  }

  ElementRef_t operator()(const Loc<Dim> &loc) const
  {
    Loc<2> L01(loc[0], loc[1]), 
           L23(loc[2], loc[3]);
    return data_m(L01)(L23);
  }

  // Element access via ints for speed.

  Element_t read(int i1, int i2, int i3, int i4) const
  {
    return data_m(i1,i2)(i3,i4);
  }

  ElementRef_t operator()(int i1, int i2, int i3, int i4) const
  {
    return data_m(i1,i2)(i3,i4);
  }

  // Get a private copy of data viewed by this Engine.

  Engine_t &makeOwnCopy()
  {

    data_m.makeOwnCopy();

    return *this;
  }

  //---------------------------------------------------------
  // Return/set the domain.

  inline const Domain_t &domain() const
  {
    return domain_m; 
  }

  //---------------------------------------------------------
  // Return the first index value for the specified direction.

  inline int first(int i) const
  {
    PAssert(i >= 0 && i < Dim);
    return firsts_m[i];
  }

  //---------------------------------------------------------
  // Return the size

  inline int size() const
  {
    return size_m;
  }

  //---------------------------------------------------------
  // Return number of free cells

  inline int free() const
  {
    int f=0;

    for (int i=domain_m[0].first(); i<=domain_m[0].last(); i++) 
      for (int j=domain_m[1].first(); j<=domain_m[1].last(); j++) {
	 f += data_m(i,j).engine().free();
      };

    return f;
  }

  //---------------------------------------------------------
  // Return tolerance:
  // (see also realization of pack() !!!)

  inline T& tolerance()
  {
    return tol_m;
  }

  //---------------------------------------------------------
  // Packing functions for the sparse array:

  void pack(void) const
  {
    for (int i=domain_m[0].first(); i<=domain_m[0].last(); i++) 
      for (int j=domain_m[1].first(); j<=domain_m[1].last(); j++) {
	 data_m(i,j).engine().tolerance()=tol_m;
	 data_m(i,j).engine().pack();
      };
  }

  //---------------------------------------------------------
  // Print sparse array:

  void print(void) const
  {
//    cout << endl << " Non-zero elements of 4D Sparse Array:" << endl;
    for (int i=domain_m[0].first(); i<=domain_m[0].last(); i++) 
      for (int j=domain_m[1].first(); j<=domain_m[1].last(); j++) {
         cout << "(" << i << "," << j << "):" << endl;
	 data_m(i,j).engine().print();
      };
  }

  //---------------------------------------------------------
  // Write sparse array:
  template<class Out>
  void write(Out& os) const
  {
    // Write the extent vector: this is separated by 'x's, e.g.
    // (1, 10) x (-4, 4) x (-5, 5)

    for (int i=0; i < Dim; ++i) {
      os << "(";
      os << domain_m[i].first(); 
      os << ", ";
      os << domain_m[i].last(); 
      os << ")";
      if (i != Dim-1) os << " x ";
    }
    os << endl << "[ ";

    for (int i=domain_m[0].first(); i<=domain_m[0].last(); i++) 
      for (int j=domain_m[1].first(); j<=domain_m[1].last(); j++) {
	 data_m(i,j).engine().write(os);
      };

    os << "]" << endl;
  }

private:

  //==========================================================
  // Private data
  //==========================================================

  T tol_m;
  int size_m;
  int firsts_m[Dim];
  Domain_t domain_m;
  Array<2, Array<2,T,Sparse> > data_m; 

};

#undef Dim

#endif // SPL_ENGINE_SPARSEENGINE4_H