Merge pull request #1992 from peastman/arrays

Simplification to CudaArray and OpenCLArray

platforms/cuda/include/CudaArray.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -57,6 +57,11 @@ public:
     static CudaArray* create(CudaContext& context, int size, const std::string& name) {
         return new CudaArray(context, size, sizeof(T), name);
     }
+    /**
+     * Create an uninitialized CudaArray object.  It does not point to any device memory,
+     * and cannot be used until initialize() is called on it.
+     */
+    CudaArray();
     /**
      * Create a CudaArray object.
      *
@@ -67,6 +72,36 @@ public:
      */
     CudaArray(CudaContext& context, int size, int elementSize, const std::string& name);
     ~CudaArray();
+    /**
+     * Initialize this object.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param elementSize       the size of each element in bytes
+     * @param name              the name of the array
+     */
+    void initialize(CudaContext& context, int size, int elementSize, const std::string& name);
+    /**
+     * Initialize this object.  The template argument is the data type of each array element.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param name              the name of the array
+     */
+    template <class T>
+    void initialize(CudaContext& context, int size, const std::string& name) {
+        initialize(context, size, sizeof(T), name);
+    }
+    /**
+     * Recreate the internal storage to have a different size.
+     */
+    void resize(int size);
+    /**
+     * Get whether this array has been initialized.
+     */
+    bool isInitialized() const {
+        return (pointer != 0);
+    }
     /**
      * Get the number of elements in the array.
      */
@@ -134,7 +169,7 @@ public:
      */
     void copyTo(CudaArray& dest) const;
 private:
-    CudaContext& context;
+    CudaContext* context;
     CUdeviceptr pointer;
     int size, elementSize;
     bool ownsMemory;

