parallel_solver.h

Go to the documentation of this file.

#ifndef PARALLEL_SOLVER_H
#define PARALLEL_SOLVER_H

#include "parallel_array.h"
#include "parallel_batch.h"
#include "parallel_math.h"
#include "parallel_reduce.h"
#include "parallel_timer.h"

#include "util.h"

namespace parallel_ode
{

using ::parallel_utils::ParallelHDArray;
using ::parallel_utils::MemManager;
using ::parallel_utils::CopyType;


namespace ParallelFlags
{
enum ParallelFlag { PARALLEL_NONE       = 0x000,
                    PARALLEL_PREPROCESS = 0x001,
                    PARALLEL_ASYNC      = 0x002,
                    PARALLEL_RANDOMIZE  = 0x004,
                    PARALLEL_ATOMICS    = 0x008,
                    PARALLEL_ALIGN      = 0x010,
                    PARALLEL_PINNED     = 0x020,
                    PARALLEL_WC         = 0x040,
                    PARALLEL_REDUCE     = 0x080,
                    PARALLEL_COMPACT    = 0x100};
}
typedef ParallelFlags::ParallelFlag ParallelFlag;

const int DEFAULT_FLAGS =
    ParallelFlags::PARALLEL_ALIGN
    | ParallelFlags::PARALLEL_PINNED
    | ParallelFlags::PARALLEL_REDUCE
    | ParallelFlags::PARALLEL_ASYNC
    | ParallelFlags::PARALLEL_RANDOMIZE
    ;//| PARALLEL_ATOMICS;

/*  ParallelPGSSolver is the abstract base class encapsulating functionality common to each of the parallel
 *  Quickstep solvers.  It takes the data from ODE, batches it appropriately, and relies on derived classes
 *  to actually drive the solver.  From there the data is transferred back to ODE.
 */
template<typename CudaT, typename ParamsT, ParallelType PType>
class ParallelPGSSolver
{
public:
  struct SolverParams;
  typedef typename vec4<CudaT>::Type Vec4T;
  typedef const CudaT* CudaTPtr;
  typedef CudaT* CudaTMutablePtr;
  typedef const ParamsT* ParamsTPtr;
  typedef ParamsT* ParamsTMutablePtr;
  typedef MemManager<CudaT,PType> PMemManager;
  typedef typename PMemManager::mem_flags MemFlags;

  ParallelPGSSolver( int parallelFlags = DEFAULT_FLAGS,
                     BatchType batchType = BatchTypes::DEFAULT_BATCH_TYPE,
                     ReduceType reduceType = ReduceTypes::DEFAULT_REDUCE_TYPE,
                     uint numBatches = ParallelOptions::MAXBATCHES )
      : parallelFlags_( parallelFlags ),
        parallelParams_( NULL ),
        batchStrategyType_( batchType ),
        batchStrategy_( NULL ),
        reduceStrategyType_( reduceType),
        reduceStrategy_( NULL ),
        batchRepetitionCount_( numBatches, 0 ),
        batchIndices_( numBatches, -1 ),
        batchSizes_( numBatches, 0 ),
        numBatches_( numBatches ),
        m_(0),
        constraintStride_(0),
        n_(0),
        bodyStride_(0),
        reduceStride_(0),
        bInit_(false) {
  }

  virtual ~ParallelPGSSolver( ) {
    delete batchStrategy_;
    delete reduceStrategy_;
  }

  virtual void worldSolve( SolverParams* params );

  virtual void initialize( );

protected:


  void syncODEToDevice( );
  void syncDeviceToODE( );
  void loadBatches( int* jb );

  void setConstraintSize( int m, bool bResizeArrays = true );
  void setBodySize( int n, bool bResizeArrays = true );


  virtual void preProcessHost( );
  virtual void preProcessDevice( const CudaT sorParam, const CudaT stepsize );
  virtual void solveAndReduce( const int offset, const int batchSize ) = 0;

  virtual void loadBodies( );
  virtual void loadConstraints( );
  virtual void loadConstants( );
  virtual void loadKernels( );
  virtual void loadSolution( );

  virtual void setMemFlags( MemFlags flags );


  inline void checkInit( ) { if(!bInit_) this->initialize(); }
  inline void setBodyStride( int bodyStride ) { bodyStride_ = bodyStride; }
  inline void setConstraintStride( int constraintStride ) { constraintStride_ = constraintStride; }
  inline void setReduceStride( int reduceStride ) { reduceStride_ = reduceStride; }

