hdiff output

r30075/cost_function.cu 2016-03-16 18:34:30.343644700 +0000 r30074/cost_function.cu 2016-03-16 18:34:31.035651762 +0000
 36:         , m_d_rigidMaxSite(0) 36:         , m_d_rigidMaxSite(0)
 37:         , m_d_coords(0) 37:         , m_d_coords(0)
 38:         , m_nSecDiag(nSecDiag) 38:         , m_nSecDiag(nSecDiag)
 39:         , m_isAtomisticNotRigid(isAtomisticNotRigid) 39:         , m_isAtomisticNotRigid(isAtomisticNotRigid)
 40:         , m_aaConvThreshold(aaConvThreshold) 40:         , m_aaConvThreshold(aaConvThreshold)
 41:         , m_nCalls(0) 41:         , m_nCalls(0)
 42:         , m_coldFusionLim(coldFusionLim) 42:         , m_coldFusionLim(coldFusionLim)
 43:         , m_shouldFreeze(shouldFreeze) 43:         , m_shouldFreeze(shouldFreeze)
 44:         , m_isAtomFrozen(isAtomFrozen) 44:         , m_isAtomFrozen(isAtomFrozen)
 45:         , m_nFreeze(nFreeze) 45:         , m_nFreeze(nFreeze)
 46:         , m_nLargeRigidBody(0) 
 47:         , m_debugPrinting(debugPrinting) 46:         , m_debugPrinting(debugPrinting)
 48:         , m_cublas(cublas) 47:         , m_cublas(cublas)
 49:         , m_isBfgsts(false) 48:         , m_isBfgsts(false)
 50:           , m_isAaConvergence(isAaConvergence) 49:           , m_isAaConvergence(isAaConvergence)
 51: { 50: {
 52:         for (int i = 0; i < m_nRigidBody; ++i) { 
 53:                 if (m_nRigidSitesPerBody[i] > 32) { 
 54:                         m_nLargeRigidBody += 1; 
 55:                 } 
 56:         } 
 57:  
 58:         m_largeRigidIndices = new int[m_nLargeRigidBody]; 
 59:  
 60:         double *zerosx; 51:         double *zerosx;
 61:         double *zerosy; 52:         double *zerosy;
 62:         double *zerosz; 53:         double *zerosz;
 63:  54: 
 64:         zerosx = new double[numDimensions]; 55:         zerosx = new double[numDimensions];
 65:         zerosy = new double[numDimensions]; 56:         zerosy = new double[numDimensions];
 66:         zerosz = new double[numDimensions]; 57:         zerosz = new double[numDimensions];
 67:  58: 
 68:         CudaSafeCall( cudaMalloc(&m_d_zerosx, m_numDimensions * sizeof(double)) ); 59:         CudaSafeCall( cudaMalloc(&m_d_zerosx, m_numDimensions * sizeof(double)) );
 69:         CudaSafeCall( cudaMalloc(&m_d_zerosy, m_numDimensions * sizeof(double)) ); 60:         CudaSafeCall( cudaMalloc(&m_d_zerosy, m_numDimensions * sizeof(double)) );
 86:         CudaSafeCall( cudaMalloc(&m_d_gkRigid, m_nDegFreedom * sizeof(double)) ); 77:         CudaSafeCall( cudaMalloc(&m_d_gkRigid, m_nDegFreedom * sizeof(double)) );
 87:  78: 
 88:         if (m_coords != NULL) { 79:         if (m_coords != NULL) {
 89:                 CudaSafeCall( cudaMalloc(&m_d_coords,   numDimensions * sizeof(double)) ); 80:                 CudaSafeCall( cudaMalloc(&m_d_coords,   numDimensions * sizeof(double)) );
 90:                 CudaSafeCall( cudaMalloc(&m_d_energy, sizeof(double)) ); 81:                 CudaSafeCall( cudaMalloc(&m_d_energy, sizeof(double)) );
 91:                 CudaSafeCall( cudaMalloc(&m_d_gradient, numDimensions * sizeof(double)) ); 82:                 CudaSafeCall( cudaMalloc(&m_d_gradient, numDimensions * sizeof(double)) );
 92:         } 83:         }
 93:  84: 
 94:         CudaSafeCall( cudaMalloc(&m_d_isAtomFrozen, numDimensions/3 * sizeof(bool)) ); 85:         CudaSafeCall( cudaMalloc(&m_d_isAtomFrozen, numDimensions/3 * sizeof(bool)) );
 95:  86: 
 96:         CudaSafeCall( cudaMalloc(&m_d_nLargeRigidBody, sizeof(int)) ); 
 97:         CudaSafeCall( cudaMalloc(&m_d_largeRigidIndices, m_nLargeRigidBody * sizeof(int)) ); 
 98:  
 99:         int counter = 0; 
100:         while (counter < m_nLargeRigidBody) { 
101:                 for (int i = 0; i < m_nRigidBody; ++i) { 
102:                         if (m_nRigidSitesPerBody[i] > 32) { 
103:                                 m_largeRigidIndices[counter] = i; 
104:                                 ++counter; 
105:                         } 
106:                 } 
107:         } 
108:  
109:         for (size_t i = 1; i < (numDimensions + 1); ++i) { 87:         for (size_t i = 1; i < (numDimensions + 1); ++i) {
110:                 if ((i + 2)%3 == 0) { 88:                 if ((i + 2)%3 == 0) {
111:                         zerosx[i-1] = 1.0; 89:                         zerosx[i-1] = 1.0;
112:                 } 90:                 }
113:                 else { 91:                 else {
114:                         zerosx[i-1] = 0.0; 92:                         zerosx[i-1] = 0.0;
115:                 } 93:                 }
116:         } 94:         }
117:  95: 
118:         for (size_t i = 1; i < (numDimensions + 1); ++i) { 96:         for (size_t i = 1; i < (numDimensions + 1); ++i) {
145:         if (m_coords != NULL) {123:         if (m_coords != NULL) {
146:                 CudaSafeCall( cudaMemcpy(m_d_coords, m_coords, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );124:                 CudaSafeCall( cudaMemcpy(m_d_coords, m_coords, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
147:         }125:         }
148: 126: 
149:         CudaSafeCall( cudaMemcpy(m_d_isAtomFrozen, m_isAtomFrozen, numDimensions/3 * sizeof(bool), cudaMemcpyHostToDevice) );127:         CudaSafeCall( cudaMemcpy(m_d_isAtomFrozen, m_isAtomFrozen, numDimensions/3 * sizeof(bool), cudaMemcpyHostToDevice) );
150: 128: 
151:         CudaSafeCall( cudaMemcpy(m_d_zerosx, zerosx, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );129:         CudaSafeCall( cudaMemcpy(m_d_zerosx, zerosx, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
152:         CudaSafeCall( cudaMemcpy(m_d_zerosy, zerosy, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );130:         CudaSafeCall( cudaMemcpy(m_d_zerosy, zerosy, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
153:         CudaSafeCall( cudaMemcpy(m_d_zerosz, zerosz, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );131:         CudaSafeCall( cudaMemcpy(m_d_zerosz, zerosz, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
154: 132: 
155:         CudaSafeCall( cudaMemcpy(m_d_nLargeRigidBody, &m_nLargeRigidBody, sizeof(int), cudaMemcpyHostToDevice) ); 
156:         CudaSafeCall( cudaMemcpy(m_d_largeRigidIndices, m_largeRigidIndices, m_nLargeRigidBody * sizeof(int), cudaMemcpyHostToDevice) ); 
157:  
158:         delete [] zerosx;133:         delete [] zerosx;
159:         delete [] zerosy;134:         delete [] zerosy;
160:         delete [] zerosz;135:         delete [] zerosz;
 136: 
161: }137: }
162: 138: 
163: 139: 
164: 140: 
165: CostFunction::~CostFunction() 141: CostFunction::~CostFunction() 
166: {142: {
167:         CudaSafeCall( cudaFree(m_d_zerosx) );143:         CudaSafeCall( cudaFree(m_d_zerosx) );
168:         CudaSafeCall( cudaFree(m_d_zerosy) );144:         CudaSafeCall( cudaFree(m_d_zerosy) );
169:         CudaSafeCall( cudaFree(m_d_zerosz) );145:         CudaSafeCall( cudaFree(m_d_zerosz) );
170: 146: 
184:         CudaSafeCall( cudaFree(m_d_nRigidBody) );160:         CudaSafeCall( cudaFree(m_d_nRigidBody) );
185:         CudaSafeCall( cudaFree(m_d_rigidMaxSite) );161:         CudaSafeCall( cudaFree(m_d_rigidMaxSite) );
186: 162: 
187:         if (m_coords != NULL) {163:         if (m_coords != NULL) {
188:                 CudaSafeCall( cudaFree(m_d_coords) );164:                 CudaSafeCall( cudaFree(m_d_coords) );
189:                 CudaSafeCall( cudaFree(m_d_energy) );165:                 CudaSafeCall( cudaFree(m_d_energy) );
190:                 CudaSafeCall( cudaFree(m_d_gradient) );166:                 CudaSafeCall( cudaFree(m_d_gradient) );
191:         }167:         }
192: 168: 
193:         CudaSafeCall( cudaFree(m_d_isAtomFrozen) );169:         CudaSafeCall( cudaFree(m_d_isAtomFrozen) );
194:  
195:         CudaSafeCall( cudaFree(m_d_nLargeRigidBody) ); 
196:         CudaSafeCall( cudaFree(m_d_largeRigidIndices) ); 
197:  
198:         delete [] m_largeRigidIndices; 
199: }170: }
200: 171: 
201: 172: 
202: 173: 
203: // Get/set functions. 174: // Get/set functions. 
204: 175: 
205: size_t CostFunction::getNumDimensions() const176: size_t CostFunction::getNumDimensions() const
206: {177: {
207:         return m_numDimensions;178:         return m_numDimensions;
208: }179: }


r30075/cost_function.h 2016-03-16 18:34:30.743648819 +0000 r30074/cost_function.h 2016-03-16 18:34:31.419655772 +0000
100:                 bool          m_isBfgsts;             // True if costFunction is being used in a Bfgsts calculation - stops aaConvergence call.100:                 bool          m_isBfgsts;             // True if costFunction is being used in a Bfgsts calculation - stops aaConvergence call.
101:                 bool          m_isRayleighRitz;       // If true, curvature of PES calculated along direction vec at point coords (R-R min).101:                 bool          m_isRayleighRitz;       // If true, curvature of PES calculated along direction vec at point coords (R-R min).
102: 102: 
103:                 const bool    m_isAtomisticNotRigid;  // If false, use rigid body coordinates.103:                 const bool    m_isAtomisticNotRigid;  // If false, use rigid body coordinates.
104:                 const bool    m_isAaConvergence;      // RBF: if true, use more accurate method of calculating RMS force.104:                 const bool    m_isAaConvergence;      // RBF: if true, use more accurate method of calculating RMS force.
105: 105: 
106:                 const double *m_sitesRigidBody;       // RBF: coordinates of the rigid body sites.106:                 const double *m_sitesRigidBody;       // RBF: coordinates of the rigid body sites.
107:                 const double *m_rigidInverse;         // RBF: inverse eigenvalues of the unweighted tensor of gyration. 107:                 const double *m_rigidInverse;         // RBF: inverse eigenvalues of the unweighted tensor of gyration. 
108: 108: 
109:                 const int     m_nDegFreedom;          // RBF: no. of degrees of freedom.109:                 const int     m_nDegFreedom;          // RBF: no. of degrees of freedom.
110:                 const int     m_nRigidBody;           // RBF: no. of rigid bodies.110:                 const int     m_nRigidBody;           // RBF: no. of rigid bodies
111:                 const int     m_rigidMaxSite;         // RBF: max. no. of sites in a rigid body.111:                 const int     m_rigidMaxSite;         // RBF: max. no. of sites in a rigid body.
112: 112: 
113:                 const int    *m_nRigidSitesPerBody;   // RBF: no. of rigid body sites.113:                 const int    *m_nRigidSitesPerBody;   // RBF: no. of rigid body sites.
114:                 const int    *m_rigidGroups;          // RBF: list of atoms in rigid bodies.114:                 const int    *m_rigidGroups;          // RBF: list of atoms in rigid bodies.
115:                 const int    *m_rigidSingles;         // RBF: list of atoms not in rigid bodies.115:                 const int    *m_rigidSingles;         // RBF: list of atoms not in rigid bodies.
116: 116: 
117:                 int           m_nLargeRigidBody;      // RBF: no. of rigid bodies with more than 32 sites. 
118:  
119:                 int          *m_largeRigidIndices;    // RBF: array containing indices for large bodies in nRigidSitesPerBody array.  
120:  
121: 117: 
122:                 double       *m_d_coords;             // The current coordinates in Bfgsts on GPU, used in R-R minimization. 118:                 double       *m_d_coords;             // The current coordinates in Bfgsts on GPU, used in R-R minimization. 
123:                 double       *m_d_energy;             // The energy at the current coordinates in Bfgsts on GPU, used in R-R minimization. 119:                 double       *m_d_energy;             // The energy at the current coordinates in Bfgsts on GPU, used in R-R minimization. 
124:                 double       *m_d_gradient;           // The gradient at the current coordinates in Bfgsts on GPU, used in R-R minimization. 120:                 double       *m_d_gradient;           // The gradient at the current coordinates in Bfgsts on GPU, used in R-R minimization. 
125:                 double       *m_d_zerosx;             // Array used in cuBLAS operations with ones at all the x positions. 121:                 double       *m_d_zerosx;             // Array used in cuBLAS operations with ones at all the x positions. 
126:                 double       *m_d_zerosy;             // Array used in cuBLAS operations with ones at all the y positions.122:                 double       *m_d_zerosy;             // Array used in cuBLAS operations with ones at all the y positions.
127:                 double       *m_d_zerosz;             // Array used in cuBLAS operations with ones at all the z positions.123:                 double       *m_d_zerosz;             // Array used in cuBLAS operations with ones at all the z positions.
128: 124: 
129:                 bool         *m_d_isAtomFrozen;       // Logical array specifying frozen atoms, on GPU.125:                 bool         *m_d_isAtomFrozen;       // Logical array specifying frozen atoms, on GPU.
130: 126: 
132:                 double       *m_d_rigidInverse;       // RBF: IINVERSE from genrigid, used in aaconvergence subroutine, on GPU.128:                 double       *m_d_rigidInverse;       // RBF: IINVERSE from genrigid, used in aaconvergence subroutine, on GPU.
133:                 double       *m_d_xRigid;             // RBF: rigid body coordinates, on GPU.129:                 double       *m_d_xRigid;             // RBF: rigid body coordinates, on GPU.
134:                 double       *m_d_gkRigid;            // RBF: rigid body gradient, on GPU.130:                 double       *m_d_gkRigid;            // RBF: rigid body gradient, on GPU.
135:                 double       *m_d_gkAtoms;            // RBF: copy of atomistic gradient for use with aaConvergence, on GPU.131:                 double       *m_d_gkAtoms;            // RBF: copy of atomistic gradient for use with aaConvergence, on GPU.
136: 132: 
137:                 int          *m_d_nRigidSitesPerBody; // RBF: no. of rigid body sites, on GPU. 133:                 int          *m_d_nRigidSitesPerBody; // RBF: no. of rigid body sites, on GPU. 
138:                 int          *m_d_rigidGroups;        // RBF: list of atoms in rigid bodies, on GPU. 134:                 int          *m_d_rigidGroups;        // RBF: list of atoms in rigid bodies, on GPU. 
139:                 int          *m_d_rigidSingles;       // RBF: list of atoms not in rigid bodies, on GPU. 135:                 int          *m_d_rigidSingles;       // RBF: list of atoms not in rigid bodies, on GPU. 
140:                 int          *m_d_nRigidBody;         // RBF: no. of rigid bodies, on GPU. 136:                 int          *m_d_nRigidBody;         // RBF: no. of rigid bodies, on GPU. 
141:                 int          *m_d_rigidMaxSite;       // RBF: maximum number of sites in a rigid body, on GPU. 137:                 int          *m_d_rigidMaxSite;       // RBF: maximum number of sites in a rigid body, on GPU. 
142:                 int          *m_d_largeRigidIndices;  // RBF: array containing indices for large bodies in nRigidSitesPerBody array, on GPU.  
143:                 int          *m_d_nLargeRigidBody;    // RBF: no. of rigid bodies with more than 32 sites, on GPU.  
144: 138: 
145: 139: 
146:                 // Rigid body gradient transformations either side of energy/gradient call.140:                 // Rigid body gradient transformations either side of energy/gradient call.
147:                 void rigidTransforms(double *d_x, double *d_f, double *d_gradf);141:                 void rigidTransforms(double *d_x, double *d_f, double *d_gradf);
148: };142: };
149: 143: 
150: #endif /* end of include guard: COST_FUNCTION_H */144: #endif /* end of include guard: COST_FUNCTION_H */


r30075/rigid_bodies.cu 2016-03-16 18:34:30.547646801 +0000 r30074/rigid_bodies.cu 2016-03-16 18:34:31.235653880 +0000
  6:   6: 
  7: #include "cost_function.h"  7: #include "cost_function.h"
  8:   8: 
  9: namespace gpu_rigid_bodies  9: namespace gpu_rigid_bodies
 10: { 10: {
 11:         // Variables on the GPU.  11:         // Variables on the GPU. 
 12:  12: 
 13:         __device__ int  singlesThreads; 13:         __device__ int  singlesThreads;
 14:         __device__ int  thisBody; 14:         __device__ int  thisBody;
 15:         __device__ int  arraySize; 15:         __device__ int  arraySize;
 16:         __device__ int  gridSize; 
 17:         __device__ int  roundMaxSite; 
 18:         __device__ int  outSize; 
 19:         __device__ int  nBlocks; 
 20:         __device__ bool shouldFindDeriv; 16:         __device__ bool shouldFindDeriv;
 21:  17: 
 22:  18: 
 23:         // Kernels.  19:         // Kernels. 
 24:  20: 
 25:         __global__ void transform(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi); 21:         __global__ void transform(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi);
 26:  22: 
 27:         __global__ void transform2(double *d_x, const double *m_d_xRigid, const int *m_d_nRigidBody,  23:         __global__ void transform2(double *d_x, const double *m_d_xRigid, const int *m_d_nRigidBody, 
 28:                         const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  24:                         const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
 29:                         const double *m_d_sitesRigidBody, const double *d_grmi, const int *m_d_rigidMaxSite); 25:                         const double *m_d_sitesRigidBody, const double *d_grmi, const int *m_d_rigidMaxSite);
 30:  26: 
 31:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody); 27:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody);
 32:  28: 
 33:         __global__ void gradTransform1a(const double *m_d_xRigid, double *d_grmi1, double *d_grmi2, double *d_grmi3,  29:         __global__ void gradTransform1a(const double *m_d_xRigid, double *d_grmi1, double *d_grmi2, double *d_grmi3, 
 34:                         const int *m_d_nRigidBody); 30:                         const int *m_d_nRigidBody);
 35:  31: 
 36:         __global__ void gradTransform1b(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi0); 32:         __global__ void gradTransform1b(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi0);
 37:  33: 
  34:         __global__ void zeros1(double *d_xZeroOne, double *d_yZeroOne, double *d_zZeroOne);
  35: 
  36:         __global__ void ones1(const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, double *d_xZeroOne, 
  37:                         double *d_yZeroOne, double *d_zZeroOne, const int *m_d_rigidMaxSite);
  38: 
 38:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1,  39:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1, 
 39:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody,  40:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody, 
 40:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody,  41:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody, 
 41:                         const int *m_d_rigidMaxSite); 42:                         const int *m_d_rigidMaxSite);
 42:  43: 
  44:         __global__ void zeros2(double *d_xZeroOne2, double *d_yZeroOne2, double *d_zZeroOne2);
  45: 
  46:         __global__ void ones2(const int *m_d_nRigidSitesPerBody, double *d_xZeroOne2, double *d_yZeroOne2, 
  47:                         double *d_zZeroOne2, const int *m_d_rigidMaxSite);
  48: 
 43:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles,  49:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 
 44:                         const int *m_d_nRigidBody); 50:                         const int *m_d_nRigidBody);
 45:  51: 
  52:         __global__ void ones3(double *d_xZeroOne2, double *d_yZeroOne2, double *d_zZeroOne2, const int *m_d_rigidMaxSite);
  53: 
 46:         __global__ void aaConvRmdrvt(double *d_rmi10, double *d_rmi20, double *d_rmi30); 54:         __global__ void aaConvRmdrvt(double *d_rmi10, double *d_rmi20, double *d_rmi30);
 47:  55: 
 48:         __global__ void intermediate2(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi0,  56:         __global__ void intermediate2(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi0, 
 49:                         const double *d_grmi10, const double *d_grmi20, const double *d_grmi30,  57:                         const double *d_grmi10, const double *d_grmi20, const double *d_grmi30, 
 50:                         const double *m_d_sitesRigidBody, const int *m_d_rigidGroups, double *d_tempArray,  58:                         const double *m_d_sitesRigidBody, const int *m_d_rigidGroups, double *d_tempArray, 
 51:                         const int *m_d_nRigidBody, const int *m_d_rigidMaxSite); 59:                         const int *m_d_nRigidBody, const int *m_d_rigidMaxSite);
 52:  60: 
 53:         __global__ void aaConvTorque(const double *d_torques, const double *m_d_gkRigid, const double *d_grmi0,  61:         __global__ void aaConvTorque(const double *d_torques, const double *m_d_gkRigid, const double *d_grmi0, 
 54:                         const int *m_d_nRigidSitesPerBody, const double *m_d_rigidInverse, double *d_rmsArray,  62:                         const int *m_d_nRigidSitesPerBody, const double *m_d_rigidInverse, double *d_rmsArray, 
 55:                         const int *m_d_nRigidBody); 63:                         const int *m_d_nRigidBody);
 56:  64: 
 57:         __global__ void warpReduce1(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  65:         __global__ void ones4(const int *m_d_nRigidBody, double *d_nRigidBodyOnes);
 58:                         const int *m_d_rigidMaxSite, const double *d_gk, double *m_d_gkRigid); 
 59:  
 60:         __global__ void warpReduce2(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  
 61:                         double *m_d_gkRigid, const double *d_tempArray); 
 62:  
 63:         __global__ void warpReduce3(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  
 64:                         double *d_torques, const double *d_tempArray); 
 65:  
 66:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  
 67:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double3 *d_outArray); 
 68:  
 69:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices,  
 70:                         double3 *d_outArray); 
 71:  
 72:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
 73:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double3 *d_outArray, 
 74:                         const double *d_tempArray); 
 75:  
 76:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices,  
 77:                         double3 *d_outArray, const int *m_d_nRigidBody); 
 78:  
 79:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *torques, const int *m_d_largeRigidIndices,  
 80:                         double3 *d_outArray, const int *m_d_nRigidBody); 
 81:  
 82:         __global__ void fullReduce(double *in, double* out, int N); 
 83:  
 84:         __global__ void fullDotReduce(double *in, double* out, int N); 
 85:  
 86: } 66: }
 87:  67: 
 88:  68: 
 89:  69: 
 90: void CostFunction::transformRigidToC(double *d_x) 70: void CostFunction::transformRigidToC(double *d_x)
 91: { 71: {
 92:         using namespace gpu_rigid_bodies; 72:         using namespace gpu_rigid_bodies;
 93:  73: 
 94:         // Copy atomistic coordinates to rigid coordinates.  74:         // Copy atomistic coordinates to rigid coordinates. 
 95:         CudaSafeCall( cudaMemcpy(m_d_xRigid, d_x, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) ); 75:         CudaSafeCall( cudaMemcpy(m_d_xRigid, d_x, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );
166:         CudaSafeCall( cudaMalloc(&d_grmi3, 9 * m_nRigidBody * sizeof(double)) );146:         CudaSafeCall( cudaMalloc(&d_grmi3, 9 * m_nRigidBody * sizeof(double)) );
167: 147: 
168:         bool shouldFindDerivatives = true;148:         bool shouldFindDerivatives = true;
169:         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::shouldFindDeriv, &shouldFindDerivatives,  sizeof(bool)) );149:         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::shouldFindDeriv, &shouldFindDerivatives,  sizeof(bool)) );
170: 150: 
171:         dim3 blockDim;151:         dim3 blockDim;
172:         blockDim.x = 256;152:         blockDim.x = 256;
173:         dim3 gridDim;153:         dim3 gridDim;
174:         gridDim.x = (m_nRigidBody + blockDim.x - 1)/blockDim.x;154:         gridDim.x = (m_nRigidBody + blockDim.x - 1)/blockDim.x;
175: 155: 
176:         // Calculation of rotation matrices - one thread per rigid body. 156:         // One thread per rigid body. 
177:         gpu_rigid_bodies::gradTransform1a<<<gridDim, blockDim>>>(m_d_xRigid, d_grmi1, d_grmi2, d_grmi3, m_d_nRigidBody);157:         gpu_rigid_bodies::gradTransform1a<<<gridDim, blockDim>>>(m_d_xRigid, d_grmi1, d_grmi2, d_grmi3, m_d_nRigidBody);
178:         CudaCheckError();158:         CudaCheckError();
179:         cudaDeviceSynchronize();159:         cudaDeviceSynchronize();
180: 160: 
181:         // Reduction - projection of atomistic forces onto translational degrees of freedom of each rigid body.161:         // Appropriate filling with zeros and ones allows summation of array sections using a dot product. 
182:         blockDim.x = 1024;162:         double *d_xZeroOne;
183:         gridDim.x = (32*m_nRigidBody + blockDim.x - 1)/blockDim.x; // 32 is warpSize163:         double *d_yZeroOne;
184:         // Reduction only takes place for rigid bodies with 32 or fewer sites. 164:         double *d_zZeroOne;
185:         gpu_rigid_bodies::warpReduce1<<<gridDim, blockDim>>>(m_d_nRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups, m_d_rigidMaxSite, 165: 
186:                         d_gk, m_d_gkRigid);166:         CudaSafeCall( cudaMalloc(&d_xZeroOne, numDimensions * sizeof(double)) );
187:         CudaCheckError();167:         CudaSafeCall( cudaMalloc(&d_yZeroOne, numDimensions * sizeof(double)) );
188:         cudaDeviceSynchronize();168:         CudaSafeCall( cudaMalloc(&d_zZeroOne, numDimensions * sizeof(double)) );
 169: 
 170:         // Projection of atomistic forces onto translational degrees of freedom of each rigid body. 
 171:         for (int hostThisBody = 0; hostThisBody < m_nRigidBody; ++hostThisBody) {
 172:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::thisBody, &hostThisBody, sizeof(int)) );
 173:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::arraySize, &numDimensions, sizeof(int)) );
