hdiff output

r29336/gpu.cpp 2015-11-17 23:34:18.477557197 +0000 r29335/gpu.cpp 2015-11-17 23:34:18.681559934 +0000
9088:                                 ((alphaD + (EthreshD - (totdih))) * temp0 * KB);9088:                                 ((alphaD + (EthreshD - (totdih))) * temp0 * KB);
9089:                         fwgtd                                   = (alphaD * alphaD)/((alphaD + EthreshD - (totdih)) * (alphaD + EthreshD - (totdih)));9089:                         fwgtd                                   = (alphaD * alphaD)/((alphaD + EthreshD - (totdih)) * (alphaD + EthreshD - (totdih)));
9090:                 }9090:                 }
9091:         }9091:         }
9092: 9092: 
9093:         gpu->sim.AMDtboost                              = tboost;9093:         gpu->sim.AMDtboost                              = tboost;
9094:         gpu->pbAMDfwgtd->_pSysData[0]                   = fwgtd;9094:         gpu->pbAMDfwgtd->_pSysData[0]                   = fwgtd;
9095:         gpu->pbAMDfwgtd->Upload();9095:         gpu->pbAMDfwgtd->Upload();
9096: }9096: }
9097: #endif9097: #endif
 9098: 
 9099: // Numerical second derivative calculation for GMIN/OPTIM
 9100: extern "C" void gminoptim_hessian(int *natoms, double *coords, double *hessian, double *delta)
 9101: {
 9102:    // Allocate necessary handle, error and status variables
 9103:    cudaError_t cudaStat;
 9104:    cublasStatus_t status;
 9105:    cublasHandle_t handle;
 9106: 
 9107:    // Also allocate 1.0 and 0.5d, as we need to pass these in
 9108:    const double minusOne = -1.0;
 9109:    const double * pToMinusOne = &minusOne;
 9110: 
 9111:    // Copy coordinates from GMIN/OPTIM to the GPU
 9112:         gminoptim_coords_cputogpu(gpu, natoms, coords);
 9113: 
 9114:    // Store location of coordinate array pointers
 9115:    PMEFloat2 * coordsxy_single = gpu->pbAtomXYSP->_pDevData;
 9116:    PMEFloat * coordsz_single = gpu->pbAtomZSP->_pDevData;
 9117:    PMEDouble * coords_double = gpu->pbAtom->_pDevData;
 9118: 
 9119:    // Allocate some storage for the minus and plus components of the gradient
 9120:    double *d_forceMinus, *d_forcePlus;
 9121:    int forceSize;
 9122:    forceSize = sizeof(gpu->pbForce->_pDevData[0]) * gpu->pbForce->_length;
 9123:    cudaMalloc(&d_forceMinus, forceSize); 
 9124:    cudaMalloc(&d_forcePlus, forceSize); 
 9125: 
 9126:    // Allocate storage for the hessian itself
 9127:    double *d_hessian;
 9128:    int hessSize = sizeof(double) * gpu->pbForce->_length * gpu->pbForce->_length;
 9129:    cudaMalloc(&d_hessian, hessSize);
 9130: 
 9131:    // Copy delta to the GPU
 9132:    const double *d_delta;
 9133:    cudaStat = cudaMalloc(&d_delta, sizeof(double));
 9134:         cudaStat = cudaMemcpy(d_delta, delta, sizeof(double), cudaMemcpyHostToDevice); 
 9135: 
 9136:    // Copy half delta to the GPU
 9137:    const double *d_halfRecipDelta;
 9138:    const double halfRecipDelta = 0.5 / *delta;
 9139:    cudaStat = cudaMalloc(&d_halfRecipDelta, sizeof(double));
 9140:         cudaStat = cudaMemcpy(d_halfRecipDelta, &halfRecipDelta, sizeof(double), cudaMemcpyHostToDevice); 
 9141:    
 9142:    // Initialise cuBLAS
 9143:    status = cublasCreate(&handle);
 9144: 
 9145:    // Should loop through all x, then all y, then all z, because of the stride.
 9146:    // x coordinates
 9147:    for (int i = 0; i < *natoms; i++) {
 9148:       // Change coordinates to plus
 9149:       coordsxy_single[i].x += *d_delta;
 9150:       coords_double[i] += *d_delta;
 9151: 
 9152:       // Calculate forcePlus
 9153:       gpu_gb_forces_();
 9154:       // At this point, the forces should be in gpu->pbForce->_pDevData
 9155:       // There are pointers from gpu->sim->pForce{X,Y,Z} that provide access to the x, y and z components.
 9156:       cudaMemcpy(d_forcePlus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9157: 
 9158:       // Change coordinates to minus
 9159:       coordsxy_single[i].x -= 2.0 * *d_delta;
 9160:       coords_double[i] -= 2.0 * *d_delta;
 9161: 
 9162:       // Calculate forceMinus
 9163:       gpu_gb_forces_();
 9164:       cudaMemcpy(d_forceMinus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9165: 
 9166:       // force = -gradient, so dg/dx = (g_+ - g_-) / 2d = (f_- - f_+) / 2d
 9167:       // Y = aX + Y
 9168:       cublasDaxpy(handle, gpu->pbForce->_length, pToMinusOne, d_forcePlus, 1, d_forceMinus, 1);
 9169:       // X = aX
 9170:       cublasDscal(handle, gpu->pbForce->_length, d_halfRecipDelta, d_forceMinus, 1)
 9171: 
 9172:       // Copy the result (at d_forceMinus) into the hessian
 9173:       cudaMemcpy((d_hessian + i * natoms), d_forceMinus, forceSize, cudaMemcpyDeviceToDevice);
 9174:    }
 9175: 
 9176:    // y coordinates
 9177:    for (int j = *natoms; j < 2 * *natoms; j++) {
 9178:       // Change coordinates to plus
 9179:       coordsxy_single[j - *natoms].y += *d_delta;
 9180:       coords_double[j] += *d_delta;
 9181: 
 9182:       // Calculate forcePlus
 9183:       gpu_gb_forces_();
 9184:       // At this point, the forces should be in gpu->pbForce->_pDevData
 9185:       // There are pointers from gpu->sim->pForce{X,Y,Z} that provide access to the x, y and z components.
 9186:       cudaMemcpy(d_forcePlus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9187: 
 9188:       // Change coordinates to minus
 9189:       coordsxy_single[j - *natoms].y -= 2.0 * *d_delta;
 9190:       coords_double[j] -= 2.0 * *d_delta;
 9191: 
 9192:       // Calculate forceMinus
 9193:       gpu_gb_forces_();
 9194:       cudaMemcpy(d_forceMinus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9195: 
 9196:       // force = -gradient, so dg/dx = (g_+ - g_-) / 2d = (f_- - f_+) / 2d
 9197:       // Y = aX + Y
 9198:       cublasDaxpy(handle, gpu->pbForce->_length, pToMinusOne, d_forcePlus, 1, d_forceMinus, 1);
 9199:       // X = aX
 9200:       cublasDscal(handle, gpu->pbForce->_length, d_halfRecipDelta, d_forceMinus, 1)
 9201: 
 9202:       // Copy the result (at d_forceMinus) into the hessian
 9203:       cudaMemcpy((d_hessian + j * natoms), d_forceMinus, forceSize, cudaMemcpyDeviceToDevice);
 9204:    }
 9205: 
 9206:    // z coordinates
 9207:    for (int k = 2 * *natoms; k < 3 * *natoms; k++) {
 9208:       // Change coordinates to plus
 9209:       coordsz_single[k - 2 * *natoms] += *d_delta;
 9210:       coords_double[k] += *d_delta;
 9211: 
 9212:       // Calculate forcePlus
 9213:       gpu_gb_forces_();
 9214:       // At this point, the forces should be in gpu->pbForce->_pDevData
 9215:       // There are pointers from gpu->sim->pForce{X,Y,Z} that provide access to the x, y and z components.
 9216:       cudaMemcpy(d_forcePlus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9217: 
 9218:       // Change coordinates to minus
 9219:       coordsz_single[j - *natoms] -= 2.0 * *d_delta;
 9220:       coords_double[k] -= 2.0 * *d_delta;
 9221: 
 9222:       // Calculate forceMinus
 9223:       gpu_gb_forces_();
 9224:       cudaMemcpy(d_forceMinus, gpu->pbForce->_pDevData, forceSize, cudaMemcpyDeviceToDevice);
 9225: 
 9226:       // force = -gradient, so dg/dx = (g_+ - g_-) / 2d = (f_- - f_+) / 2d
 9227:       // Y = aX + Y
 9228:       cublasDaxpy(handle, gpu->pbForce->_length, pToMinusOne, d_forcePlus, 1, d_forceMinus, 1);
 9229:       // X = aX
 9230:       cublasDscal(handle, gpu->pbForce->_length, d_halfRecipDelta, d_forceMinus, 1)
 9231: 
 9232:       // Copy the result (at d_forceMinus) into the hessian
 9233:       cudaMemcpy((d_hessian + k * natoms), d_forceMinus, forceSize, cudaMemcpyDeviceToDevice);
 9234:    }
 9235: 
 9236:    // Copy the hessian from the GPU to CPU
 9237:    cudaStat = cudaMemcpy(hessian, d_hessian, hessSize, cudaMemcpyDeviceToHost);
 9238: 
 9239:    // Free allocated memory
 9240:    cudaFree(d_hessian);
 9241:    cudaFree(d_forceMinus);
 9242:    cudaFree(d_forcePlus);
 9243:    cudaFree(d_delta);
 9244:    cudaFree(d_halfRecipDelta);
 9245: }


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0