ls,
here the source of Julia in CUDA C, a slightly modified version of the
one from the book "CUDA by Example", by Sanders and Kandrot. See for
the source
http://developer.nvidia.com/cuda-example-introduction-general-purpose-gpu-programming
(chapter 4).
The key idea sugested by R.E. Boss in function juliaB (see programming
forum) is the ommision of the conditional expression, which can be
harmful for the speedup.
A first naive transcription into CUDA C gives some errors. As soon as
I have something working I will inform this forum.
Jan.
/*
* Copyright 1993-2010 NVIDIA Corporation. All rights reserved.
*
* NVIDIA Corporation and its licensors retain all intellectual property and
* proprietary rights in and to this software and related documentation.
* Any use, reproduction, disclosure, or distribution of this software
* and related documentation without an express license agreement from
* NVIDIA Corporation is strictly prohibited.
*
* Please refer to the applicable NVIDIA end user license agreement (EULA)
* associated with this source code for terms and conditions that govern
* your use of this NVIDIA software.
*
*/
#include "../common/book.h"
#include "../common/cpu_bitmap.h"
#include <math.h>
#define DIM 1000
struct cuComplex {
float r;
float i;
__device__ cuComplex( float a, float b ) : r(a), i(b) {}
__device__ float magnitude2( void ) {
return r * r + i * i;
}
__device__ cuComplex operator*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
__device__ cuComplex operator+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
};
__device__ int julia( int x, int y ) {
const float scale = 1.5;
float jx = scale * (float)(DIM/2 - x)/(DIM/2);
float jy = scale * (float)(DIM/2 - y)/(DIM/2);
cuComplex c(-0.8, 0.156);
cuComplex a(jx, jy);
int i = 0;
for (i=0; i<200; i++) {
a = a * a + c;
if (a.magnitude2() > 1000)
return i;
}
return 255;
}
// declare device arrays
__device__ unsigned char dev_r[256];
__device__ unsigned char dev_g[256];
__device__ unsigned char dev_b[256];
__global__ void kernel( unsigned char *ptr ) {
// map from blockIdx to pixel position
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y * gridDim.x;
// now calculate the value at that position
int juliaValue = julia( x, y );
ptr[offset*4 + 0] = (juliaValue>>4)+(juliaValue%16<<4);
ptr[offset*4 + 1] = (juliaValue>>3)+(juliaValue%8<<5);
ptr[offset*4 + 2] = (juliaValue>>5)+(juliaValue%32<<3);
ptr[offset*4 + 3] = 255;
}
// globals needed by the update routine
struct DataBlock {
unsigned char *dev_bitmap;
};
int main( void ) {
DataBlock data;
CPUBitmap bitmap( DIM, DIM, &data );
unsigned char *dev_bitmap;
HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );
data.dev_bitmap = dev_bitmap;
dim3 grid(DIM,DIM);
kernel<<<grid,1>>>( dev_bitmap );
HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
bitmap.image_size(),
cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaFree( dev_bitmap ) );
bitmap.display_and_exit();
}
--
Jan Jacobs
Esdoornstraat 33
5995AN Kessel
W: www.sommaps.com
T: +31 77 462 1887
M: +31 6 23 82 55 21
E: [email protected]
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