189: 174: 
190:         // Reduction for rigid bodies with more than 32 sites. 175:                 blockDim.x = 1024;
191:         if (m_nLargeRigidBody != 0){176:                 gridDim.x = (numDimensions + blockDim.x - 1)/blockDim.x;
192:                 int blockSize = 1024; 
193:                 blockDim.x = blockSize; 
194:  
195:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies.  
196:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize; 
197:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
198:                 int numThreads = roundedMaxSite*m_nLargeRigidBody; 
199:  
200:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
201:                 gridDim.x = blocks; 
202:  
203:                 double3 *d_outArray; 
204:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) ); 
205:  
206:                 gpu_rigid_bodies::fullReduce1a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody,  
207:                                 m_d_rigidGroups, m_d_rigidMaxSite, d_gk, m_d_largeRigidIndices, d_outArray); 
208:  
209:                 while (blocks > m_nLargeRigidBody) { 
210:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) ); 
211:  
212:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize; 
213:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
214:                         numThreads = roundedMaxSite*m_nLargeRigidBody; 
215:  
216:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
217:                         gridDim.x = blocks; 
218:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) ); 
219: 177: 
220:                         gpu_rigid_bodies::fullReduce2a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, 178:                 // Fill array with zeros. 
221:                                         m_d_largeRigidIndices, d_outArray);179:                 gpu_rigid_bodies::zeros1<<<gridDim, blockDim>>>(d_xZeroOne, d_yZeroOne, d_zZeroOne);
222:                 }180:                 CudaCheckError();
 181:                 cudaDeviceSynchronize();
 182: 
 183:                 // Change zeros to ones up to appropriate number of rigid sites. 
 184:                 gpu_rigid_bodies::ones1<<<gridDim, blockDim>>>(m_d_nRigidSitesPerBody, m_d_rigidGroups, d_xZeroOne, 
 185:                                 d_yZeroOne, d_zZeroOne, m_d_rigidMaxSite);
 186:                 CudaCheckError();
 187:                 cudaDeviceSynchronize();
 188: 
 189:                 // Reduction - library dot product is fastest way for subsections of arrays. 
 190:                 m_cublas.dispatchDot(numDimensions, (m_d_gkRigid + 3*hostThisBody), d_xZeroOne, d_gk); // gkRigid = xZeroOne Dot gk
 191:                 m_cublas.dispatchDot(numDimensions, (m_d_gkRigid + 3*hostThisBody + 1), d_yZeroOne, d_gk); // gkRigid = yZeroOne Dot gk
 192:                 m_cublas.dispatchDot(numDimensions, (m_d_gkRigid + 3*hostThisBody + 2), d_zZeroOne, d_gk); // gkRigid = zZeroOne Dot gk