  inline int getNumBatches( ) const { return (int)batchSizes_.size(); }
  inline int getNumBodies( ) const { return n_; }
  inline int getBodyStride( ) const { return bodyStride_; }
  inline int getNumConstraints( ) const { return m_; }
  inline int getConstraintStride( ) const { return constraintStride_; }
  inline int getReduceStride( ) const { return reduceStride_; }
  inline int getDefaultAlign( ) const { return ParallelOptions::DEFAULTALIGN; }
  inline CopyType getCopyType( ) const { return asyncEnabled( ) ? parallel_utils::CopyTypes::COPY_ASYNC : parallel_utils::CopyTypes::COPY_SYNC; }

  inline int bID(const int index, const int offset) const { return index + offset * getBodyStride(); }
  inline int cID(const int index, const int offset) const { return index + offset * getConstraintStride(); }
  inline int rID(const int index, const int offset) const { return index + offset * getReduceStride(); }

  inline bool preprocessEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_PREPROCESS ); }
  inline bool asyncEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_ASYNC ); }
  inline bool randomizeEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_RANDOMIZE ); }
  inline bool atomicsEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_ATOMICS ); }
  inline bool alignEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_ALIGN ); }
  inline bool pinnedMemEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_PINNED ); }
  inline bool wcMemEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_WC ); }
  inline bool reduceEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_REDUCE ); }
  inline bool compactScalarsEnabled( ) const { return isEnabled( ParallelFlags::PARALLEL_COMPACT ); }

  inline bool isEnabled( const ParallelFlag flag ) const { return parallelFlags_ & flag; }

  void printConfig( );


  int parallelFlags_;                   
  SolverParams* parallelParams_;        
  BatchType batchStrategyType_;         
  BatchStrategy *batchStrategy_;        
  ReduceType reduceStrategyType_;       
  ReduceStrategy *reduceStrategy_;      
  IntVector constraintIndices_;         
  IntVector reductionOffsets0_;         
  IntVector reductionOffsets1_;         
  IntVector batchRepetitionCount_;      
  IntVector batchIndices_;              
  IntVector batchSizes_;                


  int numBatches_;                      
  int m_;                               
  int constraintStride_;                
  int n_;                               
  int bodyStride_;                      
  int reduceStride_;                    
  bool bInit_;                          
  ParallelHDArray<Vec4T,PType> j0;      
  ParallelHDArray<Vec4T,PType> ij0;     
  ParallelHDArray<CudaT,PType> lambda0;     
  ParallelHDArray<CudaT,PType> adcfm;       
  ParallelHDArray<CudaT,PType> rhs;         
  ParallelHDArray<CudaT,PType> lohiD;       
  ParallelHDArray<int4,PType> bodyIDs;  
  ParallelHDArray<int,PType> fIDs;      
  ParallelHDArray<CudaT,PType> iMass;       
  ParallelHDArray<Vec4T,PType> i0;      
  ParallelHDArray<Vec4T,PType> bodyFAcc; 
  ParallelHDArray<Vec4T,PType> bodyTAcc; 
  ParallelHDArray<Vec4T,PType> bodyFAccReduction; 
  ParallelHDArray<Vec4T,PType> bodyTAccReduction; 
public:

  struct SolverParams {
    SolverParams( dxWorldProcessContext *contextIn, const dxQuickStepParameters *qsIn,
                  const int mIn, const int nbIn, ParamsTMutablePtr JIn, int *jbIn, dxBody * const *bodyIn,
                  ParamsTPtr invIIn, ParamsTMutablePtr lambdaIn, ParamsTMutablePtr fcIn, ParamsTMutablePtr bIn,
                  ParamsTPtr loIn, ParamsTPtr hiIn, ParamsTPtr cfmIn, ParamsTMutablePtr iMJIn, ParamsTMutablePtr AdcfmIn,
                  const int *findexIn, const ParamsT stepsizeIn )
    : context(contextIn), qs(qsIn), m(mIn),  nb(nbIn),  J(JIn), jb(jbIn), body(bodyIn),
      invI(invIIn),  lambda(lambdaIn),  fc(fcIn),  b(bIn),  lo(loIn),  hi(hiIn),  cfm(cfmIn),
      iMJ(iMJIn),  Adcfm(AdcfmIn), findex(findexIn),  stepsize(stepsizeIn) { }

    dxWorldProcessContext *context;
    const dxQuickStepParameters *qs;
    const int m;
    const int nb;
    ParamsTMutablePtr J;
    int *jb;
    dxBody * const *body;
    ParamsTPtr invI;
    ParamsTMutablePtr lambda;
    ParamsTMutablePtr fc;
    ParamsTMutablePtr b;
    ParamsTPtr lo;
    ParamsTPtr hi;
    ParamsTPtr cfm;
    ParamsTMutablePtr iMJ;
    ParamsTMutablePtr Adcfm;
    const int *findex;
    const ParamsT stepsize;
  };

};

}

#endif