This project is read-only.

Floating Point Operations - A CUDAfy.NET Novice's Questions

Jul 28 at 10:08 PM
I have, for some time, been working with simulation of gravitational systems. I became curious about the possibility of using CUDAfy as a possible way to performing these computations quickly. The below is some code I wrote - using the CUDAfy example program right out of the box.

I know I need to read the book, but I was hoping that someone could point out what I am doing that is stupid in the below.

NOTE: The integration to find the next velocity and position of a body is grossly simplified. I figured that once I got it working, I could pull in the real algorithm.

Any help would be greatly appreciated. I am getting a run-time exception that a call to abs() - absolute value - is ambiguous. If I comment out this code, I get a run-time exception that "System.ArgumentException: Argument 1 of type System.Single[] is not on the GPU or not supported."

Thanks,

Doug
//////////////////////////////////////////

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Runtime.InteropServices;
using Cudafy;
using Cudafy.Host;
using Cudafy.Translator;

namespace CudafyByExample
{
public class Gravity
{
    private const int N = 2;
    private const int Sun = 0;

    public static int Execute()
    {
        CudafyModule km = CudafyTranslator.Cudafy();

        GPGPU gpu = CudafyHost.GetDevice(CudafyModes.Target, CudafyModes.DeviceId);

        gpu.LoadModule(km);

        float[] oldm = new float[N]; //mass in Kb
        float[] oldx = new float[N]; //distance in meters
        float[] oldy = new float[N]; //distance in meters
        float[] oldz = new float[N]; //distance in meters
        float[] oldvx = new float[N]; //veloctiy in meters/second
        float[] oldvy = new float[N]; //veloctiy in meters/second
        float[] oldvz = new float[N]; //veloctiy in meters/second

        float[] newm = new float[N]; //mass in Kb
        float[] newx = new float[N]; //distance in meters
        float[] newy = new float[N]; //distance in meters
        float[] newz = new float[N]; //distance in meters
        float[] newvx = new float[N]; //veloctiy in meters/second
        float[] newvy = new float[N]; //veloctiy in meters/second
        float[] newvz = new float[N]; //veloctiy in meters/second

        // allocate the memory on the GPU
        float[] dev_oldm = gpu.Allocate<float>(oldm);
        float[] dev_oldx = gpu.Allocate<float>(oldx);
        float[] dev_oldy = gpu.Allocate<float>(oldy);
        float[] dev_oldz = gpu.Allocate<float>(oldz);
        float[] dev_oldvx = gpu.Allocate<float>(oldvx);
        float[] dev_oldvy = gpu.Allocate<float>(oldvy);
        float[] dev_oldvz = gpu.Allocate<float>(oldvz);

        float[] dev_newm = gpu.Allocate<float>(newm);
        float[] dev_newx = gpu.Allocate<float>(newx);
        float[] dev_newy = gpu.Allocate<float>(newy);
        float[] dev_newz = gpu.Allocate<float>(newz);
        float[] dev_newvx = gpu.Allocate<float>(newvx);
        float[] dev_newvy = gpu.Allocate<float>(newvy);
        float[] dev_newvz = gpu.Allocate<float>(newvz);

        // fill the arrays 'a' and 'b' on the CPU
        for(int i = 0; i < N; i++)
        {
            oldm[i] = 0.0f; //mass in Kb
            oldx[i] = 0.0f; //distance in meters
            oldy[i] = 0.0f; //distance in meters
            oldz[i] = 0.0f; //distance in meters
            oldvx[i] = 0.0f; //veloctiy in meters/second
            oldvy[i] = 0.0f; //veloctiy in meters/second
            oldvz[i] = 0.0f; //veloctiy in meters/second

            newm[i] = 0.0f; //mass in Kb
            newx[i] = 0.0f; //distance in meters
            newy[i] = 0.0f; //distance in meters
            newz[i] = 0.0f; //distance in meters
            newvx[i] = 0.0f; //veloctiy in meters/second
            newvy[i] = 0.0f; //veloctiy in meters/second
            newvz[i] = 0.0f; //veloctiy in meters/second
        }

        // first index is the Sun
        oldm[0] = 1.988544E+30f;
        oldx[0] = 0.0f;
        oldy[0] = 0.0f;
        oldz[0] = 0.0f;
        oldvx[0] = 0.0f;
        oldvy[0] = 0.0f;
        oldvz[0] = 0.0f;

        // second index is the Earth
        oldm[1] = 5.97219E+24f;
        oldx[1] = 80571965975.3481f;
        oldy[1] = -128838109196.785f;
        oldz[1] = -36775239.7319856f;
        oldvx[1] = 24780.1934561512f;
        oldvy[1] = 15693.4332480694f;
        oldvz[1] = -1.94217087996037f;

        //print out calculations
        PrintNewState(newm, newx, newy, newz, newvx, newvy, newvz);

        // copy the arrays 'a' and 'b' to the GPU
        gpu.CopyToDevice(oldm, dev_oldm);
        gpu.CopyToDevice(oldx, dev_oldx);
        gpu.CopyToDevice(oldy, dev_oldy);
        gpu.CopyToDevice(oldz, dev_oldz);
        gpu.CopyToDevice(oldx, dev_oldvx);
        gpu.CopyToDevice(oldy, dev_oldvy);
        gpu.CopyToDevice(oldz, dev_oldvz);

        gpu.CopyToDevice(newm, dev_newm);
        gpu.CopyToDevice(newx, dev_newx);
        gpu.CopyToDevice(newy, dev_newy);
        gpu.CopyToDevice(newz, dev_newz);
        gpu.CopyToDevice(newx, dev_newvx);
        gpu.CopyToDevice(newy, dev_newvy);
        gpu.CopyToDevice(newz, dev_newvz);

        //turn the crank
        gpu.Launch(1, N).Newton(oldm, oldx, oldy, oldz, oldvx, oldvy, oldvz,
                                newm, newx, newy, newz, newvx, newvy, newvz);

        // copy the array 'c' back from the GPU to the CPU
        gpu.CopyFromDevice(dev_oldm, oldm);
        gpu.CopyFromDevice(dev_oldx, oldx);
        gpu.CopyFromDevice(dev_oldy, oldy);
        gpu.CopyFromDevice(dev_oldz, oldz);
        gpu.CopyFromDevice(dev_oldvx, oldvx);
        gpu.CopyFromDevice(dev_oldvy, oldvy);
        gpu.CopyFromDevice(dev_oldvz, oldvz);

        gpu.CopyFromDevice(dev_newm, newm);
        gpu.CopyFromDevice(dev_newx, newx);
        gpu.CopyFromDevice(dev_newy, newy);
        gpu.CopyFromDevice(dev_newz, newz);
        gpu.CopyFromDevice(dev_newvx, newvx);
        gpu.CopyFromDevice(dev_newvy, newvy);
        gpu.CopyFromDevice(dev_newvz, newvz);

        //print out calculations
        PrintNewState(newm, newx, newy, newz, newvx, newvy, newvz);

        // free the memory allocated on the GPU
        gpu.Free(dev_oldm);
        gpu.Free(dev_oldx);
        gpu.Free(dev_oldy);
        gpu.Free(dev_oldz);
        gpu.Free(dev_oldvx);
        gpu.Free(dev_oldvy);
        gpu.Free(dev_oldvz);

        gpu.Free(dev_newm);
        gpu.Free(dev_newx);
        gpu.Free(dev_newy);
        gpu.Free(dev_newz);
        gpu.Free(dev_newvx);
        gpu.Free(dev_newvy);
        gpu.Free(dev_newvz);

        return 0;
    }