223: 193: 
224:                 CudaSafeCall( cudaFree(d_outArray) ); 
225:         }194:         }
226: 195: 
 196:         CudaSafeCall( cudaFree(d_xZeroOne) );
 197:         CudaSafeCall( cudaFree(d_yZeroOne) );
 198:         CudaSafeCall( cudaFree(d_zZeroOne) );
 199: 
227:         // Temporary global memory storage for atomistic components of rigid body forces. 200:         // Temporary global memory storage for atomistic components of rigid body forces. 
228:         double *d_tempArray;201:         double *d_tempArray;
229:         int tempArraySize = 3 * m_nRigidBody * m_rigidMaxSite;202:         int tempArraySize = 3 * m_nRigidBody * m_rigidMaxSite;
230:         CudaSafeCall( cudaMalloc(&d_tempArray, tempArraySize * sizeof(double)) );203:         CudaSafeCall( cudaMalloc(&d_tempArray, tempArraySize * sizeof(double)) );
231: 204: 
232:         // Launch no. of threads equal to max. no. of sites per rigid body * number of rigid bodies for easy indexing. 205:         // Launch no. of threads equal to max. no. of sites per rigid body * number of rigid bodies for easy indexing. 
233:         // Not all threads wil be used as number of sites per rigid body can be less than maximum. 206:         // Not all threads wil be used as number of sites per rigid body can be less than maximum. 
234:         int noThreads = m_rigidMaxSite * m_nRigidBody;207:         int noThreads = m_rigidMaxSite * m_nRigidBody;
235:         blockDim.x = 1024;208:         blockDim.x = 1024;
236:         gridDim.x = (noThreads + blockDim.x - 1)/blockDim.x;209:         gridDim.x = (noThreads + blockDim.x - 1)/blockDim.x;
239:         gpu_rigid_bodies::intermediate<<<gridDim, blockDim>>>(d_gk, m_d_nRigidSitesPerBody, d_grmi1, d_grmi2, d_grmi3, 212:         gpu_rigid_bodies::intermediate<<<gridDim, blockDim>>>(d_gk, m_d_nRigidSitesPerBody, d_grmi1, d_grmi2, d_grmi3, 
240:                         m_d_sitesRigidBody, m_d_rigidGroups, d_tempArray, 213:                         m_d_sitesRigidBody, m_d_rigidGroups, d_tempArray, 
241:                         m_d_nRigidBody, m_d_rigidMaxSite);214:                         m_d_nRigidBody, m_d_rigidMaxSite);
242:         CudaCheckError();215:         CudaCheckError();
243:         cudaDeviceSynchronize();216:         cudaDeviceSynchronize();
244: 217: 
245:         CudaSafeCall( cudaFree(d_grmi1) );218:         CudaSafeCall( cudaFree(d_grmi1) );
246:         CudaSafeCall( cudaFree(d_grmi2) );219:         CudaSafeCall( cudaFree(d_grmi2) );
247:         CudaSafeCall( cudaFree(d_grmi3) );220:         CudaSafeCall( cudaFree(d_grmi3) );
248: 221: 
249:         // Reduction of atomistic components (d_tempArray) gives forces projected onto rotational d.o.f. of each rigid body.222:         // Appropriate filling with zeros and ones allows summation of array sections using a dot product. 
250:         blockDim.x = 1024;223:         double *d_xZeroOne2;
251:         gridDim.x = (32*m_nRigidBody + blockDim.x - 1)/blockDim.x; // 32 is warpSize224:         double *d_yZeroOne2;
252:         // Reduction only takes place for rigid bodies with 32 or fewer sites.225:         double *d_zZeroOne2;
253:         gpu_rigid_bodies::warpReduce2<<<gridDim, blockDim>>>(m_d_nRigidBody, m_d_nRigidSitesPerBody, m_d_rigidMaxSite, 226: 
254:                         m_d_gkRigid, d_tempArray);227:         CudaSafeCall( cudaMalloc(&d_xZeroOne2, tempArraySize * sizeof(double)) );
255:         CudaCheckError();228:         CudaSafeCall( cudaMalloc(&d_yZeroOne2, tempArraySize * sizeof(double)) );
256:         cudaDeviceSynchronize();229:         CudaSafeCall( cudaMalloc(&d_zZeroOne2, tempArraySize * sizeof(double)) );
 230: 
 231:         // Reduction of atomistic components (d_tempArray) gives forces projected onto rotational d.o.f. of each rigid body. 
 232:         for (int hostThisBody = 0; hostThisBody < m_nRigidBody; ++hostThisBody) {
 233:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::thisBody, &hostThisBody, sizeof(int)) );
 234:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::arraySize, &tempArraySize, sizeof(int)) );
