text: #ifndef __LULESH
  #define __LULESH
  
  #include <CL/sycl.hpp>
  #include <dpct/dpct.hpp>
  
  // Emulate the number of CPU threads
  #define NT 2
  
  //**************************************************
  // Allow flexibility for arithmetic representations 
  //**************************************************
  
  #define MAX(a, b) ( ((a) > (b)) ? (a) : (b))
  
  // Precision specification
  typedef float        real4 ;
  typedef double       real8 ;
  typedef long double  real10 ;  // 10 bytes on x86
  
  typedef int    Index_t ; // array subscript and loop index
  typedef real8  Real_t ;  // floating point representation
  typedef int    Int_t ;   // integer representation
  
  enum { VolumeError = -1, QStopError = -2 } ;
  
  inline real4 SQRT(real4 arg) { return sqrtf(arg); }
  inline real8 SQRT(real8 arg) { return sqrt(arg); }
  inline real10 SQRT(real10 arg) { return sqrtl(arg) ; }
  
  inline real4 CBRT(real4 arg) { return cbrtf(arg); }
  inline real8 CBRT(real8 arg) { return cbrt(arg); }
  inline real10 CBRT(real10 arg) { return cbrtl(arg) ; }
  
  inline real4 FABS(real4 arg) { return fabsf(arg); }
  inline real8 FABS(real8 arg) { return fabs(arg); }
  inline real10 FABS(real10 arg) { return fabsl(arg) ; }
  
  
  // Stuff needed for boundary conditions
  // 2 BCs on each of 6 hexahedral faces (12 bits)
  #define XI_M        0x00007
  #define XI_M_SYMM   0x00001
  #define XI_M_FREE   0x00002
  #define XI_M_COMM   0x00004
  
  #define XI_P        0x00038
  #define XI_P_SYMM   0x00008
  #define XI_P_FREE   0x00010
  #define XI_P_COMM   0x00020
  
  #define ETA_M       0x001c0
  #define ETA_M_SYMM  0x00040
  #define ETA_M_FREE  0x00080
  #define ETA_M_COMM  0x00100
  
  #define ETA_P       0x00e00
  #define ETA_P_SYMM  0x00200
  #define ETA_P_FREE  0x00400
  #define ETA_P_COMM  0x00800
  
  #define ZETA_M      0x07000
  #define ZETA_M_SYMM 0x01000
  #define ZETA_M_FREE 0x02000
  #define ZETA_M_COMM 0x04000
  
  #define ZETA_P      0x38000
  #define ZETA_P_SYMM 0x08000
  #define ZETA_P_FREE 0x10000
  #define ZETA_P_COMM 0x20000
  
  // Assume 128 byte coherence
  // Assume Real_t is an "integral power of 2" bytes wide
  #define CACHE_COHERENCE_PAD_REAL (128 / sizeof(Real_t))
  
  #define CACHE_ALIGN_REAL(n) \
     (((n) + (CACHE_COHERENCE_PAD_REAL - 1)) & ~(CACHE_COHERENCE_PAD_REAL-1))
  
  //////////////////////////////////////////////////////
  // Primary data structure
  //////////////////////////////////////////////////////
  
  /*
   * The implementation of the data abstraction used for lulesh
   * resides entirely in the Domain class below.  You can change
   * grouping and interleaving of fields here to maximize data layout
   * efficiency for your underlying architecture or compiler.
   *
   * For example, fields can be implemented as STL objects or
   * raw array pointers.  As another example, individual fields
   * m_x, m_y, m_z could be grouped into
   *
   *    struct { Real_t x, y, z ; } *m_coord ;
   *
   * allowing accessor functions such as
   *
   *  "Real_t &x(Index_t idx) { return m_coord[idx].x ; }"
   *  "Real_t &y(Index_t idx) { return m_coord[idx].y ; }"
   *  "Real_t &z(Index_t idx) { return m_coord[idx].z ; }"
   */
  class Domain {
  
     public:
  
     // Constructor
     Domain(Int_t numRanks, Index_t colLoc,
            Index_t rowLoc, Index_t planeLoc,
            Index_t nx, Int_t tp, Int_t nr, Int_t balance, Int_t cost);
  
     //
     // ALLOCATION
     //
  
     void AllocateNodePersistent(Int_t numNode) // Node-centered
     {
        m_x.resize(numNode);  // coordinates
        m_y.resize(numNode);
        m_z.resize(numNode);
  
        m_xd.resize(numNode); // velocities
        m_yd.resize(numNode);
        m_zd.resize(numNode);
  
        m_xdd.resize(numNode); // accelerations
        m_ydd.resize(numNode);
        m_zdd.resize(numNode);
  
        m_fx.resize(numNode);  // forces
        m_fy.resize(numNode);
        m_fz.resize(numNode);
  
        m_nodalMass.resize(numNode);  // mass
     }
  
     void AllocateElemPersistent(Int_t numElem) // Elem-centered
     {
        m_nodelist.resize(8*numElem);
  
        // elem connectivities through face
        m_lxim.resize(numElem);
        m_lxip.resize(numElem);
        m_letam.resize(numElem);
        m_letap.resize(numElem);
        m_lzetam.resize(numElem);
        m_lzetap.resize(numElem);
  
        m_elemBC.resize(numElem);
  
        m_e.resize(numElem);
        m_p.resize(numElem);
  
        m_q.resize(numElem);
        m_ql.resize(numElem);
        m_qq.resize(numElem);
  
        m_v.resize(numElem);
  
        m_volo.resize(numElem);
        m_delv.resize(numElem);
        m_vdov.resize(numElem);
  
        m_arealg.resize(numElem);
  
        m_ss.resize(numElem);
  
        m_elemMass.resize(numElem);
  
        m_elemRep.resize(numElem);
  
        m_elemElem.resize(numElem);
     }
  
     void AllocateGradients(Int_t numElem, Int_t allElem)
     {
        // Position gradients
        m_delx_xi.resize(numElem) ;
        m_delx_eta.resize(numElem) ;
        m_delx_zeta.resize(numElem) ;
  
        // Velocity gradients
        m_delv_xi.resize(allElem) ;
        m_delv_eta.resize(allElem);
        m_delv_zeta.resize(allElem) ;
     }
  
     void DeallocateGradients()
     {
        m_delx_zeta.clear() ;
        m_delx_eta.clear() ;
        m_delx_xi.clear() ;
  
        m_delv_zeta.clear() ;
        m_delv_eta.clear() ;
        m_delv_xi.clear() ;
     }
  
     void AllocateStrains(Int_t numElem)
     {
        m_dxx.resize(numElem) ;
        m_dyy.resize(numElem) ;
        m_dzz.resize(numElem) ;
     }
  
     void DeallocateStrains()
     {
        m_dzz.clear() ;
        m_dyy.clear() ;
        m_dxx.clear() ;
     }
     
     //
     // ACCESSORS
     //
  
     // Node-centered
  
     // Nodal coordinates
     Real_t& x(Index_t idx)    { return m_x[idx] ; }
     Real_t& y(Index_t idx)    { return m_y[idx] ; }
     Real_t& z(Index_t idx)    { return m_z[idx] ; }
  
     // Nodal velocities
     Real_t& xd(Index_t idx)   { return m_xd[idx] ; }
     Real_t& yd(Index_t idx)   { return m_yd[idx] ; }
     Real_t& zd(Index_t idx)   { return m_zd[idx] ; }
  
     // Nodal accelerations
     Real_t& xdd(Index_t idx)  { return m_xdd[idx] ; }
     Real_t& ydd(Index_t idx)  { return m_ydd[idx] ; }
     Real_t& zdd(Index_t idx)  { return m_zdd[idx] ; }
  
     // Nodal forces
     Real_t& fx(Index_t idx)   { return m_fx[idx] ; }
     Real_t& fy(Index_t idx)   { return m_fy[idx] ; }
     Real_t& fz(Index_t idx)   { return m_fz[idx] ; }
  
     // Nodal mass
     Real_t& nodalMass(Index_t idx) { return m_nodalMass[idx] ; }
  
     // Nodes on symmertry planes
     Index_t symmX(Index_t idx) { return m_symmX[idx] ; }
     Index_t symmY(Index_t idx) { return m_symmY[idx] ; }
     Index_t symmZ(Index_t idx) { return m_symmZ[idx] ; }
     bool symmXempty()          { return m_symmX.empty(); }
     bool symmYempty()          { return m_symmY.empty(); }
     bool symmZempty()          { return m_symmZ.empty(); }
  
     //
     // Element-centered
     //
     Index_t&  regElemSize(Index_t idx) { return m_regElemSize[idx] ; }
     Index_t&  regNumList(Index_t idx) { return m_regNumList[idx] ; }
     Index_t*  regNumList()            { return &m_regNumList[0] ; }
     Index_t*  regElemlist(Int_t r)    { return m_regElemlist[r] ; }
     Index_t&  regElemlist(Int_t r, Index_t idx) { return m_regElemlist[r][idx] ; }
  
     Index_t*  nodelist(Index_t idx)    { return &m_nodelist[Index_t(8)*idx] ; }
  
     // elem connectivities through face
     Index_t&  lxim(Index_t idx) { return m_lxim[idx] ; }
     Index_t&  lxip(Index_t idx) { return m_lxip[idx] ; }
     Index_t&  letam(Index_t idx) { return m_letam[idx] ; }
     Index_t&  letap(Index_t idx) { return m_letap[idx] ; }
     Index_t&  lzetam(Index_t idx) { return m_lzetam[idx] ; }
     Index_t&  lzetap(Index_t idx) { return m_lzetap[idx] ; }
  
     // elem face symm/free-surface flag
     Int_t&  elemBC(Index_t idx) { return m_elemBC[idx] ; }
  
     // Principal strains - temporary
     Real_t& dxx(Index_t idx)  { return m_dxx[idx] ; }
     Real_t& dyy(Index_t idx)  { return m_dyy[idx] ; }
     Real_t& dzz(Index_t idx)  { return m_dzz[idx] ; }
  
     // Velocity gradient - temporary
     Real_t& delv_xi(Index_t idx)    { return m_delv_xi[idx] ; }
     Real_t& delv_eta(Index_t idx)   { return m_delv_eta[idx] ; }
     Real_t& delv_zeta(Index_t idx)  { return m_delv_zeta[idx] ; }
  
     // Position gradient - temporary
     Real_t& delx_xi(Index_t idx)    { return m_delx_xi[idx] ; }
     Real_t& delx_eta(Index_t idx)   { return m_delx_eta[idx] ; }
     Real_t& delx_zeta(Index_t idx)  { return m_delx_zeta[idx] ; }
  
     // Energy
     Real_t& e(Index_t idx)          { return m_e[idx] ; }
  
     // Pressure
     Real_t& p(Index_t idx)          { return m_p[idx] ; }
  
     // Artificial viscosity
     Real_t& q(Index_t idx)          { return m_q[idx] ; }
  
     // Linear term for q
     Real_t& ql(Index_t idx)         { return m_ql[idx] ; }
     // Quadratic term for q
     Real_t& qq(Index_t idx)         { return m_qq[idx] ; }
  
     // Relative volume
     Real_t& v(Index_t idx)          { return m_v[idx] ; }
     Real_t& delv(Index_t idx)       { return m_delv[idx] ; }
  
     // Reference volume
     Real_t& volo(Index_t idx)       { return m_volo[idx] ; }
  
     // volume derivative over volume
     Real_t& vdov(Index_t idx)       { return m_vdov[idx] ; }
  
     // Element characteristic length
     Real_t& arealg(Index_t idx)     { return m_arealg[idx] ; }
  
     // Sound speed
     Real_t& ss(Index_t idx)         { return m_ss[idx] ; }
  
     // Element mass
     Real_t& elemMass(Index_t idx)  { return m_elemMass[idx] ; }
  
     // Element mass
     Index_t& elemRep(Index_t idx)  { return m_elemRep[idx] ; }
  
     // Element mass
     Index_t& elemElem(Index_t idx)  { return m_elemElem[idx] ; }
  
     Index_t nodeElemCount(Index_t idx)
     { return m_nodeElemStart[idx+1] - m_nodeElemStart[idx] ; }
  
     Index_t *nodeElemCornerList(Index_t idx)
     { return &m_nodeElemCornerList[m_nodeElemStart[idx]] ; }
  
     // Parameters 
  
     // Cutoffs
     Real_t u_cut() const               { return m_u_cut ; }
     Real_t e_cut() const               { return m_e_cut ; }
     Real_t p_cut() const               { return m_p_cut ; }
     Real_t q_cut() const               { return m_q_cut ; }
     Real_t v_cut() const               { return m_v_cut ; }
  
     // Other constants (usually are settable via input file in real codes)
     Real_t hgcoef() const              { return m_hgcoef ; }
     Real_t qstop() const               { return m_qstop ; }
     Real_t monoq_max_slope() const     { return m_monoq_max_slope ; }
     Real_t monoq_limiter_mult() const  { return m_monoq_limiter_mult ; }
     Real_t ss4o3() const               { return m_ss4o3 ; }
     Real_t qlc_monoq() const           { return m_qlc_monoq ; }
     Real_t qqc_monoq() const           { return m_qqc_monoq ; }
     Real_t qqc() const                 { return m_qqc ; }
  
     Real_t eosvmax() const             { return m_eosvmax ; }
     Real_t eosvmin() const             { return m_eosvmin ; }
     Real_t pmin() const                { return m_pmin ; }
     Real_t emin() const                { return m_emin ; }
     Real_t dvovmax() const             { return m_dvovmax ; }
     Real_t refdens() const             { return m_refdens ; }
  
     // Timestep controls, etc...
     Real_t& time()                 { return m_time ; }
     Real_t& deltatime()            { return m_deltatime ; }
     Real_t& deltatimemultlb()      { return m_deltatimemultlb ; }
     Real_t& deltatimemultub()      { return m_deltatimemultub ; }
     Real_t& stoptime()             { return m_stoptime ; }
     Real_t& dtcourant()            { return m_dtcourant ; }
     Real_t& dthydro()              { return m_dthydro ; }
     Real_t& dtmax()                { return m_dtmax ; }
     Real_t& dtfixed()              { return m_dtfixed ; }
  
     Int_t&  cycle()                { return m_cycle ; }
     Index_t&  numRanks()           { return m_numRanks ; }
  
     Index_t&  colLoc()             { return m_colLoc ; }
     Index_t&  rowLoc()             { return m_rowLoc ; }
     Index_t&  planeLoc()           { return m_planeLoc ; }
     Index_t&  tp()                 { return m_tp ; }
  
     Index_t&  sizeX()              { return m_sizeX ; }
     Index_t&  sizeY()              { return m_sizeY ; }
     Index_t&  sizeZ()              { return m_sizeZ ; }
     Index_t&  numReg()             { return m_numReg ; }
     Int_t&  cost()             { return m_cost ; }
     Index_t&  numElem()            { return m_numElem ; }
     Index_t&  numNode()            { return m_numNode ; }
     
     Index_t&  maxPlaneSize()       { return m_maxPlaneSize ; }
     Index_t&  maxEdgeSize()        { return m_maxEdgeSize ; }
     
    // private:
  
     void BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems);
     void SetupThreadSupportStructures();
     void CreateRegionIndexSets(Int_t nreg, Int_t balance);
     void SetupCommBuffers(Int_t edgeNodes);
     void SetupSymmetryPlanes(Int_t edgeNodes);
     void SetupElementConnectivities(Int_t edgeElems);
     void SetupBoundaryConditions(Int_t edgeElems);
  
     //
     // IMPLEMENTATION
     //
  
     /* Node-centered */
     std::vector<Real_t> m_x ;  /* coordinates */
     std::vector<Real_t> m_y ;
     std::vector<Real_t> m_z ;
  
     std::vector<Real_t> m_xd ; /* velocities */
     std::vector<Real_t> m_yd ;
     std::vector<Real_t> m_zd ;
  
     std::vector<Real_t> m_xdd ; /* accelerations */
     std::vector<Real_t> m_ydd ;
     std::vector<Real_t> m_zdd ;
  
     std::vector<Real_t> m_fx ;  /* forces */
     std::vector<Real_t> m_fy ;
     std::vector<Real_t> m_fz ;
  
     std::vector<Real_t> m_nodalMass ;  /* mass */
  
     std::vector<Index_t> m_symmX ;  /* symmetry plane nodesets */
     std::vector<Index_t> m_symmY ;
     std::vector<Index_t> m_symmZ ;
  
     // Element-centered
  
     // Region information
     Int_t    m_numReg ;
     Int_t    m_cost; //imbalance cost
     Index_t *m_regElemSize ;   // Size of region sets
     Index_t *m_regNumList ;    // Region number per domain element
     Index_t **m_regElemlist ;  // region indexset 
  
     std::vector<Index_t>  m_nodelist ;     /* elemToNode connectivity */
  
     std::vector<Index_t>  m_lxim ;  /* element connectivity across each face */
     std::vector<Index_t>  m_lxip ;
     std::vector<Index_t>  m_letam ;
     std::vector<Index_t>  m_letap ;
     std::vector<Index_t>  m_lzetam ;
     std::vector<Index_t>  m_lzetap ;
  
     std::vector<Int_t>    m_elemBC ;  /* symmetry/free-surface flags for each elem face */
  
     std::vector<Real_t> m_dxx ;  /* principal strains -- temporary */
     std::vector<Real_t> m_dyy ;
     std::vector<Real_t> m_dzz ;
  
     std::vector<Real_t> m_delv_xi ;    /* velocity gradient -- temporary */
     std::vector<Real_t> m_delv_eta ;
     std::vector<Real_t> m_delv_zeta ;
  
     std::vector<Real_t> m_delx_xi ;    /* coordinate gradient -- temporary */
     std::vector<Real_t> m_delx_eta ;
     std::vector<Real_t> m_delx_zeta ;
     
     std::vector<Real_t> m_e ;   /* energy */
  
     std::vector<Real_t> m_p ;   /* pressure */
     std::vector<Real_t> m_q ;   /* q */
     std::vector<Real_t> m_ql ;  /* linear term for q */
     std::vector<Real_t> m_qq ;  /* quadratic term for q */
  
     std::vector<Real_t> m_v ;     /* relative volume */
     std::vector<Real_t> m_volo ;  /* reference volume */
     std::vector<Real_t> m_vnew ;  /* new relative volume -- temporary */
     std::vector<Real_t> m_delv ;  /* m_vnew - m_v */
     std::vector<Real_t> m_vdov ;  /* volume derivative over volume */
  
     std::vector<Real_t> m_arealg ;  /* characteristic length of an element */
     
     std::vector<Real_t> m_ss ;      /* "sound speed" */
  
     std::vector<Real_t> m_elemMass ;  /* mass */
  
     std::vector<Index_t> m_elemRep ;  /* reps */
  
     std::vector<Index_t> m_elemElem ;  /* elems */
  
     // Cutoffs (treat as constants)
     const Real_t  m_e_cut ;             // energy tolerance 
     const Real_t  m_p_cut ;             // pressure tolerance 
     const Real_t  m_q_cut ;             // q tolerance 
     const Real_t  m_v_cut ;             // relative volume tolerance 
     const Real_t  m_u_cut ;             // velocity tolerance 
  
     // Other constants (usually setable, but hardcoded in this proxy app)
  
     const Real_t  m_hgcoef ;            // hourglass control 
     const Real_t  m_ss4o3 ;
     const Real_t  m_qstop ;             // excessive q indicator 
     const Real_t  m_monoq_max_slope ;
     const Real_t  m_monoq_limiter_mult ;
     const Real_t  m_qlc_monoq ;         // linear term coef for q 
     const Real_t  m_qqc_monoq ;         // quadratic term coef for q 
     const Real_t  m_qqc ;
     const Real_t  m_eosvmax ;
     const Real_t  m_eosvmin ;
     const Real_t  m_pmin ;              // pressure floor 
     const Real_t  m_emin ;              // energy floor 
     const Real_t  m_dvovmax ;           // maximum allowable volume change 
     const Real_t  m_refdens ;           // reference density 
  
     // Variables to keep track of timestep, simulation time, and cycle
     Real_t  m_dtcourant ;         // courant constraint 
     Real_t  m_dthydro ;           // volume change constraint 
     Int_t   m_cycle ;             // iteration count for simulation 
     Real_t  m_dtfixed ;           // fixed time increment 
     Real_t  m_time ;              // current time 
     Real_t  m_deltatime ;         // variable time increment 
     Real_t  m_deltatimemultlb ;
     Real_t  m_deltatimemultub ;
     Real_t  m_dtmax ;             // maximum allowable time increment 
     Real_t  m_stoptime ;          // end time for simulation 
  
  
     Int_t   m_numRanks ;
  
     Index_t m_colLoc ;
     Index_t m_rowLoc ;
     Index_t m_planeLoc ;
     Index_t m_tp ;
  
     Index_t m_sizeX ;
     Index_t m_sizeY ;
     Index_t m_sizeZ ;
     Index_t m_numElem ;
     Index_t m_numNode ;
  
     Index_t m_maxPlaneSize ;
     Index_t m_maxEdgeSize ;
  
     // OMP hack 
     Index_t *m_nodeElemStart ;
     Index_t *m_nodeElemCornerList ;
  
     // Used in setup
     Index_t m_rowMin, m_rowMax;
     Index_t m_colMin, m_colMax;
     Index_t m_planeMin, m_planeMax ;
  } ;
  
  typedef Real_t &(Domain::* Domain_member )(Index_t) ;
  
  struct cmdLineOpts {
     Int_t its; // -i 
     Int_t nx;  // -s 
     Int_t numReg; // -r 
     Int_t numFiles; // -f
     Int_t showProg; // -p
     Int_t quiet; // -q
     Int_t viz; // -v 
     Int_t cost; // -c
     Int_t balance; // -b
     Int_t iteration_cap; // -z
  };
  
  
  
  // Function Prototypes
  
  Real_t CalcElemVolume( const Real_t x[8], const Real_t y[8], const Real_t z[8] );
  
  
  // lulesh-util
  void ParseCommandLineOptions(int argc, char *argv[],
                               Int_t myRank, struct cmdLineOpts *opts);
  void VerifyAndWriteFinalOutput(Real_t elapsed_time,
                                 Domain& locDom,
                                 Int_t nx,
                                 Int_t numRanks);
  
  // lulesh-viz
  void DumpToVisit(Domain& domain, int numFiles, int myRank, int numRanks);
  
  // lulesh-init
  void InitMeshDecomp(Int_t numRanks, Int_t myRank,
                      Int_t *col, Int_t *row, Int_t *plane, Int_t *side);
  
  #endif
  