Continuing to implement multi-GPU support

platforms/opencl/include/OpenCLPlatform.h

@@ -55,7 +55,7 @@ public:
     void contextCreated(ContextImpl& context, const std::map<std::string, std::string>& properties) const;
     void contextDestroyed(ContextImpl& context) const;
     /**
-     * This is the name of the parameter for selecting which OpenCL device to use.
+     * This is the name of the parameter for selecting which OpenCL device or devices to use.
      */
     static const std::string& OpenCLDeviceIndex() {
         static const std::string key = "OpenCLDeviceIndex";
@@ -65,7 +65,7 @@ public:
 
 class OpenCLPlatform::PlatformData {
 public:
-    PlatformData(int numParticles, int deviceIndex);
+    PlatformData(int numParticles, const std::string& deviceIndexProperty);
     ~PlatformData();
     void initializeContexts(const System& system);
     std::vector<OpenCLContext*> contexts;

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -56,6 +56,10 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
             return new OpenCLParallelCalcCustomTorsionForceKernel(name, platform, data, context.getSystem());
         if (name == CalcNonbondedForceKernel::Name())
             return new OpenCLParallelCalcNonbondedForceKernel(name, platform, data, context.getSystem());
+        if (name == CalcCustomExternalForceKernel::Name())
+            return new OpenCLParallelCalcCustomExternalForceKernel(name, platform, data, context.getSystem());
+        if (name == CalcCustomHbondForceKernel::Name())
+            return new OpenCLParallelCalcCustomHbondForceKernel(name, platform, data, context.getSystem());
     }
     OpenCLContext& cl = *data.contexts[0];
     if (name == CalcForcesAndEnergyKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -1685,6 +1685,8 @@ OpenCLCalcGBSAOBCForceKernel::~OpenCLCalcGBSAOBCForceKernel() {
 }
 
 void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
+    if (cl.getPlatformData().contexts.size() > 1)
+        throw OpenMMException("GBSAOBCForce does not support using multiple OpenCL devices");
     OpenCLNonbondedUtilities& nb = cl.getNonbondedUtilities();
     params = new OpenCLArray<mm_float2>(cl, cl.getPaddedNumAtoms(), "gbsaObcParams");
     bornRadii = new OpenCLArray<cl_float>(cl, cl.getPaddedNumAtoms(), "bornRadii");
@@ -1868,6 +1870,8 @@ OpenCLCalcCustomGBForceKernel::~OpenCLCalcCustomGBForceKernel() {
 }
 
 void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
+    if (cl.getPlatformData().contexts.size() > 1)
+        throw OpenMMException("CustomGBForce does not support using multiple OpenCL devices");
     bool useExclusionsForValue = false;
     vector<string> computedValueNames(force.getNumComputedValues());
     vector<string> computedValueExpressions(force.getNumComputedValues());
@@ -2660,7 +2664,12 @@ OpenCLCalcCustomExternalForceKernel::~OpenCLCalcCustomExternalForceKernel() {
 }
 
 void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
-    numParticles = force.getNumParticles();
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumParticles()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumParticles()/numContexts;
+    numParticles = endIndex-startIndex;
+    if (numParticles == 0)
+        return;
     params = new OpenCLParameterSet(cl, force.getNumPerParticleParameters(), numParticles, "customExternalParams");
     indices = new OpenCLArray<cl_int>(cl, numParticles, "customExternalIndices");
     string extraArguments;
@@ -2672,7 +2681,7 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
     vector<cl_int> indicesVector(numParticles);
     for (int i = 0; i < numParticles; i++) {
         vector<double> parameters;
-        force.getParticleParameters(i, indicesVector[i], parameters);
+        force.getParticleParameters(startIndex+i, indicesVector[i], parameters);
         paramVector[i].resize(parameters.size());
         for (int j = 0; j < (int) parameters.size(); j++)
             paramVector[i][j] = (cl_float) parameters[j];
@@ -2732,6 +2741,8 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
 }
 
 double OpenCLCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (numParticles == 0)
+        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -2879,8 +2890,13 @@ static void applyDonorAndAcceptorForces(stringstream& applyToDonor, stringstream
 void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
     // Record the lists of donors and acceptors, and the parameters for each one.
 
-    numDonors = force.getNumDonors();
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumDonors()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumDonors()/numContexts;
+    numDonors = endIndex-startIndex;
     numAcceptors = force.getNumAcceptors();
+    if (numDonors == 0 || numAcceptors == 0)
+        return;
     int numParticles = system.getNumParticles();
     donors = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondDonors");
     acceptors = new OpenCLArray<mm_int4>(cl, numAcceptors, "customHbondAcceptors");
@@ -2892,7 +2908,7 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     vector<mm_int4> donorVector(numDonors);
     for (int i = 0; i < numDonors; i++) {
         vector<double> parameters;
-        force.getDonorParameters(i, donorVector[i].x, donorVector[i].y, donorVector[i].z, parameters);
+        force.getDonorParameters(startIndex+i, donorVector[i].x, donorVector[i].y, donorVector[i].z, parameters);
         donorParamVector[i].resize(parameters.size());
         for (int j = 0; j < (int) parameters.size(); j++)
             donorParamVector[i][j] = (cl_float) parameters[j];
@@ -2943,6 +2959,9 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     for (int i = 0; i < force.getNumExclusions(); i++) {
         int donor, acceptor;
         force.getExclusionParticles(i, donor, acceptor);
+        if (donor < startIndex || donor >= endIndex)
+            continue;
+        donor -= startIndex;
         if (donorExclusionVector[donor].x == -1)
             donorExclusionVector[donor].x = acceptor;
         else if (donorExclusionVector[donor].y == -1)
@@ -3172,8 +3191,8 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
     map<string, string> defines;
     defines["PADDED_NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
-    defines["NUM_DONORS"] = intToString(force.getNumDonors());
-    defines["NUM_ACCEPTORS"] = intToString(force.getNumAcceptors());
+    defines["NUM_DONORS"] = intToString(numDonors);
+    defines["NUM_ACCEPTORS"] = intToString(numAcceptors);
     defines["M_PI"] = doubleToString(M_PI);
     if (force.getNonbondedMethod() != CustomHbondForce::NoCutoff) {
         defines["USE_CUTOFF"] = "1";
@@ -3189,6 +3208,8 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
 }
 
 double OpenCLCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (numDonors == 0 || numAcceptors == 0)
+        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -240,3 +240,39 @@ double OpenCLParallelCalcNonbondedForceKernel::execute(ContextImpl& context, boo
         energy += getKernel(i).execute(context, includeForces, includeEnergy);
     return energy;
 }
+
+OpenCLParallelCalcCustomExternalForceKernel::OpenCLParallelCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) :
+        CalcCustomExternalForceKernel(name, platform), data(data) {
+    for (int i = 0; i < (int) data.contexts.size(); i++)
+        kernels.push_back(Kernel(new OpenCLCalcCustomExternalForceKernel(name, platform, *data.contexts[i], system)));
+}
+
+void OpenCLParallelCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).initialize(system, force);
+}
+
+double OpenCLParallelCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    double energy = 0.0;
+    for (int i = 0; i < (int) kernels.size(); i++)
+        energy += getKernel(i).execute(context, includeForces, includeEnergy);
+    return energy;
+}
+
+OpenCLParallelCalcCustomHbondForceKernel::OpenCLParallelCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) :
+        CalcCustomHbondForceKernel(name, platform), data(data) {
+    for (int i = 0; i < (int) data.contexts.size(); i++)
+        kernels.push_back(Kernel(new OpenCLCalcCustomHbondForceKernel(name, platform, *data.contexts[i], system)));
+}
+
+void OpenCLParallelCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).initialize(system, force);
+}
+
+double OpenCLParallelCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    double energy = 0.0;
+    for (int i = 0; i < (int) kernels.size(); i++)
+        energy += getKernel(i).execute(context, includeForces, includeEnergy);
+    return energy;
+}