257: 235: 
258:         // Reduction for rigid bodies with more than 32 sites.236:                 blockDim.x = 1024;
259:         if (m_nLargeRigidBody != 0){237:                 gridDim.x = (tempArraySize + blockDim.x - 1)/blockDim.x;
260:                 int blockSize = 1024; 
261:                 blockDim.x = blockSize; 
262:  
263:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies. 
264:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize; 
265:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
266:                 int numThreads = roundedMaxSite*m_nLargeRigidBody; 
267:  
268:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
269:                 gridDim.x = blocks; 
270:  
271:                 double3 *d_outArray; 
272:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) ); 
273:  
274:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups, 
275:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray); 
276:  
277:                 while (blocks > m_nLargeRigidBody) { 
278:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) ); 
279:  
280:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize; 
281:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
282:                         numThreads = roundedMaxSite*m_nLargeRigidBody; 
283:  
284:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
285:                         gridDim.x = blocks; 
286:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) ); 
287: 238: 
288:                         gpu_rigid_bodies::fullReduce2b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, m_d_largeRigidIndices, 239:                 // Fill array with zeros. 
289:                                         d_outArray, m_d_nRigidBody);240:                 gpu_rigid_bodies::zeros2<<<gridDim, blockDim>>>(d_xZeroOne2, d_yZeroOne2, d_zZeroOne2);
290:                 }241:                 CudaCheckError();
 242:                 cudaDeviceSynchronize();
 243: 
 244:                 // Change zeros to ones up to appropriate number of rigid sites. 
 245:                 gpu_rigid_bodies::ones2<<<gridDim, blockDim>>>(m_d_nRigidSitesPerBody, d_xZeroOne2, d_yZeroOne2, 
 246:                                 d_zZeroOne2, m_d_rigidMaxSite);
 247:                 CudaCheckError();
 248:                 cudaDeviceSynchronize();
291: 249: 
292:                 CudaSafeCall( cudaFree(d_outArray) );250:                 // Reduction - library dot product is fastest way for subsections of arrays. 
 251:                 m_cublas.dispatchDot(tempArraySize, (m_d_gkRigid + 3*m_nRigidBody + 3*hostThisBody),     d_xZeroOne2, d_tempArray); // gkRigid = xZeroOne2 Dot tempArray
 252:                 m_cublas.dispatchDot(tempArraySize, (m_d_gkRigid + 3*m_nRigidBody + 3*hostThisBody + 1), d_yZeroOne2, d_tempArray); // gkRigid = yZeroOne2 Dot tempArray
 253:                 m_cublas.dispatchDot(tempArraySize, (m_d_gkRigid + 3*m_nRigidBody + 3*hostThisBody + 2), d_zZeroOne2, d_tempArray); // gkRigid = zZeroOne2 Dot tempArray
293:         }254:         }
294: 255: 
 256:         CudaSafeCall( cudaFree(d_xZeroOne2) );
 257:         CudaSafeCall( cudaFree(d_yZeroOne2) );
 258:         CudaSafeCall( cudaFree(d_zZeroOne2) );