platforms/cuda/include/CudaBondedUtilities.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -81,7 +81,6 @@ namespace OpenMM {
 class OPENMM_EXPORT_CUDA CudaBondedUtilities {
 public:
     CudaBondedUtilities(CudaContext& context);
-    ~CudaBondedUtilities();
     /**
      * Add a bonded interaction.
      *
@@ -136,7 +135,7 @@ private:
     std::vector<int> forceGroup;
     std::vector<CUdeviceptr> arguments;
     std::vector<std::string> argTypes;
-    std::vector<std::vector<CudaArray*> > atomIndices;
+    std::vector<std::vector<CudaArray> > atomIndices;
     std::vector<std::string> prefixCode;
     std::vector<std::string> energyParameterDerivatives;
     std::vector<void*> kernelArgs;

platforms/cuda/include/CudaContext.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -41,6 +41,7 @@
 #include <builtin_types.h>
 #include <vector_functions.h>
 #include "windowsExportCuda.h"
+#include "CudaArray.h"
 #include "CudaPlatform.h"
 #include "openmm/Kernel.h"
 
@@ -48,7 +49,6 @@ typedef unsigned int tileflags;
 
 namespace OpenMM {
 
-class CudaArray;
 class CudaForceInfo;
 class CudaExpressionUtilities;
 class CudaIntegrationUtilities;
@@ -152,37 +152,37 @@ public:
      * Get the array which contains the position (the xyz components) and charge (the w component) of each atom.
      */
     CudaArray& getPosq() {
-        return *posq;
+        return posq;
     }
     /**
      * Get the array which contains a correction to the position of each atom.  This only exists if getUseMixedPrecision() returns true.
      */
     CudaArray& getPosqCorrection() {
-        return *posqCorrection;
+        return posqCorrection;
     }
     /**
      * Get the array which contains the velocity (the xyz components) and inverse mass (the w component) of each atom.
      */
     CudaArray& getVelm() {
-        return *velm;
+        return velm;
     }
     /**
      * Get the array which contains the force on each atom (represented as three long longs in 64 bit fixed point).
      */
     CudaArray& getForce() {
-        return *force;
+        return force;
     }
     /**
      * Get the array which contains the buffer in which energy is computed.
      */
     CudaArray& getEnergyBuffer() {
-        return *energyBuffer;
+        return energyBuffer;
     }
     /**
      * Get the array which contains the buffer in which derivatives of the energy with respect to parameters are computed.
      */
     CudaArray& getEnergyParamDerivBuffer() {
-        return *energyParamDerivBuffer;
+        return energyParamDerivBuffer;
     }
     /**
      * Get a pointer to a block of pinned memory that can be used for efficient transfers between host and device.
@@ -201,7 +201,7 @@ public:
      * Get the array which contains the index of each atom.
      */
     CudaArray& getAtomIndexArray() {
-        return *atomIndexDevice;
+        return atomIndexDevice;
     }
     /**
      * Get the number of cells by which the positions are offset.
@@ -649,15 +649,15 @@ private:
     std::vector<MoleculeGroup> moleculeGroups;
     std::vector<int4> posCellOffsets;
     void* pinnedBuffer;
-    CudaArray* posq;
-    CudaArray* posqCorrection;
-    CudaArray* velm;
-    CudaArray* force;
-    CudaArray* energyBuffer;
-    CudaArray* energySum;
-    CudaArray* energyParamDerivBuffer;
-    CudaArray* atomIndexDevice;
-    CudaArray* chargeBuffer;
+    CudaArray posq;
+    CudaArray posqCorrection;
+    CudaArray velm;
+    CudaArray force;
+    CudaArray energyBuffer;
+    CudaArray energySum;
+    CudaArray energyParamDerivBuffer;
+    CudaArray atomIndexDevice;
+    CudaArray chargeBuffer;
     std::vector<std::string> energyParamDerivNames;
     std::map<std::string, double> energyParamDerivWorkspace;
     std::vector<int> atomIndex;

platforms/cuda/include/CudaIntegrationUtilities.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,20 +47,20 @@ public:
      * Get the array which contains position deltas.
      */
     CudaArray& getPosDelta() {
-        return *posDelta;
+        return posDelta;
     }
     /**
      * Get the array which contains random values.  Each element is a float4, whose components
      * are independent, normally distributed random numbers with mean 0 and variance 1.
      */
     CudaArray& getRandom() {
-        return *random;
+        return random;
     }
     /**
      * Get the array which contains the current step size.
      */
     CudaArray& getStepSize() {
-        return *stepSize;
+        return stepSize;
     }
     /**
      * Set the size to use for the next step.
@@ -131,38 +131,38 @@ private:
     CUfunction ccmaUpdateKernel;
     CUfunction vsitePositionKernel, vsiteForceKernel;
     CUfunction randomKernel, timeShiftKernel;
-    CudaArray* posDelta;
-    CudaArray* settleAtoms;
-    CudaArray* settleParams;
-    CudaArray* shakeAtoms;
-    CudaArray* shakeParams;
-    CudaArray* random;
-    CudaArray* randomSeed;
-    CudaArray* stepSize;
-    CudaArray* ccmaAtoms;
-    CudaArray* ccmaDistance;
-    CudaArray* ccmaReducedMass;
-    CudaArray* ccmaAtomConstraints;
-    CudaArray* ccmaNumAtomConstraints;
-    CudaArray* ccmaConstraintMatrixColumn;
-    CudaArray* ccmaConstraintMatrixValue;
-    CudaArray* ccmaDelta1;
-    CudaArray* ccmaDelta2;
-    CudaArray* ccmaConverged;
+    CudaArray posDelta;
+    CudaArray settleAtoms;
+    CudaArray settleParams;
+    CudaArray shakeAtoms;
+    CudaArray shakeParams;
+    CudaArray random;
+    CudaArray randomSeed;
+    CudaArray stepSize;
+    CudaArray ccmaAtoms;
+    CudaArray ccmaDistance;
+    CudaArray ccmaReducedMass;
+    CudaArray ccmaAtomConstraints;
+    CudaArray ccmaNumAtomConstraints;
+    CudaArray ccmaConstraintMatrixColumn;
+    CudaArray ccmaConstraintMatrixValue;
+    CudaArray ccmaDelta1;
+    CudaArray ccmaDelta2;
+    CudaArray ccmaConverged;
     int* ccmaConvergedMemory;
     CUdeviceptr ccmaConvergedDeviceMemory;
     CUevent ccmaEvent;
-    CudaArray* vsite2AvgAtoms;
-    CudaArray* vsite2AvgWeights;
-    CudaArray* vsite3AvgAtoms;
-    CudaArray* vsite3AvgWeights;
-    CudaArray* vsiteOutOfPlaneAtoms;
-    CudaArray* vsiteOutOfPlaneWeights;
-    CudaArray* vsiteLocalCoordsIndex;
-    CudaArray* vsiteLocalCoordsAtoms;
-    CudaArray* vsiteLocalCoordsWeights;
-    CudaArray* vsiteLocalCoordsPos;
-    CudaArray* vsiteLocalCoordsStartIndex;
+    CudaArray vsite2AvgAtoms;
+    CudaArray vsite2AvgWeights;
+    CudaArray vsite3AvgAtoms;
+    CudaArray vsite3AvgWeights;
+    CudaArray vsiteOutOfPlaneAtoms;
+    CudaArray vsiteOutOfPlaneWeights;
+    CudaArray vsiteLocalCoordsIndex;
+    CudaArray vsiteLocalCoordsAtoms;
+    CudaArray vsiteLocalCoordsWeights;
+    CudaArray vsiteLocalCoordsPos;
+    CudaArray vsiteLocalCoordsStartIndex;
     int randomPos;
     int lastSeed, numVsites;
     double2 lastStepSize;

platforms/cuda/include/CudaNonbondedUtilities.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -164,61 +164,61 @@ public:
      * Get the array containing the center of each atom block.
      */
     CudaArray& getBlockCenters() {
-        return *blockCenter;
+        return blockCenter;
     }
     /**
      * Get the array containing the dimensions of each atom block.
      */
     CudaArray& getBlockBoundingBoxes() {
-        return *blockBoundingBox;
+        return blockBoundingBox;
     }
     /**
      * Get the array whose first element contains the number of tiles with interactions.
      */
     CudaArray& getInteractionCount() {
-        return *interactionCount;
+        return interactionCount;
     }
     /**
      * Get the array containing tiles with interactions.
      */
     CudaArray& getInteractingTiles() {
-        return *interactingTiles;
+        return interactingTiles;
     }
     /**
      * Get the array containing the atoms in each tile with interactions.
      */
     CudaArray& getInteractingAtoms() {
-        return *interactingAtoms;
+        return interactingAtoms;
     }
     /**
      * Get the array containing single pairs in the neighbor list.
      */
     CudaArray& getSinglePairs() {
-        return *singlePairs;
+        return singlePairs;
     }
     /**
      * Get the array containing exclusion flags.
      */
     CudaArray& getExclusions() {
-        return *exclusions;
+        return exclusions;
     }
     /**
      * Get the array containing tiles with exclusions.
      */
     CudaArray& getExclusionTiles() {
-        return *exclusionTiles;
+        return exclusionTiles;
     }
     /**
      * Get the array containing the index into the exclusion array for each tile.
      */
     CudaArray& getExclusionIndices() {
-        return *exclusionIndices;
+        return exclusionIndices;
     }
     /**
      * Get the array listing where the exclusion data starts for each row.
      */
     CudaArray& getExclusionRowIndices() {
-        return *exclusionRowIndices;
+        return exclusionRowIndices;
     }
     /**
      * Get the index of the first tile this context is responsible for processing.
@@ -270,22 +270,22 @@ private:
     class BlockSortTrait;
     CudaContext& context;
     std::map<int, KernelSet> groupKernels;
-    CudaArray* exclusionTiles;
-    CudaArray* exclusions;
-    CudaArray* exclusionIndices;
-    CudaArray* exclusionRowIndices;
-    CudaArray* interactingTiles;
-    CudaArray* interactingAtoms;
-    CudaArray* interactionCount;
-    CudaArray* singlePairs;
-    CudaArray* singlePairCount;
-    CudaArray* blockCenter;
-    CudaArray* blockBoundingBox;
-    CudaArray* sortedBlocks;
-    CudaArray* sortedBlockCenter;
-    CudaArray* sortedBlockBoundingBox;
-    CudaArray* oldPositions;
-    CudaArray* rebuildNeighborList;
+    CudaArray exclusionTiles;
+    CudaArray exclusions;
+    CudaArray exclusionIndices;
+    CudaArray exclusionRowIndices;
+    CudaArray interactingTiles;
+    CudaArray interactingAtoms;
+    CudaArray interactionCount;
+    CudaArray singlePairs;
+    CudaArray singlePairCount;
+    CudaArray blockCenter;
+    CudaArray blockBoundingBox;
+    CudaArray sortedBlocks;
+    CudaArray sortedBlockCenter;
+    CudaArray sortedBlockBoundingBox;
+    CudaArray oldPositions;
+    CudaArray rebuildNeighborList;
     CudaSort* blockSorter;
     CUevent downloadCountEvent;
     int* pinnedCountBuffer;

platforms/cuda/include/CudaParallelKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -84,7 +84,7 @@ private:
     std::vector<long long> completionTimes;
     std::vector<double> contextNonbondedFractions;
     int2* interactionCounts;
-    CudaArray* contextForces;
+    CudaArray contextForces;
     void* pinnedPositionBuffer;
     long long* pinnedForceBuffer;
     CUfunction sumKernel;

platforms/cuda/include/CudaSort.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -87,11 +87,11 @@ public:
 private:
     CudaContext& context;
     SortTrait* trait;
-    CudaArray* dataRange;
-    CudaArray* bucketOfElement;
-    CudaArray* offsetInBucket;
-    CudaArray* bucketOffset;
-    CudaArray* buckets;
+    CudaArray dataRange;
+    CudaArray bucketOfElement;
+    CudaArray offsetInBucket;
+    CudaArray bucketOffset;
+    CudaArray buckets;
     CUfunction shortListKernel, computeRangeKernel, assignElementsKernel, computeBucketPositionsKernel, copyToBucketsKernel, sortBucketsKernel;
     unsigned int dataLength, rangeKernelSize, positionsKernelSize, sortKernelSize;
     bool isShortList;

platforms/cuda/src/CudaArray.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2012 Stanford University and the Authors.           *
+ * Portions copyright (c) 2012-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,18 +32,15 @@
 
 using namespace OpenMM;
 
-CudaArray::CudaArray(CudaContext& context, int size, int elementSize, const std::string& name) :
-        context(context), size(size), elementSize(elementSize), name(name), ownsMemory(true) {
-    CUresult result = cuMemAlloc(&pointer, size*elementSize);
-    if (result != CUDA_SUCCESS) {
-        std::stringstream str;
-        str<<"Error creating array "<<name<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";
-        throw OpenMMException(str.str());
-    }
+CudaArray::CudaArray() : pointer(0), ownsMemory(false) {
+}
+
+CudaArray::CudaArray(CudaContext& context, int size, int elementSize, const std::string& name) : pointer(0) {
+    initialize(context, size, elementSize, name);
 }
 
 CudaArray::~CudaArray() {
-    if (ownsMemory && context.getContextIsValid()) {
+    if (pointer != 0 && ownsMemory && context->getContextIsValid()) {
         CUresult result = cuMemFree(pointer);
         if (result != CUDA_SUCCESS) {
             std::stringstream str;
@@ -53,12 +50,45 @@ CudaArray::~CudaArray() {
     }
 }
 
+void CudaArray::initialize(CudaContext& context, int size, int elementSize, const std::string& name) {
+    if (this->pointer != 0)
+        throw OpenMMException("CudaArray has already been initialized");
+    this->context = &context;
+    this->size = size;
+    this->elementSize = elementSize;
+    this->name = name;
+    ownsMemory = true;
+    CUresult result = cuMemAlloc(&pointer, size*elementSize);
+    if (result != CUDA_SUCCESS) {
+        std::stringstream str;
+        str<<"Error creating array "<<name<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";
+        throw OpenMMException(str.str());
+    }
+}
+
+void CudaArray::resize(int size) {
+    if (pointer == 0)
+        throw OpenMMException("CudaArray has not been initialized");
+    if (!ownsMemory)
+        throw OpenMMException("Cannot resize an array that does not own its storage");
+    CUresult result = cuMemFree(pointer);
+    if (result != CUDA_SUCCESS) {
+        std::stringstream str;
+        str<<"Error deleting array "<<name<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";
+        throw OpenMMException(str.str());
+    }
+    pointer = 0;
+    initialize(*context, size, elementSize, name);
+}
+
 void CudaArray::upload(const void* data, bool blocking) {
+    if (pointer == 0)
+        throw OpenMMException("CudaArray has not been initialized");
     CUresult result;
     if (blocking)
         result = cuMemcpyHtoD(pointer, data, size*elementSize);
     else
-        result = cuMemcpyHtoDAsync(pointer, data, size*elementSize, context.getCurrentStream());
+        result = cuMemcpyHtoDAsync(pointer, data, size*elementSize, context->getCurrentStream());
     if (result != CUDA_SUCCESS) {
         std::stringstream str;
         str<<"Error uploading array "<<name<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";
@@ -67,11 +97,13 @@ void CudaArray::upload(const void* data, bool blocking) {
 }
 
 void CudaArray::download(void* data, bool blocking) const {
+    if (pointer == 0)
+        throw OpenMMException("CudaArray has not been initialized");
     CUresult result;
     if (blocking)
         result = cuMemcpyDtoH(data, pointer, size*elementSize);
     else
-        result = cuMemcpyDtoHAsync(data, pointer, size*elementSize, context.getCurrentStream());
+        result = cuMemcpyDtoHAsync(data, pointer, size*elementSize, context->getCurrentStream());
     if (result != CUDA_SUCCESS) {
         std::stringstream str;
         str<<"Error downloading array "<<name<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";
@@ -80,9 +112,11 @@ void CudaArray::download(void* data, bool blocking) const {
 }
 
 void CudaArray::copyTo(CudaArray& dest) const {
+    if (pointer == 0)
+        throw OpenMMException("CudaArray has not been initialized");
     if (dest.getSize() != size || dest.getElementSize() != elementSize)
         throw OpenMMException("Error copying array "+name+" to "+dest.getName()+": The destination array does not match the size of the array");
-    CUresult result = cuMemcpyDtoDAsync(dest.getDevicePointer(), pointer, size*elementSize, context.getCurrentStream());
+    CUresult result = cuMemcpyDtoDAsync(dest.getDevicePointer(), pointer, size*elementSize, context->getCurrentStream());
     if (result != CUDA_SUCCESS) {
         std::stringstream str;
         str<<"Error copying array "<<name<<" to "<<dest.getName()<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")";

platforms/cuda/src/CudaBondedUtilities.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -37,12 +37,6 @@ using namespace std;
 CudaBondedUtilities::CudaBondedUtilities(CudaContext& context) : context(context), numForceBuffers(0), maxBonds(0), allGroups(0), hasInitializedKernels(false) {
 }
 
-CudaBondedUtilities::~CudaBondedUtilities() {
-    for (int i = 0; i < (int) atomIndices.size(); i++)
-        for (int j = 0; j < (int) atomIndices[i].size(); j++)
-            delete atomIndices[i][j];
-}
-
 void CudaBondedUtilities::addInteraction(const vector<vector<int> >& atoms, const string& source, int group) {
     if (atoms.size() > 0) {
         forceAtoms.push_back(atoms);
@@ -90,8 +84,10 @@ void CudaBondedUtilities::initialize(const System& system) {
     for (int i = 0; i < numForces; i++) {
         int numBonds = forceAtoms[i].size();
         int numAtoms = forceAtoms[i][0].size();
+        int numArrays = (numAtoms+3)/4;
         int startAtom = 0;
-        while (startAtom < numAtoms) {
+        atomIndices[i].resize(numArrays);
+        for (int j = 0; j < numArrays; j++) {
             int width = min(numAtoms-startAtom, 4);
             int paddedWidth = (width == 3 ? 4 : width);
             vector<unsigned int> indexVec(paddedWidth*numBonds);
@@ -99,9 +95,8 @@ void CudaBondedUtilities::initialize(const System& system) {
                 for (int atom = 0; atom < width; atom++)
                     indexVec[bond*paddedWidth+atom] = forceAtoms[i][bond][startAtom+atom];
             }
-            CudaArray* indices = new CudaArray(context, numBonds, 4*paddedWidth, "bondedIndices");
-            indices->upload(&indexVec[0]);
-            atomIndices[i].push_back(indices);
+            atomIndices[i][j].initialize(context, numBonds, 4*paddedWidth, "bondedIndices");
+            atomIndices[i][j].upload(&indexVec[0]);
             startAtom += width;
         }
     }
@@ -115,7 +110,7 @@ void CudaBondedUtilities::initialize(const System& system) {
     s<<"extern \"C\" __global__ void computeBondedForces(unsigned long long* __restrict__ forceBuffer, mixed* __restrict__ energyBuffer, const real4* __restrict__ posq, int groups, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ";
     for (int force = 0; force < numForces; force++) {
         for (int i = 0; i < (int) atomIndices[force].size(); i++) {
-            int indexWidth = atomIndices[force][i]->getElementSize()/4;
+            int indexWidth = atomIndices[force][i].getElementSize()/4;
             string indexType = "uint"+context.intToString(indexWidth);
             s<<", const "<<indexType<<"* __restrict__ atomIndices"<<force<<"_"<<i;
         }
@@ -154,7 +149,7 @@ string CudaBondedUtilities::createForceSource(int forceIndex, int numBonds, int
     s<<"for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < "<<numBonds<<"; index += blockDim.x*gridDim.x) {\n";
     int startAtom = 0;
     for (int i = 0; i < (int) atomIndices[forceIndex].size(); i++) {
-        int indexWidth = atomIndices[forceIndex][i]->getElementSize()/4;
+        int indexWidth = atomIndices[forceIndex][i].getElementSize()/4;
         string indexType = "uint"+context.intToString(indexWidth);
         s<<"    "<<indexType<<" atoms"<<i<<" = atomIndices"<<forceIndex<<"_"<<i<<"[index];\n";
         int atomsToLoad = min(indexWidth, numAtoms-startAtom);
@@ -191,7 +186,7 @@ void CudaBondedUtilities::computeInteractions(int groups) {
         kernelArgs.push_back(context.getPeriodicBoxVecZPointer());
         for (int i = 0; i < (int) atomIndices.size(); i++)
             for (int j = 0; j < (int) atomIndices[i].size(); j++)
-                kernelArgs.push_back(&atomIndices[i][j]->getDevicePointer());
+                kernelArgs.push_back(&atomIndices[i][j].getDevicePointer());
         for (int i = 0; i < (int) arguments.size(); i++)
             kernelArgs.push_back(&arguments[i]);
         if (energyParameterDerivatives.size() > 0)

platforms/cuda/src/CudaContext.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -108,8 +108,7 @@ static int executeInWindows(const string &command) {
 CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlockingSync, const string& precision, const string& compiler,
         const string& tempDir, const std::string& hostCompiler, CudaPlatform::PlatformData& platformData, CudaContext* originalContext) : system(system), currentStream(0),
         time(0.0), platformData(platformData), stepCount(0), computeForceCount(0), stepsSinceReorder(99999), contextIsValid(false), atomsWereReordered(false), hasCompilerKernel(false), isNvccAvailable(false),
-        pinnedBuffer(NULL), posq(NULL), posqCorrection(NULL), velm(NULL), force(NULL), energyBuffer(NULL), energySum(NULL), energyParamDerivBuffer(NULL), atomIndexDevice(NULL), chargeBuffer(NULL),
-        integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
+        pinnedBuffer(NULL), integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
     // Determine what compiler to use.
     
     this->compiler = "\""+compiler+"\"";
@@ -268,8 +267,8 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
         compilationDefines["BALLOT(var)"] = "__ballot(var);";
     }
     if (useDoublePrecision) {
-        posq = CudaArray::create<double4>(*this, paddedNumAtoms, "posq");
-        velm = CudaArray::create<double4>(*this, paddedNumAtoms, "velm");
+        posq.initialize<double4>(*this, paddedNumAtoms, "posq");
+        velm.initialize<double4>(*this, paddedNumAtoms, "velm");
         compilationDefines["USE_DOUBLE_PRECISION"] = "1";
         compilationDefines["make_real2"] = "make_double2";
         compilationDefines["make_real3"] = "make_double3";
@@ -279,9 +278,9 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
         compilationDefines["make_mixed4"] = "make_double4";
     }
     else if (useMixedPrecision) {
-        posq = CudaArray::create<float4>(*this, paddedNumAtoms, "posq");
-        posqCorrection = CudaArray::create<float4>(*this, paddedNumAtoms, "posqCorrection");
-        velm = CudaArray::create<double4>(*this, paddedNumAtoms, "velm");
+        posq.initialize<float4>(*this, paddedNumAtoms, "posq");
+        posqCorrection.initialize<float4>(*this, paddedNumAtoms, "posqCorrection");
+        velm.initialize<double4>(*this, paddedNumAtoms, "velm");
         compilationDefines["USE_MIXED_PRECISION"] = "1";
         compilationDefines["make_real2"] = "make_float2";
         compilationDefines["make_real3"] = "make_float3";
@@ -291,8 +290,8 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
         compilationDefines["make_mixed4"] = "make_double4";
     }
     else {
-        posq = CudaArray::create<float4>(*this, paddedNumAtoms, "posq");
-        velm = CudaArray::create<float4>(*this, paddedNumAtoms, "velm");
+        posq.initialize<float4>(*this, paddedNumAtoms, "posq");
+        velm.initialize<float4>(*this, paddedNumAtoms, "velm");
         compilationDefines["make_real2"] = "make_float2";
         compilationDefines["make_real3"] = "make_float3";
         compilationDefines["make_real4"] = "make_float4";
@@ -415,24 +414,6 @@ CudaContext::~CudaContext() {
         delete computation;
     if (pinnedBuffer != NULL)
         cuMemFreeHost(pinnedBuffer);
-    if (posq != NULL)
-        delete posq;
-    if (posqCorrection != NULL)
-        delete posqCorrection;
-    if (velm != NULL)
-        delete velm;
-    if (force != NULL)
-        delete force;
-    if (energyBuffer != NULL)
-        delete energyBuffer;
-    if (energySum != NULL)
-        delete energySum;
-    if (energyParamDerivBuffer != NULL)
-        delete energyParamDerivBuffer;
-    if (atomIndexDevice != NULL)
-        delete atomIndexDevice;
-    if (chargeBuffer != NULL)
-        delete chargeBuffer;
     if (integration != NULL)
         delete integration;
     if (expression != NULL)
@@ -456,20 +437,20 @@ void CudaContext::initialize() {
     string errorMessage = "Error initializing Context";
     int numEnergyBuffers = max(numThreadBlocks*ThreadBlockSize, nonbonded->getNumEnergyBuffers());
     if (useDoublePrecision) {
-        energyBuffer = CudaArray::create<double>(*this, numEnergyBuffers, "energyBuffer");
-        energySum = CudaArray::create<double>(*this, 1, "energySum");
+        energyBuffer.initialize<double>(*this, numEnergyBuffers, "energyBuffer");
+        energySum.initialize<double>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
     }
     else if (useMixedPrecision) {
-        energyBuffer = CudaArray::create<double>(*this, numEnergyBuffers, "energyBuffer");
-        energySum = CudaArray::create<double>(*this, 1, "energySum");
+        energyBuffer.initialize<double>(*this, numEnergyBuffers, "energyBuffer");
+        energySum.initialize<double>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
     }
     else {
-        energyBuffer = CudaArray::create<float>(*this, numEnergyBuffers, "energyBuffer");
-        energySum = CudaArray::create<float>(*this, 1, "energySum");
+        energyBuffer.initialize<float>(*this, numEnergyBuffers, "energyBuffer");
+        energySum.initialize<float>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*6, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(float), 0));
     }
@@ -480,24 +461,24 @@ void CudaContext::initialize() {
         else
             ((float4*) pinnedBuffer)[i] = make_float4(0.0f, 0.0f, 0.0f, mass == 0.0 ? 0.0f : (float) (1.0/mass));
     }
-    velm->upload(pinnedBuffer);
+    velm.upload(pinnedBuffer);
     bonded->initialize(system);
-    force = CudaArray::create<long long>(*this, paddedNumAtoms*3, "force");
-    addAutoclearBuffer(force->getDevicePointer(), force->getSize()*force->getElementSize());
-    addAutoclearBuffer(energyBuffer->getDevicePointer(), energyBuffer->getSize()*energyBuffer->getElementSize());
+    force.initialize<long long>(*this, paddedNumAtoms*3, "force");
+    addAutoclearBuffer(force.getDevicePointer(), force.getSize()*force.getElementSize());
+    addAutoclearBuffer(energyBuffer.getDevicePointer(), energyBuffer.getSize()*energyBuffer.getElementSize());
     int numEnergyParamDerivs = energyParamDerivNames.size();
     if (numEnergyParamDerivs > 0) {
         if (useDoublePrecision || useMixedPrecision)
-            energyParamDerivBuffer = CudaArray::create<double>(*this, numEnergyParamDerivs*numEnergyBuffers, "energyParamDerivBuffer");
+            energyParamDerivBuffer.initialize<double>(*this, numEnergyParamDerivs*numEnergyBuffers, "energyParamDerivBuffer");
         else
-            energyParamDerivBuffer = CudaArray::create<float>(*this, numEnergyParamDerivs*numEnergyBuffers, "energyParamDerivBuffer");
-        addAutoclearBuffer(*energyParamDerivBuffer);
+            energyParamDerivBuffer.initialize<float>(*this, numEnergyParamDerivs*numEnergyBuffers, "energyParamDerivBuffer");
+        addAutoclearBuffer(energyParamDerivBuffer);
     }
-    atomIndexDevice = CudaArray::create<int>(*this, paddedNumAtoms, "atomIndex");
+    atomIndexDevice.initialize<int>(*this, paddedNumAtoms, "atomIndex");
     atomIndex.resize(paddedNumAtoms);
     for (int i = 0; i < paddedNumAtoms; ++i)
         atomIndex[i] = i;
-    atomIndexDevice->upload(atomIndex);
+    atomIndexDevice.upload(atomIndex);
     findMoleculeGroups();
     nonbonded->initialize(system);
 }
@@ -890,11 +871,11 @@ void CudaContext::clearAutoclearBuffers() {
 }
 
 double CudaContext::reduceEnergy() {
-    int bufferSize = energyBuffer->getSize();
+    int bufferSize = energyBuffer.getSize();
     int workGroupSize  = 512;
-    void* args[] = {&energyBuffer->getDevicePointer(), &energySum->getDevicePointer(), &bufferSize, &workGroupSize};
-    executeKernel(reduceEnergyKernel, args, workGroupSize, workGroupSize, workGroupSize*energyBuffer->getElementSize());
-    energySum->download(pinnedBuffer);
+    void* args[] = {&energyBuffer.getDevicePointer(), &energySum.getDevicePointer(), &bufferSize, &workGroupSize};
+    executeKernel(reduceEnergyKernel, args, workGroupSize, workGroupSize, workGroupSize*energyBuffer.getElementSize());
+    energySum.download(pinnedBuffer);
     if (getUseDoublePrecision() || getUseMixedPrecision())
         return *((double*) pinnedBuffer);
     else
@@ -902,21 +883,21 @@ double CudaContext::reduceEnergy() {
 }
 
 void CudaContext::setCharges(const vector<double>& charges) {
-    if (chargeBuffer == NULL)
-        chargeBuffer = new CudaArray(*this, numAtoms, useDoublePrecision ? sizeof(double) : sizeof(float), "chargeBuffer");
+    if (!chargeBuffer.isInitialized())
+        chargeBuffer.initialize(*this, numAtoms, useDoublePrecision ? sizeof(double) : sizeof(float), "chargeBuffer");
     if (getUseDoublePrecision()) {
         double* c = (double*) getPinnedBuffer();
         for (int i = 0; i < charges.size(); i++)
             c[i] = charges[i];
-        chargeBuffer->upload(c);
+        chargeBuffer.upload(c);
     }
     else {
         float* c = (float*) getPinnedBuffer();
         for (int i = 0; i < charges.size(); i++)
             c[i] = (float) charges[i];
-        chargeBuffer->upload(c);
+        chargeBuffer.upload(c);
     }
-    void* args[] = {&chargeBuffer->getDevicePointer(), &posq->getDevicePointer(), &atomIndexDevice->getDevicePointer(), &numAtoms};
+    void* args[] = {&chargeBuffer.getDevicePointer(), &posq.getDevicePointer(), &atomIndexDevice.getDevicePointer(), &numAtoms};
     executeKernel(setChargesKernel, args, numAtoms);
 }
 
@@ -1178,16 +1159,16 @@ bool CudaContext::invalidateMolecules(CudaForceInfo* force) {
         vector<double4> newPosq(paddedNumAtoms, make_double4(0, 0, 0, 0));
         vector<double4> oldVelm(paddedNumAtoms);
         vector<double4> newVelm(paddedNumAtoms, make_double4(0, 0, 0, 0));
-        posq->download(oldPosq);
-        velm->download(oldVelm);
+        posq.download(oldPosq);
+        velm.download(oldVelm);
         for (int i = 0; i < numAtoms; i++) {
             int index = atomIndex[i];
             newPosq[index] = oldPosq[i];
             newVelm[index] = oldVelm[i];
             newCellOffsets[index] = posCellOffsets[i];
         }
-        posq->upload(newPosq);
-        velm->upload(newVelm);
+        posq.upload(newPosq);
+        velm.upload(newVelm);
     }
     else if (useMixedPrecision) {
         vector<float4> oldPosq(paddedNumAtoms);
@@ -1196,8 +1177,8 @@ bool CudaContext::invalidateMolecules(CudaForceInfo* force) {
         vector<float4> newPosqCorrection(paddedNumAtoms, make_float4(0, 0, 0, 0));
         vector<double4> oldVelm(paddedNumAtoms);
         vector<double4> newVelm(paddedNumAtoms, make_double4(0, 0, 0, 0));
-        posq->download(oldPosq);
-        velm->download(oldVelm);
+        posq.download(oldPosq);
+        velm.download(oldVelm);
         for (int i = 0; i < numAtoms; i++) {
             int index = atomIndex[i];
             newPosq[index] = oldPosq[i];
@@ -1205,31 +1186,31 @@ bool CudaContext::invalidateMolecules(CudaForceInfo* force) {
             newVelm[index] = oldVelm[i];
             newCellOffsets[index] = posCellOffsets[i];
         }
-        posq->upload(newPosq);
-        posqCorrection->upload(newPosqCorrection);
-        velm->upload(newVelm);
+        posq.upload(newPosq);
+        posqCorrection.upload(newPosqCorrection);
+        velm.upload(newVelm);
     }
     else {
         vector<float4> oldPosq(paddedNumAtoms);
         vector<float4> newPosq(paddedNumAtoms, make_float4(0, 0, 0, 0));
         vector<float4> oldVelm(paddedNumAtoms);
         vector<float4> newVelm(paddedNumAtoms, make_float4(0, 0, 0, 0));
-        posq->download(oldPosq);
-        velm->download(oldVelm);
+        posq.download(oldPosq);
+        velm.download(oldVelm);
         for (int i = 0; i < numAtoms; i++) {
             int index = atomIndex[i];
             newPosq[index] = oldPosq[i];
             newVelm[index] = oldVelm[i];
             newCellOffsets[index] = posCellOffsets[i];
         }
-        posq->upload(newPosq);
-        velm->upload(newVelm);
+        posq.upload(newPosq);
+        velm.upload(newVelm);
     }
     for (int i = 0; i < numAtoms; i++) {
         atomIndex[i] = i;
         posCellOffsets[i] = newCellOffsets[i];
     }
-    atomIndexDevice->upload(atomIndex);
+    atomIndexDevice.upload(atomIndex);
     findMoleculeGroups();
     for (auto listener : reorderListeners)
         listener->execute();
@@ -1262,10 +1243,10 @@ void CudaContext::reorderAtomsImpl() {
     vector<Real4> oldPosqCorrection(paddedNumAtoms, padding);
     Mixed4 paddingMixed = {0, 0, 0, 0};
     vector<Mixed4> oldVelm(paddedNumAtoms, paddingMixed);
-    posq->download(oldPosq);
-    velm->download(oldVelm);
+    posq.download(oldPosq);
+    velm.download(oldVelm);
     if (useMixedPrecision)
-        posqCorrection->download(oldPosqCorrection);
+        posqCorrection.download(oldPosqCorrection);
     Real minx = oldPosq[0].x, maxx = oldPosq[0].x;
     Real miny = oldPosq[0].y, maxy = oldPosq[0].y;
     Real minz = oldPosq[0].z, maxz = oldPosq[0].z;
@@ -1409,11 +1390,11 @@ void CudaContext::reorderAtomsImpl() {
         atomIndex[i] = originalIndex[i];
         posCellOffsets[i] = newCellOffsets[i];
     }
-    posq->upload(newPosq);
+    posq.upload(newPosq);
     if (useMixedPrecision)
-        posqCorrection->upload(newPosqCorrection);
-    velm->upload(newVelm);
-    atomIndexDevice->upload(atomIndex);
+        posqCorrection.upload(newPosqCorrection);
+    velm.upload(newVelm);
+    atomIndexDevice.upload(atomIndex);
     for (auto listener : reorderListeners)
         listener->execute();
 }