platforms/opencl/src/OpenCLParallelKernels.h

@@ -346,102 +346,15 @@ private:
     std::vector<Kernel> kernels;
 };
 
-/**
- * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
- */
-class OpenCLParallelCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
-public:
-    OpenCLParallelCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcCustomNonbondedForceKernel(name, platform),
-            data(data) {
-    }
-    ~OpenCLParallelCalcCustomNonbondedForceKernel();
-    /**
-     * Initialize the kernel.
-     *
-     * @param system     the System this kernel will be applied to
-     * @param force      the CustomNonbondedForce this kernel will be used for
-     */
-    void initialize(const System& system, const CustomNonbondedForce& force);
-    /**
-     * Execute the kernel to calculate the forces and/or energy.
-     *
-     * @param context        the context in which to execute this kernel
-     * @param includeForces  true if forces should be calculated
-     * @param includeEnergy  true if the energy should be calculated
-     * @return the potential energy due to the force
-     */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-private:
-    OpenCLPlatform::PlatformData& data;
-};
-
-/**
- * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
- */
-class OpenCLParallelCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
-public:
-    OpenCLParallelCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : CalcGBSAOBCForceKernel(name, platform),
-            data(data) {
-    }
-    ~OpenCLParallelCalcGBSAOBCForceKernel();
-    /**
-     * Initialize the kernel.
-     *
-     * @param system     the System this kernel will be applied to
-     * @param force      the GBSAOBCForce this kernel will be used for
-     */
-    void initialize(const System& system, const GBSAOBCForce& force);
-    /**
-     * Execute the kernel to calculate the forces and/or energy.
-     *
-     * @param context        the context in which to execute this kernel
-     * @param includeForces  true if forces should be calculated
-     * @param includeEnergy  true if the energy should be calculated
-     * @return the potential energy due to the force
-     */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-private:
-    OpenCLPlatform::PlatformData& data;
-};
-
-/**
- * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
- */
-class OpenCLParallelCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
-public:
-    OpenCLParallelCalcCustomGBForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcCustomGBForceKernel(name, platform),
-            data(data) {
-    }
-    ~OpenCLParallelCalcCustomGBForceKernel();
-    /**
-     * Initialize the kernel.
-     *
-     * @param system     the System this kernel will be applied to
-     * @param force      the CustomGBForce this kernel will be used for
-     */
-    void initialize(const System& system, const CustomGBForce& force);
-    /**
-     * Execute the kernel to calculate the forces and/or energy.
-     *
-     * @param context        the context in which to execute this kernel
-     * @param includeForces  true if forces should be calculated
-     * @param includeEnergy  true if the energy should be calculated
-     * @return the potential energy due to the force
-     */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-private:
-    OpenCLPlatform::PlatformData& data;
-};
-
 /**
  * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
  */
 class OpenCLParallelCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
 public:
-    OpenCLParallelCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcCustomExternalForceKernel(name, platform),
-            data(data) {
+    OpenCLParallelCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system);
+    OpenCLCalcCustomExternalForceKernel& getKernel(int index) {
+        return dynamic_cast<OpenCLCalcCustomExternalForceKernel&>(kernels[index].getImpl());
     }
-    ~OpenCLParallelCalcCustomExternalForceKernel();
     /**
      * Initialize the kernel.
      *
@@ -460,6 +373,7 @@ public:
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
     OpenCLPlatform::PlatformData& data;
+    std::vector<Kernel> kernels;
 };
 
 /**
@@ -467,10 +381,10 @@ private:
  */
 class OpenCLParallelCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
 public:
-    OpenCLParallelCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcCustomHbondForceKernel(name, platform),
-            data(data) {
+    OpenCLParallelCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system);
+    OpenCLCalcCustomHbondForceKernel& getKernel(int index) {
+        return dynamic_cast<OpenCLCalcCustomHbondForceKernel&>(kernels[index].getImpl());
     }
-    ~OpenCLParallelCalcCustomHbondForceKernel();
     /**
      * Initialize the kernel.
      *
@@ -489,6 +403,7 @@ public:
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
     OpenCLPlatform::PlatformData& data;
+    std::vector<Kernel> kernels;
 };
 
 } // namespace OpenMM

platforms/opencl/src/OpenCLPlatform.cpp

@@ -37,6 +37,7 @@ using namespace OpenMM;
 using std::map;
 using std::string;
 using std::stringstream;
+using std::vector;
 
 extern "C" OPENMM_EXPORT void registerPlatforms() {
     Platform::registerPlatform(new OpenCLPlatform());
@@ -91,13 +92,10 @@ void OpenCLPlatform::setPropertyValue(Context& context, const string& property,
 }
 
 void OpenCLPlatform::contextCreated(ContextImpl& context, const map<string, string>& properties) const {
-    unsigned int deviceIndex = -1;
     const string& devicePropValue = (properties.find(OpenCLDeviceIndex()) == properties.end() ?
             getPropertyDefaultValue(OpenCLDeviceIndex()) : properties.find(OpenCLDeviceIndex())->second);
-    if (devicePropValue.length() > 0)
-        stringstream(devicePropValue) >> deviceIndex;
     int numParticles = context.getSystem().getNumParticles();
-    context.setPlatformData(new PlatformData(numParticles, deviceIndex));
+    context.setPlatformData(new PlatformData(numParticles, devicePropValue));
 }
 
 void OpenCLPlatform::contextDestroyed(ContextImpl& context) const {
@@ -105,10 +103,29 @@ void OpenCLPlatform::contextDestroyed(ContextImpl& context) const {
     delete data;
 }
 
-OpenCLPlatform::PlatformData::PlatformData(int numParticles, int deviceIndex) : removeCM(false), stepCount(0), computeForceCount(0), time(0.0)  {
-    contexts.push_back(new OpenCLContext(numParticles, deviceIndex, *this));
+OpenCLPlatform::PlatformData::PlatformData(int numParticles, const string& deviceIndexProperty) : removeCM(false), stepCount(0), computeForceCount(0), time(0.0)  {
+    vector<string> devices;
+    size_t searchPos = 0, nextPos;
+    while ((nextPos = deviceIndexProperty.find(',', searchPos)) != string::npos) {
+        devices.push_back(deviceIndexProperty.substr(searchPos, nextPos-searchPos));
+        searchPos = nextPos+1;
+    }
+    devices.push_back(deviceIndexProperty.substr(searchPos));
+    for (int i = 0; i < (int) devices.size(); i++) {
+        if (devices[i].length() > 0) {
+            unsigned int deviceIndex;
+            stringstream(deviceIndexProperty) >> deviceIndex;
+            contexts.push_back(new OpenCLContext(numParticles, deviceIndex, *this));
+        }
+    }
+    if (contexts.size() == 0)
+        contexts.push_back(new OpenCLContext(numParticles, -1, *this));
     stringstream device;
-    device << contexts[0]->getDeviceIndex();
+    for (int i = 0; i < (int) contexts.size(); i++) {
+        if (i > 0)
+            device << ',';
+        device << contexts[i]->getDeviceIndex();
+    }
     propertyValues[OpenCLPlatform::OpenCLDeviceIndex()] = device.str();
 }
 

platforms/opencl/tests/TestOpenCLRandom.cpp

@@ -48,7 +48,7 @@ void testGaussian() {
     System system;
     for (int i = 0; i < numAtoms; i++)
         system.addParticle(1.0);
-    OpenCLPlatform::PlatformData platformData(numAtoms, -1);
+    OpenCLPlatform::PlatformData platformData(numAtoms, "");
     OpenCLContext& context = *platformData.contexts[0];
     context.initialize(system);
     context.getIntegrationUtilities().initRandomNumberGenerator(0);

platforms/opencl/tests/TestOpenCLSort.cpp

@@ -51,7 +51,7 @@ void verifySorting(vector<float> array) {
 
     System system;
     system.addParticle(0.0);
-    OpenCLPlatform::PlatformData platformData(1, -1);
+    OpenCLPlatform::PlatformData platformData(1, "");
     OpenCLContext& context = *platformData.contexts[0];
     context.initialize(system);
     OpenCLArray<float> data(context, array.size(), "sortData");