295:         CudaSafeCall( cudaFree(d_tempArray) );259:         CudaSafeCall( cudaFree(d_tempArray) );
296: 260: 
297:         if (m_nDegFreedom > 6*m_nRigidBody) {261:         if (m_nDegFreedom > 6*m_nRigidBody) {
298:                 int hostSinglesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;262:                 int hostSinglesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;
299:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::singlesThreads, &hostSinglesThreads, sizeof(int)) );263:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::singlesThreads, &hostSinglesThreads, sizeof(int)) );
300: 264: 
301:                 blockDim.x = 1024;265:                 blockDim.x = 1024;
302:                 gridDim.x = (hostSinglesThreads + blockDim.x - 1)/blockDim.x;266:                 gridDim.x = (hostSinglesThreads + blockDim.x - 1)/blockDim.x;
303: 267: 
304:                 // Treatment of single atoms not in rigid bodies. 268:                 // Treatment of single atoms not in rigid bodies. 
309:         }273:         }
310: 274: 
311:         // Copy rigid coordinates/gradient to atomistic coordinates/gradient and pad with zeros. 275:         // Copy rigid coordinates/gradient to atomistic coordinates/gradient and pad with zeros. 
312:         CudaSafeCall( cudaMemcpy(d_x,  zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );276:         CudaSafeCall( cudaMemcpy(d_x,  zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
313:         CudaSafeCall( cudaMemcpy(d_gk, zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );277:         CudaSafeCall( cudaMemcpy(d_gk, zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
314: 278: 
315:         CudaSafeCall( cudaMemcpy(d_x,  m_d_xRigid,  m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );279:         CudaSafeCall( cudaMemcpy(d_x,  m_d_xRigid,  m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );
316:         CudaSafeCall( cudaMemcpy(d_gk, m_d_gkRigid, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );280:         CudaSafeCall( cudaMemcpy(d_gk, m_d_gkRigid, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );
317: 281: 
318:         delete [] zeros;282:         delete [] zeros;
319:  
320: }283: }
321: 284: 
322: 285: 
323: 286: 
324: void CostFunction::aaConvergence(const double *d_gk, double *outRms)287: void CostFunction::aaConvergence(const double *d_gk, double *outRms)
325: {288: {
326:         using namespace gpu_rigid_bodies;289:         using namespace gpu_rigid_bodies;
327: 290: 
328:         *outRms = 0.0;291:         *outRms = 0.0;
329: 292: 
369:         gpu_rigid_bodies::intermediate2<<<gridDim, blockDim>>>(d_gk, m_d_nRigidSitesPerBody, d_grmi0, d_rmi10, d_rmi20, 332:         gpu_rigid_bodies::intermediate2<<<gridDim, blockDim>>>(d_gk, m_d_nRigidSitesPerBody, d_grmi0, d_rmi10, d_rmi20, 
370:                         d_rmi30, m_d_sitesRigidBody, m_d_rigidGroups, d_tempArray, 333:                         d_rmi30, m_d_sitesRigidBody, m_d_rigidGroups, d_tempArray, 
371:                         m_d_nRigidBody, m_d_rigidMaxSite);334:                         m_d_nRigidBody, m_d_rigidMaxSite);
372:         CudaCheckError();335:         CudaCheckError();
373:         cudaDeviceSynchronize();336:         cudaDeviceSynchronize();
374: 337: 
375:         CudaSafeCall( cudaFree(d_rmi10) );338:         CudaSafeCall( cudaFree(d_rmi10) );
376:         CudaSafeCall( cudaFree(d_rmi20) );339:         CudaSafeCall( cudaFree(d_rmi20) );
377:         CudaSafeCall( cudaFree(d_rmi30) );340:         CudaSafeCall( cudaFree(d_rmi30) );
378: 341: 
 342:         double *d_xZeroOne2;
 343:         double *d_yZeroOne2;
 344:         double *d_zZeroOne2;
379:         double *d_torques;345:         double *d_torques;
 346: 
 347:         CudaSafeCall( cudaMalloc(&d_xZeroOne2,  tempArraySize * sizeof(double)) );
 348:         CudaSafeCall( cudaMalloc(&d_yZeroOne2,  tempArraySize * sizeof(double)) );
 349:         CudaSafeCall( cudaMalloc(&d_zZeroOne2,  tempArraySize * sizeof(double)) );
380:         CudaSafeCall( cudaMalloc(&d_torques, 3 * m_nRigidBody * sizeof(double)) );350:         CudaSafeCall( cudaMalloc(&d_torques, 3 * m_nRigidBody * sizeof(double)) );
381: 351: 
382:         // Reduction. 352:         for (int hostThisBody = 0; hostThisBody < m_nRigidBody; ++hostThisBody) {
383:         blockDim.x = 1024;353:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::thisBody,   &hostThisBody, sizeof(int)) );
384:         gridDim.x = (32*m_nRigidBody + blockDim.x - 1)/blockDim.x; // 32 is warpSize354:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::arraySize, &tempArraySize, sizeof(int)) );
385:         // Reduction only takes place for rigid bodies with 32 or fewer sites. 
386:         gpu_rigid_bodies::warpReduce3<<<gridDim, blockDim>>>(m_d_nRigidBody, m_d_nRigidSitesPerBody, m_d_rigidMaxSite, d_torques,  
387:                         d_tempArray); 
388:         CudaCheckError(); 
389:         cudaDeviceSynchronize(); 
390: 355: 
391:         // Reduction for rigid bodies with more than 32 sites.356:                 blockDim.x = 1024;
392:         if (m_nLargeRigidBody != 0){357:                 gridDim.x = (tempArraySize + blockDim.x - 1)/blockDim.x;
393:                 int blockSize = 1024; 
394:                 blockDim.x = blockSize; 
395:  
396:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies. 
397:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize; 
398:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
399:                 int numThreads = roundedMaxSite*m_nLargeRigidBody; 
400:  
401:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
402:                 gridDim.x = blocks; 
403:  
404:                 double3 *d_outArray; 
405:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) ); 
406:  
407:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups, 
408:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray); 
409:  
410:                 while (blocks > m_nLargeRigidBody) { 
411:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) ); 
412:  
413:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize; 
414:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) ); 
415:                         numThreads = roundedMaxSite*m_nLargeRigidBody; 
416:  
417:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
418:                         gridDim.x = blocks; 
419:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) ); 
420: 358: 
421:                         gpu_rigid_bodies::fullReduce2c<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, d_torques, 359:                 gpu_rigid_bodies::zeros2<<<gridDim, blockDim>>>(d_xZeroOne2, d_yZeroOne2, d_zZeroOne2);
422:                                         m_d_largeRigidIndices, d_outArray, m_d_nRigidBody);360:                 CudaCheckError();
423:                 }361:                 cudaDeviceSynchronize();
424: 362: 
425:                 CudaSafeCall( cudaFree(d_outArray) );363:                 gpu_rigid_bodies::ones3<<<gridDim, blockDim>>>(d_xZeroOne2, d_yZeroOne2, d_zZeroOne2, m_d_rigidMaxSite);
 364:                 CudaCheckError();
 365:                 cudaDeviceSynchronize();
 366: 
 367:                 m_cublas.dispatchDot(tempArraySize, d_torques + 3*hostThisBody,     d_xZeroOne2, d_tempArray); // torques =  xZeroOne2 Dot tempArray
 368:                 m_cublas.dispatchDot(tempArraySize, d_torques + 3*hostThisBody + 1, d_yZeroOne2, d_tempArray); // torques =  yZeroOne2 Dot tempArray
 369:                 m_cublas.dispatchDot(tempArraySize, d_torques + 3*hostThisBody + 2, d_zZeroOne2, d_tempArray); // torques =  zZeroOne2 Dot twmpArray
426:         }370:         }
427: 371: 
428:         CudaSafeCall( cudaFree(d_tempArray) );372:         CudaSafeCall( cudaFree(d_tempArray) );
 373:         CudaSafeCall( cudaFree(d_xZeroOne2) );
 374:         CudaSafeCall( cudaFree(d_yZeroOne2) );
 375:         CudaSafeCall( cudaFree(d_zZeroOne2) );
429: 376: 
430:         double *d_rmsArray;377:         double *d_rmsArray;
431:         CudaSafeCall( cudaMalloc(&d_rmsArray, m_nRigidBody * sizeof(double)) );378:         double *d_nRigidBodyOnes;
 379: 
 380:         CudaSafeCall( cudaMalloc(&d_rmsArray,       m_nRigidBody * sizeof(double)) );
 381:         CudaSafeCall( cudaMalloc(&d_nRigidBodyOnes, m_nRigidBody * sizeof(double)) );
