Quantum Mechanics on the GPU

Richard Gass

Department of Physics
University of Cincinnati
Cincinnati, OH 45211

Adam Simpson

Oak Ridge National Laboratory
Oak Ridge, TN 37831

Why Run on a GPU?

Speed

•	GPU’s have many more cores than CPU’s

On this machine I have two cores on my CPU but 16 cores on my GPU

• Only works for some problems

Problem needs to be highly parallel and each piece should have a low memory footprint

S[p_,f_]=1/(1-f+f/p);
Manipulate[Plot[S[n,f],{f,0,1},FrameLabel->{"fraction of code run in parallel","Speedup"},PlotLabel->"Amdahl's Law",RotateLabel->True,Frame->True,PlotRange->{1,16}],{{n,2,"Processors"},2,16,2,Appearance->"Labeled"},SaveDefinitions->True]

Why you should not run a GPU

Programming GPU’s is much harder

There maybe other ways to get the speed you need

• Compile into C
• Load a library function

Matrix Addition

In Mathematica

C=A+B

In OpenCL

/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */

// *********************************************************************
// oclVectorAdd Notes:  
//
// A simple OpenCL API demo application that implements 
// element by element vector addition between 2 float arrays. 
//
// Runs computations with OpenCL on the GPU device and then checks results 
// against basic host CPU/C++ computation.
//
// Uses some 'shr' and 'ocl' functions from oclUtils and shrUtils libraries for 
// compactness, but these are NOT required libs for OpenCL developement in general.
// *********************************************************************

// common SDK header for standard utilities and system libs 
#include <oclUtils.h>
#include <shrQATest.h>

// Name of the file with the source code for the computation kernel
// *********************************************************************
const char* cSourceFile = "VectorAdd.cl";

// Host buffers for demo
// *********************************************************************
void *srcA, *srcB, *dst;        // Host buffers for OpenCL test
void* Golden;                   // Host buffer for host golden processing cross check

// OpenCL Vars
cl_context cxGPUContext;        // OpenCL context
cl_command_queue cqCommandQueue;// OpenCL command que
cl_platform_id cpPlatform;      // OpenCL platform
cl_device_id cdDevice;          // OpenCL device
cl_program cpProgram;           // OpenCL program
cl_kernel ckKernel;             // OpenCL kernel
cl_mem cmDevSrcA;               // OpenCL device source buffer A
cl_mem cmDevSrcB;               // OpenCL device source buffer B 
cl_mem cmDevDst;                // OpenCL device destination buffer 
size_t szGlobalWorkSize;        // 1D var for Total # of work items
size_t szLocalWorkSize;		    // 1D var for # of work items in the work group	
size_t szParmDataBytes;			// Byte size of context information
size_t szKernelLength;			// Byte size of kernel code
cl_int ciErr1, ciErr2;			// Error code var
char* cPathAndName = NULL;      // var for full paths to data, src, etc.
char* cSourceCL = NULL;         // Buffer to hold source for compilation
const char* cExecutableName = NULL;

// demo config vars
int iNumElements = 11444777;	// Length of float arrays to process (odd # for illustration)
shrBOOL bNoPrompt = shrFALSE;  

// Forward Declarations
// *********************************************************************
void VectorAddHost(const float* pfData1, const float* pfData2, float* pfResult, int iNumElements);
void Cleanup (int argc, char **argv, int iExitCode);

