Performing operations on CUDA matrices while reading from a global Point

Question

Hey there, I have a mathematical function (multidimensional which means that there's an index which I pass to the C++-function on which single mathematical function I want to return. E.g. let's say I have a mathematical function like that:

f = Vector(x^2*y^2 / y^2 / x^2*z^2)

I would implement it like that:

double myFunc(int function_index)
{       
    switch(function_index)
    {
        case 1:
            return PNT[0]*PNT[0]*PNT[1]*PNT[1];
        case 2:
            return PNT[1]*PNT[1];
        case 3:
            return PNT[2]*PNT[2]*PNT[1]*PNT[1];
    }
}

whereas PNT is defined globally like that: double PNT[ NUM_COORDINATES ]. Now I want to implement the derivatives of each function for each coordinate thus generating the derivative matrix (columns = coordinates; rows = single functions). I wrote my kernel already which works so far and which call's myFunc().

The Problem is: For calculating the derivative of the mathematical sub-function i concerning coordinate j, I would use in sequential mode (on CPUs e.g.) the following code (whereas this is simplified because usually you would decrease h until you reach a certain precision of your derivative):

f0 = myFunc(i);
PNT[ j ] += h;
derivative = (myFunc(j)-f0)/h;
PNT[ j ] -= h;

now as I want to do this on the GPU in parallel, the problem is coming up: What to do with PNT? As I have to increase certain coordinates by h, calculate the value and than decrease it again, there's a problem coming up: How to do it without 'disturbing' the other threads? I can't modify PNT because other threads need the 'original' point to modify their own coordinate.

The second idea I had was to save one modified point for each thread but I discarded this idea quite fast because when using some thousand threads in parallel, this is a quite bad and probably slow (perhaps not realizable at all because of memory limits) idea.

'FINAL' SOLUTION So how I do it currently is the following, which adds the value 'add' on runtime (without storing it somewhere) via preprocessor macro to the coordinate identified by coord_index.

#define X(n) ((coordinate_index == n) ? (PNT[n]+add) : PNT[n])
__device__ double myFunc(int function_index, int coordinate_index, double add)
{           
    //*// Example: f[i] = x[i]^3
    return (X(function_index)*X(function_index)*X(function_index));
    // */
}

That works quite nicely and fast. When using a derivative matrix with 10000 functions and 10000 coordinates, it just takes like 0.5seks. PNT is defined either globally or as constant memory like __constant__ double PNT[ NUM_COORDINATES ];, depending on the preprocessor variable USE_CONST. The line return (X(function_index)*X(function_index)*X(function_index)); is just an example where every sub-function looks the same scheme, mathematically spoken:

f = Vector(x0^3 / x1^3 / ... / xN^3)

NOW THE BIG PROBLEM ARISES:

myFunc is a mathematical function which the user should be able to implement as he likes to. E.g. he could also implement the following mathematical function:

f = Vector(x0^2*x1^2*...*xN^2 / x0^2*x1^2*...*xN^2 / ... / x0^2*x1^2*...*xN^2)

thus every function looking the same. You as a programmer should only code once and not depending on the implemented mathematical function. So when the above function is being implemented in C++, it looks like the following:

__device__ double myFunc(int function_index, int coordinate_index, double add)
{           
    double ret = 1.0;
    for(int i = 0; i < NUM_COORDINATES; i++)
        ret *= X(i)*X(i);
    return ret; 
}

And now the memory accesses are very 'weird' and bad for performance issues because each thread needs access to each element of PNT twice. Surely, in such a case where each function looks the same, I could rewrite the complete algorithm which surrounds the calls to myFunc, but as I stated already: I don't want to code depending on the user-implemented function myFunc...

Could anybody come up with an idea how to solve this problem?? Thanks!

Answer 1

Rewinding back to the beginning and starting with a clean sheet, it seems you want to be able to do two things

• compute an arbitrary scalar valued function over an input array
• approximate the partial derivative of an arbitrary scalar valued function over the input array using first order accurate finite differencing

While the function is scalar valued and arbitrary, it seems that there are, in fact, two clear forms which this function can take:

1. A scalar valued function with scalar arguments
2. A scalar valued function with vector arguments

You appeared to have started with the first type of function and have put together code to deal with computing both the function and the approximate derivative, and are now wrestling with the problem of how to deal with the second case using the same code.

If this is a reasonable summary of the problem, then please indicate so in a comment and I will continue to expand it with some code samples and concepts. If it isn't, I will delete it in a few days.

In comments, I have been trying to suggest that conflating the first type of function with the second is not a good approach. The requirements for correctness in parallel execution, and the best way of extracting parallelism and performance on the GPU are very different. You would be better served by treating both types of functions separately in two different code frameworks with different usage models. When a given mathematical expression needs to be implemented, the "user" should make a basic classification as to whether that expression is like the model of the first type of function, or the second. The act of classification is what drives algorithmic selection in your code. This type of "classification by algorithm" is almost universal in well designed libraries - you can find it in C++ template libraries like Boost and the STL, and you can find it in legacy Fortran codes like the BLAS.