    [Cudafy]
    public static void Newton(GThread thread, float[] oldm, float[] oldx, float[] oldy, float[] oldz, float[] oldvx, float[] oldvy, float[] oldvz,
                                            float[] newm, float[] newx, float[] newy, float[] newz, float[] newvx, float[] newvy, float[] newvz)
    {
        int tid = thread.blockIdx.x;
        if (tid < N)
        {
            float G = 6.674E-11f;
            float t = 1.0f;
            float fx = 0.0f;
            float fy = 0.0f;
            float fz = 0.0f;
            for (int i = 0; i < N; i++)
            {
                if (i != tid)
                {
                    float dx = GMath.Abs(oldx[tid] - oldx[i]);
                    fx += (G * oldm[i] * oldm[tid]) / (dx * dx);
                    float dy = GMath.Abs(oldy[tid] - oldy[i]);
                    fy += (G * oldm[i] * oldm[tid]) / (dy * dy);
                    float dz = GMath.Abs(oldz[tid] - oldz[i]);
                    fz += (G * oldm[i] * oldm[tid]) / (dz * dz);
                }
            }

            //f = ma or a = F/m
            float ax = fx / oldm[tid];
            float ay = fy / oldm[tid];
            float az = fz / oldm[tid];

            newm[tid] = oldm[tid];
            //v = at
            newvx[tid] = oldvx[tid] + (ax * t);
            newvy[tid] = oldvy[tid] + (ay * t);
            newvz[tid] = oldvz[tid] + (az * t);
            //d = 0.5at^2
            newx[tid] = oldx[tid] + (0.5f * ax * t * t);
            newy[tid] = oldy[tid] + (0.5f * ay * t * t);
            newz[tid] = oldz[tid] + (0.5f * az * t * t);
        }
    }

    public static void PrintNewState(float[] newm, float[] newx, float[] newy, float[] newz, float[] newvx, float[] newvy, float[] newvz)
    {
        for (int i = 0; i < N; i++)
        {
            Console.Write("Body " + i.ToString("##"));
            Console.Write("P(" + newx[i].ToString("E3") + ", " + newy[i].ToString("E3") + ", " + newz[i].ToString("E3") + ")");
            Console.Write("V(" + newvx[i].ToString("E3") + ", " + newvy[i].ToString("E3") + ", " + newvz[i].ToString("E3") + ")");
            Console.WriteLine("");
        }
    }
}
}