// Main function 
// *********************************************************************
int main(int argc, char **argv)
{
    shrQAStart(argc, argv);

    // get command line arg for quick test, if provided
    bNoPrompt = shrCheckCmdLineFlag(argc, (const char**)argv, "noprompt");
    
    // start logs 
	cExecutableName = argv[0];
    shrSetLogFileName ("oclVectorAdd.txt");
    shrLog("%s Starting...\n\n# of float elements per Array \t= %i\n", argv[0], iNumElements); 

    // set and log Global and Local work size dimensions
    szLocalWorkSize = 256;
    szGlobalWorkSize = shrRoundUp((int)szLocalWorkSize, iNumElements);  // rounded up to the nearest multiple of the LocalWorkSize
    shrLog("Global Work Size \t\t= %u\nLocal Work Size \t\t= %u\n# of Work Groups \t\t= %u\n\n", 
           szGlobalWorkSize, szLocalWorkSize, (szGlobalWorkSize % szLocalWorkSize + szGlobalWorkSize/szLocalWorkSize)); 

    // Allocate and initialize host arrays 
    shrLog( "Allocate and Init Host Mem...\n"); 
    srcA = (void *)malloc(sizeof(cl_float) * szGlobalWorkSize);
    srcB = (void *)malloc(sizeof(cl_float) * szGlobalWorkSize);
    dst = (void *)malloc(sizeof(cl_float) * szGlobalWorkSize);
    Golden = (void *)malloc(sizeof(cl_float) * iNumElements);
    shrFillArray((float*)srcA, iNumElements);
    shrFillArray((float*)srcB, iNumElements);

    //Get an OpenCL platform
    ciErr1 = clGetPlatformIDs(1, &cpPlatform, NULL);

    shrLog("clGetPlatformID...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clGetPlatformID, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    //Get the devices
    ciErr1 = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &cdDevice, NULL);
    shrLog("clGetDeviceIDs...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clGetDeviceIDs, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    //Create the context
    cxGPUContext = clCreateContext(0, 1, &cdDevice, NULL, NULL, &ciErr1);
    shrLog("clCreateContext...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clCreateContext, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Create a command-queue
    cqCommandQueue = clCreateCommandQueue(cxGPUContext, cdDevice, 0, &ciErr1);
    shrLog("clCreateCommandQueue...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clCreateCommandQueue, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Allocate the OpenCL buffer memory objects for source and result on the device GMEM
    cmDevSrcA = clCreateBuffer(cxGPUContext, CL_MEM_READ_ONLY, sizeof(cl_float) * szGlobalWorkSize, NULL, &ciErr1);
    cmDevSrcB = clCreateBuffer(cxGPUContext, CL_MEM_READ_ONLY, sizeof(cl_float) * szGlobalWorkSize, NULL, &ciErr2);
    ciErr1 |= ciErr2;
    cmDevDst = clCreateBuffer(cxGPUContext, CL_MEM_WRITE_ONLY, sizeof(cl_float) * szGlobalWorkSize, NULL, &ciErr2);
    ciErr1 |= ciErr2;
    shrLog("clCreateBuffer...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clCreateBuffer, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }
    
    // Read the OpenCL kernel in from source file
    shrLog("oclLoadProgSource (%s)...\n", cSourceFile); 
    cPathAndName = shrFindFilePath(cSourceFile, argv[0]);
    cSourceCL = oclLoadProgSource(cPathAndName, "", &szKernelLength);

    // Create the program
    cpProgram = clCreateProgramWithSource(cxGPUContext, 1, (const char **)&cSourceCL, &szKernelLength, &ciErr1);
    shrLog("clCreateProgramWithSource...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clCreateProgramWithSource, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Build the program with 'mad' Optimization option
    #ifdef MAC
        char* flags = "-cl-fast-relaxed-math -DMAC";
    #else
        char* flags = "-cl-fast-relaxed-math";
    #endif
    ciErr1 = clBuildProgram(cpProgram, 0, NULL, NULL, NULL, NULL);
    shrLog("clBuildProgram...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clBuildProgram, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Create the kernel
    ckKernel = clCreateKernel(cpProgram, "VectorAdd", &ciErr1);
    shrLog("clCreateKernel (VectorAdd)...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clCreateKernel, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Set the Argument values
    ciErr1 = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), (void*)&cmDevSrcA);
    ciErr1 |= clSetKernelArg(ckKernel, 1, sizeof(cl_mem), (void*)&cmDevSrcB);
    ciErr1 |= clSetKernelArg(ckKernel, 2, sizeof(cl_mem), (void*)&cmDevDst);
    ciErr1 |= clSetKernelArg(ckKernel, 3, sizeof(cl_int), (void*)&iNumElements);
    shrLog("clSetKernelArg 0 - 3...\n\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clSetKernelArg, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // --------------------------------------------------------
    // Start Core sequence... copy input data to GPU, compute, copy results back

    // Asynchronous write of data to GPU device
    ciErr1 = clEnqueueWriteBuffer(cqCommandQueue, cmDevSrcA, CL_FALSE, 0, sizeof(cl_float) * szGlobalWorkSize, srcA, 0, NULL, NULL);
    ciErr1 |= clEnqueueWriteBuffer(cqCommandQueue, cmDevSrcB, CL_FALSE, 0, sizeof(cl_float) * szGlobalWorkSize, srcB, 0, NULL, NULL);
    shrLog("clEnqueueWriteBuffer (SrcA and SrcB)...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clEnqueueWriteBuffer, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Launch kernel
    ciErr1 = clEnqueueNDRangeKernel(cqCommandQueue, ckKernel, 1, NULL, &szGlobalWorkSize, &szLocalWorkSize, 0, NULL, NULL);
    shrLog("clEnqueueNDRangeKernel (VectorAdd)...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clEnqueueNDRangeKernel, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }

    // Synchronous/blocking read of results, and check accumulated errors
    ciErr1 = clEnqueueReadBuffer(cqCommandQueue, cmDevDst, CL_TRUE, 0, sizeof(cl_float) * szGlobalWorkSize, dst, 0, NULL, NULL);
    shrLog("clEnqueueReadBuffer (Dst)...\n\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        shrLog("Error in clEnqueueReadBuffer, Line %u in file %s !!!\n\n", __LINE__, __FILE__);
        Cleanup(argc, argv, EXIT_FAILURE);
    }
    //--------------------------------------------------------

    // Compute and compare results for golden-host and report errors and pass/fail
    shrLog("Comparing against Host/C++ computation...\n\n"); 
    VectorAddHost ((const float*)srcA, (const float*)srcB, (float*)Golden, iNumElements);
    shrBOOL bMatch = shrComparefet((const float*)Golden, (const float*)dst, (unsigned int)iNumElements, 0.0f, 0);

    // Cleanup and leave
    Cleanup (argc, argv, (bMatch == shrTRUE) ? EXIT_SUCCESS : EXIT_FAILURE);
}

void Cleanup (int argc, char **argv, int iExitCode)
{
    // Cleanup allocated objects
    shrLog("Starting Cleanup...\n\n");
    if(cPathAndName)free(cPathAndName);
    if(cSourceCL)free(cSourceCL);
	if(ckKernel)clReleaseKernel(ckKernel);  
    if(cpProgram)clReleaseProgram(cpProgram);
    if(cqCommandQueue)clReleaseCommandQueue(cqCommandQueue);
    if(cxGPUContext)clReleaseContext(cxGPUContext);
    if(cmDevSrcA)clReleaseMemObject(cmDevSrcA);
    if(cmDevSrcB)clReleaseMemObject(cmDevSrcB);
    if(cmDevDst)clReleaseMemObject(cmDevDst);

    // Free host memory
    free(srcA); 
    free(srcB);
    free (dst);
    free(Golden);

    // finalize logs and leave
    shrQAFinishExit(argc, (const char **)argv, (iExitCode == EXIT_SUCCESS) ? QA_PASSED : QA_FAILED);
}

// "Golden" Host processing vector addition function for comparison purposes
// *********************************************************************
void VectorAddHost(const float* pfData1, const float* pfData2, float* pfResult, int iNumElements)
{
    int i;
    for (i = 0; i < iNumElements; i++) 
    {
        pfResult[i] = pfData1[i] + pfData2[i]; 
    }
}

But it’s mostly “boilerplate” the kernel is just

/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */
 
 // OpenCL Kernel Function for element by element vector addition
__kernel void VectorAdd(__global const float* a, __global const float* b, __global float* c, int iNumElements)
{
    // get index into global data array
    int iGID = get_global_id(0);

    // bound check (equivalent to the limit on a 'for' loop for standard/serial C code
    if (iGID >= iNumElements)
    {   
        return; 
    }
    
    // add the vector elements
    c[iGID] = a[iGID] + b[iGID];
}

Mathematica good at writing boilerplate so just write the kernel and let CUDALink or OpenCL link do the rest.

Eigenfunction expansions in QM

We want the eigenvalues and eigenvectors A ψ=λ ψ.

Expand ψ as ψ=(∑_n)(Subscript[a, n])(Subscript[ϕ, n]) .

Then eigenvalues (Subscript[λ, n]) are given by the eigenvalues of the matrix M defined by (Subscript[M, m,n])= <(Subscript[ϕ, m])|A (Subscript[ϕ, n]) > and the eigenfunctions ofA are given by the inner product of the eigenfunctions of M with the basis functions ϕ.

If we have n basis functions then we need to do (n^2) integrals to compute M. Computing the integrals is were the code spends it time. Finding the eigenvalues is fast.

Timing[Eigenvalues[Table[RandomReal[],{i,1,100},{j,1,100}]];]

All the integrals are independent so compute them on the GPU.

Development Cycle

We want to compute (on a GPU)

(∫_(-L/2)^(L/2))(Subscript[ϕ, i])V(x)(Subscript[ϕ, j])dx

where the (Subscript[ϕ, i]) are particle in a box basis states: (Subscript[ϕ, i])=(Sqrt[2/L])Sin(i π (x+L/2)/L) and V(x) is the potential.

Prototype in Mathematica

Our basis functions are

Clear[L,ϕ,i]
ϕ[i_,x_]=(Sqrt[2/L]) Sin[i π (x+L/2)/L]

and the integrand we want to evaluate is:

Clear[f,n,n,𝒱,L]
f[x_,n_,m_]=(2/L) Sin[n π (x-L/2)/L] Sin[m π (x-L/2)/ L] 𝒱[x]

As a first example we will work with a harmonic potential since we know the exact result for this case.

𝒱[x_]=(x^2)/2

Since our C code will use the Composite Simpson’s rule

 ((∫_a)^b)f(x) dx=(h/6)(f(a)+f(b)+2(∑_(i=1)^(N-1))f((Subscript[x, i]))+4(∑_(i=1)^N)f((Subscript[x, i-1/2]))) 
 
 we implement this in Mathematica to compute s[1,1] which is the integral of f(x,1,1).

Clear[s,i,j]
s[i_,j_]=f[-L/2,i,j]+f[L/2,i,j];
h=0.1;
L=1.0;
For[i=1,i≤9,i++,s[1,1]=s[1,1]+4 f[-L/2+(i-0.5) h,1,1]+2.0 f[-L/2+i h,1,1]];
result=h(s[1,1]+4.0 f[-L/2+9.5 h,1,1])/6.0
Clear[i]

We can check this using NIntegrate which shows that the results are very close.

NIntegrate[f[x,1,1],{x,-L/2,L/2}]

Write the code in C

To compute the general case we define f(x,n,m) and use two For loops to loop through the allowed values of n and m. In the example below we have restricted m and n to be less than two.

/* Simpson.c */
# include <stdio.h>
# include <math.h>
# define V(x) ((x)*(x)/2)
# define PI 3.14159
# define f(x,a,b) (2.0*sin(a*PI*(x+0.5))*V(x)*sin(b*PI*(x+0.5)))
main()
{int i,m,n;
float h,s;
for(m=1;m<=2;m++){for(n=1;n<=2;n++){s=f(-0.5,n,m)+f(0.5,n,m);
h=1.0/10.0;
for (i=1;i<=9;i++){
s=s+4.0*f(-0.5+(i-0.5)*h,n,m)+2.0*f(-0.5+i*h,n,m);
}
s=s+4.0*f(-0.5+9.5*h,n,m);
s=s*h/6.0;
printf("integral=%f\n",s);}}
}

Running the code produces the following output.

integral=0.016341
integral=-0.000000
integral=-0.000000
integral=0.035352

We get good agreement with NIntegrate so we are ready to port the code to the GPU.

Table[NIntegrate[2Sin[n π (x+L/2)]((x^2)/2) Sin[m π (x+L/2)],{x,0-L/2,L/2},AccuracyGoal->6],{n,1,2},{m,1,2}]

Port to OpenCL or CUDA

Needs["OpenCLLink`"]

The script that we want is

Clear[L]
srcintegrate="

#define V(x) (x)*(x)/2
#define f(x,i,j) (2.0/L) * sin((float)(i * p *( x+L/2.0 )/L )) * sin((float)(j * p * (x+L/2.0)/L))*V(x)
constant float p=3.14159;

__kernel void integrals(__global float * C , float L, int mesh, int n, __local float * s_vals) 
{
	// i and j identify our work groups(thread blocks) and range from 1 to n
	// these are the basis indices
	// each thread block calculates one basis integral
	// total number of grid points is limited by the largest block size supported by the GPU
	int i = get_group_id(0) + 1;
	int j = get_group_id(1) + 1;

	// l runs from 1 to mesh and represents a 1D mesh
	int l = get_local_id(0) + 1;

	float s = 0.0;
	// add 1 so that a thread in each block isn't sitting around idle
	mesh = mesh+1;
	float h = L/mesh;

    // Each thread calculates a piece of series
    s = 4.0*f(-L/2.0+(l-0.5)*h,i,j)+2.0*f(-L/2.0+l*h,i,j);	

    // Each thread writes the value to shared memory
    s_vals[l-1] = s;

    // Sync up all threads in a given block to ensure all have finished writting their values to shared memory
    barrier(CLK_LOCAL_MEM_FENCE);
	
	// Thread 1 in each block sums up values in shared memory and writes to global memory
	if(l == 1){
		float sum = f(-L/2.0,i,j)+f(L/2.0,i,j) + 4.0*f(-L/2.0+(mesh-0.5)*h,i,j);
	// get_local_size(dimensionx) gives the number of local work-items in //dimension dimensionx
		for(uint x = 0; x<get_local_size(0); x++) {
			sum += s_vals[x];
		}
		C[get_group_id(1)*get_num_groups(0) + get_group_id(0)] = sum*h/6.0;
	} 
}
";

gridPoints = 16;
numberOfBasisFunctions= 16;
blockSize=gridPoints;
threadsInX = gridPoints*numberOfBasisFunctions;
threadsInY = numberOfBasisFunctions;
SharedMemoryBytes = blockSize*4;

mem=OpenCLMemoryAllocate["Float",{numberOfBasisFunctions,numberOfBasisFunctions}];
integrals=OpenCLFunctionLoad[srcintegrate,"integrals",{{"Float"},"Float",_Integer,_Integer,"Shared"},{blockSize,1}];
t1=AbsoluteTime[];
ints = integrals[mem,10.0,gridPoints,numberOfBasisFunctions,SharedMemoryBytes,{threadsInX,threadsInY}];
V=OpenCLMemoryGet[mem];
AbsoluteTime[]-t1

Check the results

check=With[{L=10},Table[NIntegrate[(2/L)Sin[n π (x+L/2)/L]((x^2)/2) Sin[m π (x+L/2)/L],{x,0-L/2,L/2},AccuracyGoal->6],{n,1,numberOfBasisFunctions},{m,1,numberOfBasisFunctions}]];

Chop[V-check,(10^(-4))]

MatrixPlot[Chop[V-check,(10^(-4))]]

• The last term looks like it could be a problem.

Highly oscillatory function gives Simpson’s rule but not NIntegrate trouble

With[{L=10},Plot[(x^2)((Sin[16 π x/L])^2),{x,0,L}]]

Integrate the function into Mathematica

• Compute the kinetic energy terms

Clear[L,i,j]
a[i_,i_,L_]=(((π^2) (i^2) (ℏ^2))/(2 (L^2) m));
a[i_,j_,L_]=0;

• Compute M and find it’s eigenvalues

Clear[M]
ℏ=1;
m=1;
M[L_]:=Module[{},integrals[mem,L,gridPoints,numberOfBasisFunctions,SharedMemoryBytes,{threadsInX,threadsInY}];
V=OpenCLMemoryGet[mem];
Table[a[i,j,L]+V[[i,j]],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions}]]

Eigenvalues[M[10]]

exact=Table[(n+1/2),{n,15,0,-1}];

percentError=100(Eigenvalues[M[10]]-exact)/exact

One Dimensional Quantum Wells

A more interesting  well

v[x_]=(x^4)+.7 (x^3)-2.0(x^2)-1.13x-2.0

Plot[v[x],{x,-2,2},AxesLabel->{"x","V(x)"}]

But it would be better still if we could change the well

Manipulate[Plot[Tooltip[(x^4)+.7 (x^3)-c (x^2)-1.13x-2.0],{x,-2,2},AxesLabel->{"x","V(x)"}],{{c,2.0},3.0,1.0},SaveDefinitions->True]

Just change the definition of V in our script.

srcintegrate2="
#define V(x,c) ((x)*(x)*(x)*(x)+0.7*(x)*(x)*(x)-c*(x)*(x)-1.13*(x)-2.0)
#define f(x,i,j) (2.0/L) * sin((float)(i * p *( x+L/2.0 )/L )) * sin((float)(j * p * (x+L/2.0)/L))*V(x,c)
constant float p=3.14159;

__kernel void integrals(__global float * C , float L,float c, int mesh, int n, __local float * s_vals) 
{
	// i and j identify our work groups(thread blocks) and range from 1 to n
	// these are the basis indices
	// each thread block calculates one basis integral
	// total number of grid points is limited by the largest block size supported by the GPU
	int i = get_group_id(0) + 1;
	int j = get_group_id(1) + 1;

	// l runs from 1 to mesh and represents a 1D mesh
	int l = get_local_id(0) + 1;

	float s = 0.0;
	// add 1 so that a thread in each block isn't sitting around idle
	mesh = mesh+1;
	float h = L/mesh;

    // Each thread calculates a piece of series
    s = 4.0*f(-L/2.0+(l-0.5)*h,i,j)+2.0*f(-L/2.0+l*h,i,j);	

    // Each thread writes the value to shared memory
    s_vals[l-1] = s;

    // Sync up all threads in a given block to ensure all have finished writting their values to shared memory
    barrier(CLK_LOCAL_MEM_FENCE);
	
	// Thread 1 in each block sums up values in shared memory and writes to global memory
	if(l == 1){
		float sum = f(-L/2.0,i,j)+f(L/2.0,i,j) + 4.0*f(-L/2.0+(mesh-0.5)*h,i,j);
	// get_local_size(dimensionx) gives the number of local work-items in
    //dimension dimensionx
		for(uint x = 0; x<get_local_size(0); x++) {
			sum += s_vals[x];
		}
		C[get_group_id(1)*get_num_groups(0) + get_group_id(0)] = sum*h/6.0;
	} 
}
";

Load the function and check the result

gridPoints = 16;
numberOfBasisFunctions= 16;
blockSize=gridPoints;
threadsInX = gridPoints*numberOfBasisFunctions;
threadsInY =numberOfBasisFunctions;
SharedMemoryBytes = blockSize*4;

mem=OpenCLMemoryAllocate["Float",{numberOfBasisFunctions,numberOfBasisFunctions}];
integrals=OpenCLFunctionLoad[srcintegrate2,"integrals",{{"Float"},"Float",_Real,_Integer,_Integer,"Shared"},{blockSize,1},ShellOutputFunction->Print,ShellCommandFunction->Print];

Clear[M]
ℏ=1;
m=1;
M[L_,c_]:=Module[{},integrals[mem,L,c,gridPoints,numberOfBasisFunctions,SharedMemoryBytes,{threadsInX,threadsInY}];
V=OpenCLMemoryGet[mem];
Table[a[i,j,L]+V[[i,j]],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions}]]

Sort[Eigenvalues[M[6.0,2.0]]]

Make it interactive

Manipulate[Show[Plot[Tooltip[(x^4)+.7 (x^3)-c (x^2)-1.13x-2.0],{x,-3,3},PlotRange->{-5,3}],Graphics[Transpose[{{Blue,Green,Red},Map[Line[{{-2.1,#},{2.1,#}}]&,Take[Sort[Eigenvalues[M[6.0,c]]],3]]}]]],{{c,2.0},3.0,1.0}]

Clear[L,ϕ]
ϕ[i_,x_]=(Sqrt[2/L]) Sin[i π (x+L/2)/L];basisFunctions=Table[ϕ[i,x]/.L->6.0,{i,1,numberOfBasisFunctions}];
Manipulate[Show[Plot[Tooltip[(x^4)+.7 (x^3)-c (x^2)-1.13x-2.0],{x,-3,3},PlotRange->{-5,3},PlotStyle->Black],Plot[Evaluate[2p Power[Take[Map[Part[#,2]&,Sort[Transpose[Eigensystem[M[6,c]]],#1[[1]]<#2[[1]]&]].basisFunctions,{1,3}],p]],{x,-3,3},PlotStyle->{Blue,Green,Red}],Graphics[Transpose[{{Blue,Green,Red},Map[Line[{{-2.1,#},{2.1,#}}]&,Take[Sort[Eigenvalues[M[6,c]]],3]]}]],ImageSize->7 72],{{c,2},3,1},{p,1,2,1,ControlType->SetterBar}]

Summary

• Even with OpenCLLink or CUDALink GPU programming is harder than programming in C
• Much harder than in Mathematica and you give up the sophistication of functions such as NIntegrate
• But you can get significant gains in speed