432: 382: 
433:         blockDim.x = 1024;383:         blockDim.x = 1024;
434:         gridDim.x = (m_nRigidBody + blockDim.x - 1)/blockDim.x;384:         gridDim.x = (m_nRigidBody + blockDim.x - 1)/blockDim.x;
435: 385: 
436:         gpu_rigid_bodies::aaConvTorque<<<gridDim, blockDim>>>(d_torques, m_d_gkRigid, d_grmi0, m_d_nRigidSitesPerBody, 386:         gpu_rigid_bodies::aaConvTorque<<<gridDim, blockDim>>>(d_torques, m_d_gkRigid, d_grmi0, m_d_nRigidSitesPerBody, 
437:                         m_d_rigidInverse, d_rmsArray, m_d_nRigidBody);387:                         m_d_rigidInverse, d_rmsArray, m_d_nRigidBody);
438:         CudaCheckError();388:         CudaCheckError();
439:         cudaDeviceSynchronize();389:         cudaDeviceSynchronize();
440: 390: 
441:         CudaSafeCall( cudaFree(d_grmi0) );391:         CudaSafeCall( cudaFree(d_grmi0) );
442:         CudaSafeCall( cudaFree(d_torques) );392:         CudaSafeCall( cudaFree(d_torques) );
443: 393: 
444:         double temp1 = 0.0;394:         gpu_rigid_bodies::ones4<<<gridDim, blockDim>>>(m_d_nRigidBody, d_nRigidBodyOnes);
445:         double temp2 = 0.0; 
446:  
447:         // Reduction.  
448:         blockDim.x = 512; 
449:         gridDim.x = min((m_nRigidBody + blockDim.x - 1) / blockDim.x, 1024); 
450:  
451:         double *d_outArray; 
452:         CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) ); 
453:  
454:         gpu_rigid_bodies::fullReduce<<<gridDim, blockDim>>>(d_rmsArray, d_outArray, m_nRigidBody); 
455:         gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x); 
456:         CudaCheckError();395:         CudaCheckError();
457:         cudaDeviceSynchronize();396:         cudaDeviceSynchronize();
458: 397: 
459:         CudaSafeCall( cudaMemcpy(&temp1, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) );398:         double temp1 = 0.0;
 399:         double temp2 = 0.0;
 400: 
 401:         m_cublas.dispatchDot(m_nRigidBody, &temp1, d_rmsArray, d_nRigidBodyOnes, false); // temp1 =  rmsArray Dot nRigidBodyOnes
460: 402: 
461:         CudaSafeCall( cudaFree(d_outArray) ); 
462:         CudaSafeCall( cudaFree(d_rmsArray) );403:         CudaSafeCall( cudaFree(d_rmsArray) );
 404:         CudaSafeCall( cudaFree(d_nRigidBodyOnes) );
463: 405: 
464:         if (m_nDegFreedom > 6*m_nRigidBody) {406:         if (m_nDegFreedom > 6*m_nRigidBody) {
465:                 int arraysize = m_nDegFreedom - 6*m_nRigidBody;407:                 int arraysize = m_nDegFreedom - 6*m_nRigidBody;
466: 408: 
467:                 // Dot product.  409:                 m_cublas.dispatchDot(arraysize, &temp2, (m_d_gkRigid + 6*m_nRigidBody), (m_d_gkRigid + 6*m_nRigidBody), false); // temp2 =  gkRigid Dot gkRigid
468:                 blockDim.x = 512; 
469:                 gridDim.x = min((arraysize + blockDim.x - 1) / blockDim.x, 1024); 
470:  
471:                 CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) ); 
472:  
473:                 gpu_rigid_bodies::fullDotReduce<<<gridDim, blockDim>>>((m_d_gkRigid + 6*m_nRigidBody), d_outArray, arraysize); 
474:                 gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x); 
475:                 CudaCheckError(); 
476:                 cudaDeviceSynchronize(); 
477:  
478:                 CudaSafeCall( cudaMemcpy(&temp2, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) ); 
479:  
480:                 CudaSafeCall( cudaFree(d_outArray) ); 
481:         }410:         }
482: 411: 
483:         *outRms = temp1 + temp2;412:         *outRms = temp1 + temp2;
484: }413: }
485: 414: 
486: 415: 
487: 416: 
488: namespace gpu_rigid_bodies417: namespace gpu_rigid_bodies
489: {418: {
490:         __device__ void rmdrvt(double3 p, double rmi[9], double drmi1[9], double drmi2[9], double drmi3[9], 419:         __device__ void rmdrvt(double3 p, double rmi[9], double drmi1[9], double drmi2[9], double drmi3[9], 
896:                         d_grmi0[5+9*tid]= drmi0[5];825:                         d_grmi0[5+9*tid]= drmi0[5];
897: 826: 
898:                         d_grmi0[6+9*tid]= drmi0[6];827:                         d_grmi0[6+9*tid]= drmi0[6];
899:                         d_grmi0[7+9*tid]= drmi0[7];828:                         d_grmi0[7+9*tid]= drmi0[7];
900:                         d_grmi0[8+9*tid]= drmi0[8];829:                         d_grmi0[8+9*tid]= drmi0[8];
901: 830: 
902:                         tid += blockDim.x * gridDim.x;831:                         tid += blockDim.x * gridDim.x;
903:                 }832:                 }
904:         }833:         }
905: 834: 
 835:         __global__ void zeros1(double *d_xZeroOne, double *d_yZeroOne, double *d_zZeroOne)
 836:         {
 837:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 838: 
 839:                 while (tid < arraySize) {
 840:                         d_xZeroOne[tid] = 0.0;
 841:                         d_yZeroOne[tid] = 0.0;
 842:                         d_zZeroOne[tid] = 0.0;
 843: 
 844:                         tid += blockDim.x * gridDim.x;
 845:                 }
 846:         }
 847: 
 848:         __global__ void ones1(const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, double *d_xZeroOne, 
 849:                         double *d_yZeroOne, double *d_zZeroOne, const int *m_d_rigidMaxSite)
 850:         {
 851:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 852: 
 853:                 while (tid < m_d_nRigidSitesPerBody[thisBody]) {
 854:                         int myAtom = m_d_rigidGroups[tid+(*m_d_rigidMaxSite)*thisBody];
 855: 
 856:                         d_xZeroOne[3*myAtom-3] = 1.0;
 857:                         d_yZeroOne[3*myAtom-2] = 1.0;
 858:                         d_zZeroOne[3*myAtom-1] = 1.0;
 859: 
 860:                         tid += blockDim.x * gridDim.x;
 861:                 }
 862:         }
 863: 
906:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1, 864:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1, 
907:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody, 865:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody, 
908:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody, 866:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody, 
909:                         const int *m_d_rigidMaxSite)867:                         const int *m_d_rigidMaxSite)
910:         {868:         {
911:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;869:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
912: 870: 
913:                 int thisRigidBody = tid / (*m_d_rigidMaxSite);871:                 int thisRigidBody = tid / (*m_d_rigidMaxSite);
914: 872: 
915:                 while (tid < (*m_d_nRigidBody) * (*m_d_rigidMaxSite)) {873:                 while (tid < (*m_d_nRigidBody) * (*m_d_rigidMaxSite)) {
997:                                 d_tempArray[2+3*i+(*m_d_rigidMaxSite)*3*thisRigidBody] = d_gk[3*myAtom-3]*dr3.x + 955:                                 d_tempArray[2+3*i+(*m_d_rigidMaxSite)*3*thisRigidBody] = d_gk[3*myAtom-3]*dr3.x + 
998:                                         d_gk[3*myAtom-2]*dr3.y + 956:                                         d_gk[3*myAtom-2]*dr3.y + 
999:                                         d_gk[3*myAtom-1]*dr3.z;957:                                         d_gk[3*myAtom-1]*dr3.z;
1000: 958: 
1001:                         }959:                         }
1002: 960: 
1003:                         tid += blockDim.x * gridDim.x;961:                         tid += blockDim.x * gridDim.x;
1004:                 }962:                 }
1005:         }963:         }
1006: 964: 
 965:         __global__ void zeros2(double *d_xZeroOne2, double *d_yZeroOne2, double *d_zZeroOne2)
 966:         {
 967:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 968: 
 969:                 while (tid < arraySize) {
 970:                         d_xZeroOne2[tid] = 0.0;
 971:                         d_yZeroOne2[tid] = 0.0;
 972:                         d_zZeroOne2[tid] = 0.0;
 973: 
 974:                         tid += blockDim.x * gridDim.x;
 975:                 }
 976:         }
 977: 
 978:         __global__ void ones2(const int *m_d_nRigidSitesPerBody, double *d_xZeroOne2, double *d_yZeroOne2, 
 979:                         double *d_zZeroOne2, const int *m_d_rigidMaxSite)
 980:         {
 981:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 982: 
 983:                 while (tid < m_d_nRigidSitesPerBody[thisBody]) {
 984:                         d_xZeroOne2[0+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 985:                         d_yZeroOne2[1+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 986:                         d_zZeroOne2[2+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 987: 
 988:                         tid += blockDim.x * gridDim.x;
 989:                 }
 990:         }
 991: 
 992:         __global__ void ones3(double *d_xZeroOne2, double *d_yZeroOne2, double *d_zZeroOne2, const int *m_d_rigidMaxSite)
 993:         {
 994:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 995: 
 996:                 while (tid < (*m_d_rigidMaxSite)) {
 997:                         d_xZeroOne2[0+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 998:                         d_yZeroOne2[1+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 999:                         d_zZeroOne2[2+3*tid+3*(*m_d_rigidMaxSite)*thisBody] = 1.0;
 1000: 
 1001:                         tid += blockDim.x * gridDim.x;
 1002:                 }
 1003:         }
 1004: 
 1005:         __global__ void ones4(const int *m_d_nRigidBody, double *d_nRigidBodyOnes)
 1006:         {
 1007:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
 1008: 
 1009:                 while (tid < (*m_d_nRigidBody)) {
 1010:                         d_nRigidBodyOnes[tid] = 1.0;
 1011: 
 1012:                         tid += blockDim.x * gridDim.x;
 1013:                 }
 1014:         }
 1015: 
1007:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 1016:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 
1008:                         const int *m_d_nRigidBody)1017:                         const int *m_d_nRigidBody)
1009:         {1018:         {
1010:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1019:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1011: 1020: 
1012:                 while (tid < singlesThreads) {1021:                 while (tid < singlesThreads) {
1013:                         int thisAtom = m_d_rigidSingles[tid];1022:                         int thisAtom = m_d_rigidSingles[tid];
1014: 1023: 
1015:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+0] = d_gk[3*thisAtom-3];1024:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+0] = d_gk[3*thisAtom-3];
1016:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+1] = d_gk[3*thisAtom-2];1025:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+1] = d_gk[3*thisAtom-2];
1179:                         double rmscontrib2 = (1.0/m_d_nRigidSitesPerBody[tid])*(m_d_gkRigid[3*tid+0]*m_d_gkRigid[3*tid+0] + 1188:                         double rmscontrib2 = (1.0/m_d_nRigidSitesPerBody[tid])*(m_d_gkRigid[3*tid+0]*m_d_gkRigid[3*tid+0] + 
1180:                                         m_d_gkRigid[3*tid+1]*m_d_gkRigid[3*tid+1] + 1189:                                         m_d_gkRigid[3*tid+1]*m_d_gkRigid[3*tid+1] + 
1181:                                         m_d_gkRigid[3*tid+2]*m_d_gkRigid[3*tid+2]);1190:                                         m_d_gkRigid[3*tid+2]*m_d_gkRigid[3*tid+2]);
1182: 1191: 
1183: 1192: 
1184:                         d_rmsArray[tid] = rmscontrib1 + rmscontrib2;1193:                         d_rmsArray[tid] = rmscontrib1 + rmscontrib2;
1185: 1194: 
1186:                         tid += blockDim.x * gridDim.x;1195:                         tid += blockDim.x * gridDim.x;
1187:                 }1196:                 }
1188:         }1197:         }
1189:  
1190:         // Undocumented CUDA 6.5 function copied and pasted here as not present in CUDA 6.0 and below.  
1191:         static __device__ __inline__  
1192:                 double __myshfl_down(double var, unsigned int delta, int width=warpSize) { 
1193:                         float lo, hi; 
1194:                         asm volatile("mov.b64 {%0,%1}, %2;" : "=f"(lo), "=f"(hi) : "d"(var)); 
1195:                         hi = __shfl_down(hi, delta, width); 
1196:                         lo = __shfl_down(lo, delta, width); 
1197:                         asm volatile("mov.b64 %0, {%1,%2};" : "=d"(var) : "f"(lo), "f"(hi)); 
1198:                         return var; 
1199:                 } 
1200:  
1201:         // See 'Faster Parallel Reductions on Kepler' from the NVIDIA dev blog.  
1202:         __inline__ __device__ 
1203:                 double warpReduceSum(double val) { 
1204:                         for (int offset = warpSize/2; offset > 0; offset /= 2) { 
1205:                                 val += __myshfl_down(val, offset); 
1206:                         } 
1207:                         return val; 
1208:                 } 
1209:  
1210:         __global__ void warpReduce1(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  
1211:                         const int *m_d_rigidMaxSite, const double *d_gk, double *m_d_gkRigid) 
1212:         { 
1213:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1214:  
1215:                 while (tid < (warpSize*(*m_d_nRigidBody))) { 
1216:  
1217:                         double3 elements; 
1218:                         elements.x = 0.0; 
1219:                         elements.y = 0.0; 
1220:                         elements.z = 0.0; 
1221:  
1222:                         int index = ((tid + warpSize) % warpSize); 
1223:                         int currentBody = tid/warpSize; 
1224:  
1225:                         if ((index < m_d_nRigidSitesPerBody[currentBody]) && (m_d_nRigidSitesPerBody[currentBody] <= warpSize)) { 
1226:                                 int myAtom = m_d_rigidGroups[index+(*m_d_rigidMaxSite)*currentBody]; 
1227:  
1228:                                 elements.x = d_gk[3*myAtom-3]; 
1229:                                 elements.y = d_gk[3*myAtom-2]; 
1230:                                 elements.z = d_gk[3*myAtom-1]; 
1231:                         } 
1232:  
1233:                         elements.x = warpReduceSum(elements.x); 
1234:                         elements.y = warpReduceSum(elements.y); 
1235:                         elements.z = warpReduceSum(elements.z); 
1236:  
1237:                         if (tid % warpSize == 0) { 
1238:                                 m_d_gkRigid[3*currentBody] = elements.x; 
1239:                                 m_d_gkRigid[3*currentBody+1] = elements.y; 
1240:                                 m_d_gkRigid[3*currentBody+2] = elements.z; 
1241:                         } 
1242:  
1243:                         tid += blockDim.x * gridDim.x; 
1244:                 } 
1245:         } 
1246:  
1247:         __global__ void warpReduce2(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  
1248:                         double *m_d_gkRigid, const double *d_tempArray) 
1249:         { 
1250:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1251:  
1252:                 while (tid < (warpSize*(*m_d_nRigidBody))) { 
1253:  
1254:                         double3 elements; 
1255:                         elements.x = 0.0; 
1256:                         elements.y = 0.0; 
1257:                         elements.z = 0.0; 
1258:  
1259:                         int index = ((tid + warpSize) % warpSize); 
1260:                         int currentBody = tid/warpSize; 
1261:  
1262:                         if ((index < m_d_nRigidSitesPerBody[currentBody]) && (m_d_nRigidSitesPerBody[currentBody] <= warpSize)) { 
1263:                                 elements.x = d_tempArray[0+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1264:                                 elements.y = d_tempArray[1+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1265:                                 elements.z = d_tempArray[2+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1266:                         } 
1267:  
1268:                         elements.x = warpReduceSum(elements.x); 
1269:                         elements.y = warpReduceSum(elements.y); 
1270:                         elements.z = warpReduceSum(elements.z); 
1271:  
1272:                         if (tid % warpSize == 0) { 
1273:                                 m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody] = elements.x; 
1274:                                 m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+1] = elements.y; 
1275:                                 m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+2] = elements.z; 
1276:                         } 
1277:  
1278:                         tid += blockDim.x * gridDim.x; 
1279:                 } 
1280:         } 
1281:  
1282:         __global__ void warpReduce3(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  
1283:                         double *d_torques, const double *d_tempArray) 
1284:         { 
1285:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1286:  
1287:                 while (tid < (warpSize*(*m_d_nRigidBody))) { 
1288:  
1289:                         double3 elements; 
1290:                         elements.x = 0.0; 
1291:                         elements.y = 0.0; 
1292:                         elements.z = 0.0; 
1293:  
1294:                         int index = ((tid + warpSize) % warpSize); 
1295:                         int currentBody = tid/warpSize; 
1296:  
1297:                         if ((index < m_d_nRigidSitesPerBody[currentBody]) && (m_d_nRigidSitesPerBody[currentBody] <= warpSize)) { 
1298:                                 elements.x = d_tempArray[0+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1299:                                 elements.y = d_tempArray[1+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1300:                                 elements.z = d_tempArray[2+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1301:                         } 
1302:  
1303:                         elements.x = warpReduceSum(elements.x); 
1304:                         elements.y = warpReduceSum(elements.y); 
1305:                         elements.z = warpReduceSum(elements.z); 
1306:  
1307:                         if (tid % warpSize == 0) { 
1308:                                 d_torques[3*currentBody] = elements.x; 
1309:                                 d_torques[3*currentBody+1] = elements.y; 
1310:                                 d_torques[3*currentBody+2] = elements.z; 
1311:                         } 
1312:  
1313:                         tid += blockDim.x * gridDim.x; 
1314:                 } 
1315:         } 
1316:  
1317:         __inline__ __device__ 
1318:                 double blockReduceSum(double val) { 
1319:                         static __shared__ double shared[32]; // Shared mem for 32 partial sums 
1320:                         int lane = threadIdx.x % warpSize; 
1321:                         int wid = threadIdx.x / warpSize; 
1322:  
1323:                         val = warpReduceSum(val);     // Each warp performs partial reduction 
1324:  
1325:                         // Write reduced value to shared memory 
1326:                         if (lane==0) { 
1327:                                 shared[wid]=val; 
1328:                         }  
1329:  
1330:                         __syncthreads();              // Wait for all partial reductions 
1331:  
1332:                         //read from shared memory only if that warp existed 
1333:                         val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0; 
1334:  
1335:                         //Final reduce within first warp 
1336:                         if (wid==0) { 
1337:                                 val = warpReduceSum(val); 
1338:                         }  
1339:  
1340:                         return val; 
1341:                 } 
1342:  
1343:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  
1344:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double3 *d_outArray)  
1345:         { 
1346:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1347:  
1348:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) { 
1349:  
1350:                         double3 elements; 
1351:                         elements.x = 0.0; 
1352:                         elements.y = 0.0; 
1353:                         elements.z = 0.0; 
1354:  
1355:                         int index = ((tid + roundMaxSite) % roundMaxSite); 
1356:                         int thisBody = tid/roundMaxSite; 
1357:                         int currentBody = m_d_largeRigidIndices[thisBody]; 
1358:  
1359:                         if (index < m_d_nRigidSitesPerBody[currentBody]) { 
1360:                                 int myAtom = m_d_rigidGroups[index+(*m_d_rigidMaxSite)*currentBody]; 
1361:                                 elements.x = d_gk[3*myAtom-3]; 
1362:                                 elements.y = d_gk[3*myAtom-2]; 
1363:                                 elements.z = d_gk[3*myAtom-1]; 
1364:                         } 
1365:  
1366:                         elements.x = blockReduceSum(elements.x); 
1367:                         elements.y = blockReduceSum(elements.y); 
1368:                         elements.z = blockReduceSum(elements.z); 
1369:  
1370:                         if (threadIdx.x==0) { 
1371:                                 d_outArray[blockIdx.x].x = elements.x; 
1372:                                 d_outArray[blockIdx.x].y = elements.y; 
1373:                                 d_outArray[blockIdx.x].z = elements.z; 
1374:                         } 
1375:  
1376:                         tid += blockDim.x * gridDim.x; 
1377:                 } 
1378:         } 
1379:  
1380:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  
1381:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double3 *d_outArray,  
1382:                         const double *d_tempArray)  
1383:         { 
1384:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1385:  
1386:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) { 
1387:  
1388:                         double3 elements; 
1389:                         elements.x = 0.0; 
1390:                         elements.y = 0.0; 
1391:                         elements.z = 0.0; 
1392:  
1393:                         int index = ((tid + roundMaxSite) % roundMaxSite); 
1394:                         int thisBody = tid/roundMaxSite; 
1395:                         int currentBody = m_d_largeRigidIndices[thisBody]; 
1396:  
1397:                         if (index < m_d_nRigidSitesPerBody[currentBody]) { 
1398:                                 elements.x = d_tempArray[0+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1399:                                 elements.y = d_tempArray[1+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1400:                                 elements.z = d_tempArray[2+3*index+3*(*m_d_rigidMaxSite)*currentBody]; 
1401:                         } 
1402:  
1403:                         elements.x = blockReduceSum(elements.x); 
1404:                         elements.y = blockReduceSum(elements.y); 
1405:                         elements.z = blockReduceSum(elements.z); 
1406:  
1407:                         if (threadIdx.x==0) { 
1408:                                 d_outArray[blockIdx.x].x = elements.x; 
1409:                                 d_outArray[blockIdx.x].y = elements.y; 
1410:                                 d_outArray[blockIdx.x].z = elements.z; 
1411:                         } 
1412:  
1413:                         tid += blockDim.x * gridDim.x; 
1414:                 } 
1415:         } 
1416:  
1417:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, double3 *d_outArray)  
1418:         { 
1419:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1420:  
1421:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) { 
1422:  
1423:                         double3 elements; 
1424:                         elements.x = 0.0; 
1425:                         elements.y = 0.0; 
1426:                         elements.z = 0.0; 
1427:  
1428:                         int index = ((tid + roundMaxSite) % roundMaxSite); 
1429:                         int thisBody = tid/roundMaxSite; 
1430:                         int currentBody = m_d_largeRigidIndices[thisBody]; 
1431:  
1432:                         int sectionSize = outSize/(*m_d_nLargeRigidBody); 
1433:                         if (index < sectionSize) { 
1434:                                 elements.x = d_outArray[index+thisBody*sectionSize].x; 
1435:                                 elements.y = d_outArray[index+thisBody*sectionSize].y; 
1436:                                 elements.z = d_outArray[index+thisBody*sectionSize].z; 
1437:                         } 
1438:  
1439:                         elements.x = blockReduceSum(elements.x); 
1440:                         elements.y = blockReduceSum(elements.y); 
1441:                         elements.z = blockReduceSum(elements.z); 
1442:  
1443:                         if (threadIdx.x==0) { 
1444:                                 if (nBlocks == (*m_d_nLargeRigidBody)) { 
1445:                                         m_d_gkRigid[3*currentBody] = elements.x; 
1446:                                         m_d_gkRigid[3*currentBody+1] = elements.y; 
1447:                                         m_d_gkRigid[3*currentBody+2] = elements.z; 
1448:                                 } 
1449:                                 else { 
1450:                                         d_outArray[blockIdx.x].x = elements.x; 
1451:                                         d_outArray[blockIdx.x].y = elements.y; 
1452:                                         d_outArray[blockIdx.x].z = elements.z; 
1453:                                 } 
1454:                         } 
1455:  
1456:                         tid += blockDim.x * gridDim.x; 
1457:                 } 
1458:         } 
1459:  
1460:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices,  
1461:                         double3 *d_outArray, const int *m_d_nRigidBody)  
1462:         { 
1463:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1464:  
1465:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) { 
1466:  
1467:                         double3 elements; 
1468:                         elements.x = 0.0; 
1469:                         elements.y = 0.0; 
1470:                         elements.z = 0.0; 
1471:  
1472:                         int index = ((tid + roundMaxSite) % roundMaxSite); 
1473:                         int thisBody = tid/roundMaxSite; 
1474:                         int currentBody = m_d_largeRigidIndices[thisBody]; 
1475:  
1476:                         int sectionSize = outSize/(*m_d_nLargeRigidBody); 
1477:                         if (index < sectionSize) { 
1478:                                 elements.x = d_outArray[index+thisBody*sectionSize].x; 
1479:                                 elements.y = d_outArray[index+thisBody*sectionSize].y; 
1480:                                 elements.z = d_outArray[index+thisBody*sectionSize].z; 
1481:                         } 
1482:  
1483:                         elements.x = blockReduceSum(elements.x); 
1484:                         elements.y = blockReduceSum(elements.y); 
1485:                         elements.z = blockReduceSum(elements.z); 
1486:  
1487:                         if (threadIdx.x==0) { 
1488:                                 if (nBlocks == (*m_d_nLargeRigidBody)) { 
1489:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody] = elements.x; 
1490:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+1] = elements.y; 
1491:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+2] = elements.z; 
1492:                                 } 
1493:                                 else { 
1494:                                         d_outArray[blockIdx.x].x = elements.x; 
1495:                                         d_outArray[blockIdx.x].y = elements.y; 
1496:                                         d_outArray[blockIdx.x].z = elements.z; 
1497:                                 } 
1498:                         } 
1499:  
1500:                         tid += blockDim.x * gridDim.x; 
1501:                 } 
1502:         } 
1503:  
1504:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *d_torques, const int *m_d_largeRigidIndices, double3 *d_outArray, const int *m_d_nRigidBody)  
1505:         { 
1506:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
1507:  
1508:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) { 
1509:  
1510:                         double3 elements; 
1511:                         elements.x = 0.0; 
1512:                         elements.y = 0.0; 
1513:                         elements.z = 0.0; 
1514:  
1515:                         int index = ((tid + roundMaxSite) % roundMaxSite); // Technically unnecessary with one body at a time - equivalent to tid 
1516:                         int thisBody = tid/roundMaxSite; // Also unnecessary with one body - could just pass as argument 
1517:                         int currentBody = m_d_largeRigidIndices[thisBody]; 
1518:  
1519:                         int sectionSize = outSize/(*m_d_nLargeRigidBody); 
1520:                         if (index < sectionSize) { 
1521:                                 elements.x = d_outArray[index+thisBody*sectionSize].x; 
1522:                                 elements.y = d_outArray[index+thisBody*sectionSize].y; 
1523:                                 elements.z = d_outArray[index+thisBody*sectionSize].z; 
1524:                         } 
1525:  
1526:                         elements.x = blockReduceSum(elements.x); 
1527:                         elements.y = blockReduceSum(elements.y); 
1528:                         elements.z = blockReduceSum(elements.z); 
1529:  
1530:                         if (threadIdx.x==0) { 
1531:                                 if (nBlocks == (*m_d_nLargeRigidBody)) { 
1532:                                         d_torques[3*currentBody] = elements.x; 
1533:                                         d_torques[3*currentBody+1] = elements.y; 
1534:                                         d_torques[3*currentBody+2] = elements.z; 
1535:                                 } 
1536:                                 else { 
1537:                                         d_outArray[blockIdx.x].x = elements.x; 
1538:                                         d_outArray[blockIdx.x].y = elements.y; 
1539:                                         d_outArray[blockIdx.x].z = elements.z; 
1540:                                 } 
1541:                         } 
1542:  
1543:                         tid += blockDim.x * gridDim.x; 
1544:                 } 
1545:  
1546:         } 
1547:  
1548:         __global__ void fullReduce(double *in, double* out, int N)  
1549:         { 
1550:                 double sum = 0; 
1551:                 //reduce multiple elements per thread 
1552:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) { 
1553:                         sum += in[i]; 
1554:                 } 
1555:                 sum = blockReduceSum(sum); 
1556:                 if (threadIdx.x==0) { 
1557:                         out[blockIdx.x]=sum; 
1558:                 } 
1559:         } 
1560:  
1561:         __global__ void fullDotReduce(double *in, double* out, int N) 
1562:         { 
1563:                 double sum = 0; 
1564:                 //reduce multiple elements per thread 
1565:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) { 
1566:                         sum += in[i]*in[i]; 
1567:                 } 
1568:                 sum = blockReduceSum(sum); 
1569:                 if (threadIdx.x==0) { 
1570:                         out[blockIdx.x]=sum; 
1571:                 } 
1572:         } 
1573:  
1574: }1198: }


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0