First commit [09/03/2018]
This commit is contained in:
Pierre8rTeam
@@ -0,0 +1,217 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort.mq5 |
|
||||
//| Copyright 2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
//--- COpenCL class
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
#resource "Kernels/bitonicsort.cl" as string cl_program
|
||||
//+------------------------------------------------------------------+
|
||||
//| QuickSortAscending |
|
||||
//+------------------------------------------------------------------+
|
||||
//| The function sorts array[] QuickSort algorithm. |
|
||||
//| |
|
||||
//| Arguments: |
|
||||
//| array : Array with values to sort |
|
||||
//| first : First element index |
|
||||
//| last : Last element index |
|
||||
//| |
|
||||
//| Return value: None |
|
||||
//+------------------------------------------------------------------+
|
||||
void QuickSortAscending(double &array[],int first,int last)
|
||||
{
|
||||
int i,j;
|
||||
double p_double,t_double;
|
||||
if(first<0 || last<0)
|
||||
return;
|
||||
i=first;
|
||||
j=last;
|
||||
while(i<last)
|
||||
{
|
||||
p_double=array[(first+last)>>1];
|
||||
while(i<j)
|
||||
{
|
||||
while(array[i]<p_double)
|
||||
{
|
||||
if(i==ArraySize(array)-1)
|
||||
break;
|
||||
i++;
|
||||
}
|
||||
while(array[j]>p_double)
|
||||
{
|
||||
if(j==0)
|
||||
break;
|
||||
j--;
|
||||
}
|
||||
if(i<=j)
|
||||
{
|
||||
//-- swap elements i and j
|
||||
t_double=array[i];
|
||||
array[i]=array[j];
|
||||
array[j]=t_double;
|
||||
i++;
|
||||
if(j==0)
|
||||
break;
|
||||
j--;
|
||||
}
|
||||
}
|
||||
if(first<j)
|
||||
QuickSortAscending(array,first,j);
|
||||
first=i;
|
||||
j=last;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| QuickSort_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool QuickSort_CPU(double &data_array[],ulong &time_cpu)
|
||||
{
|
||||
int data_count=ArraySize(data_array);
|
||||
if(data_count<=1)
|
||||
return(false);
|
||||
//--- sort values on CPU
|
||||
time_cpu=GetMicrosecondCount();
|
||||
QuickSortAscending(data_array,0,data_count-1);
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu)/1000);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool BitonicSort_GPU(COpenCL &OpenCL,double &data_array[],ulong &time_gpu)
|
||||
{
|
||||
int data_count=ArraySize(data_array);
|
||||
if(data_count<=1)
|
||||
return(false);
|
||||
//--- check support working with double
|
||||
if(!OpenCL.SupportDouble())
|
||||
{
|
||||
PrintFormat("Working with double (cl_khr_fp64) is not supported on the device.");
|
||||
return(false);
|
||||
}
|
||||
OpenCL.SetKernelsCount(1);
|
||||
OpenCL.KernelCreate(0,"BitonicSort_GPU");
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(1);
|
||||
if(!OpenCL.BufferFromArray(0,data_array,0,data_count,CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
OpenCL.SetArgumentBuffer(0,0,0);
|
||||
//---
|
||||
uint work_offset[1]={0};
|
||||
uint global_size[1];
|
||||
global_size[0]=data_count>>1;
|
||||
//---
|
||||
uint passes_total=0;
|
||||
uint stages_total=0;
|
||||
//---
|
||||
for(uint temp=data_count; temp>1; temp>>=1)
|
||||
stages_total++;
|
||||
//--- GPU calculation start
|
||||
time_gpu=GetMicrosecondCount();
|
||||
for(uint stage=0; stage<stages_total; stage++)
|
||||
{
|
||||
//--- set stage of the algorithm
|
||||
OpenCL.SetArgument(0,1,stage);
|
||||
for(uint pass=0; pass<stage+1; pass++)
|
||||
{
|
||||
//--- set pass of the current stage
|
||||
OpenCL.SetArgument(0,2,pass);
|
||||
//--- execute kernel
|
||||
if(!OpenCL.Execute(0,1,work_offset,global_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
else
|
||||
passes_total++;
|
||||
}
|
||||
}
|
||||
//---
|
||||
if(!OpenCL.BufferRead(0,data_array,0,0,data_count))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for data array C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finish
|
||||
time_gpu=ulong((GetMicrosecondCount()-time_gpu)/1000);
|
||||
PrintFormat("Bitonic sort finished. Total stages=%d, total passes=%d",stages_total,passes_total);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareDataArray |
|
||||
//+------------------------------------------------------------------+
|
||||
bool PrepareDataArray(long global_memory_size,double &data[],int &data_count)
|
||||
{
|
||||
int pwr_max=(int)(MathLog(global_memory_size/sizeof(float))/MathLog(2));
|
||||
int pwr=(int)MathMax(15,pwr_max-4);
|
||||
data_count=(int)MathPow(2,pwr);
|
||||
//--- prepare array and generate random data
|
||||
if(ArrayResize(data,data_count)<data_count)
|
||||
{
|
||||
Print("Error in ArrayResize. Error code=",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
for(int i=0; i<data_count; i++)
|
||||
data[i]=(double)(100000000*MathRand()/32767.0);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//--- OpenCL
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
long global_memory_size=0;
|
||||
if(!OpenCL.GetGlobalMemorySize(global_memory_size))
|
||||
{
|
||||
Print("Error in request of global memory size. Error code=",GetLastError());
|
||||
return;
|
||||
}
|
||||
double data_cpu[];
|
||||
int data_count;
|
||||
//--- prepare array with random values
|
||||
if(PrepareDataArray(global_memory_size,data_cpu,data_count)==false)
|
||||
return;
|
||||
//--- copy array values for sorting on GPU
|
||||
double data_gpu[];
|
||||
if(ArrayCopy(data_gpu,data_cpu,0,0,data_count)!=data_count)
|
||||
return;
|
||||
//--- Quick sort values using CPU
|
||||
ulong time_cpu=0;
|
||||
if(!QuickSort_CPU(data_cpu,time_cpu))
|
||||
return;
|
||||
//--- Bitonic sort values using GPU
|
||||
ulong time_gpu=0;
|
||||
if(!BitonicSort_GPU(OpenCL,data_gpu,time_gpu))
|
||||
return;
|
||||
//--- remove OpenCL objects
|
||||
OpenCL.Shutdown();
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
PrintFormat("time CPU=%d ms, time GPU =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- check calculations
|
||||
double total_error=0;
|
||||
for(int i=0; i<data_count; i++)
|
||||
{
|
||||
total_error+=MathAbs(data_gpu[i]-data_cpu[i]);
|
||||
}
|
||||
PrintFormat("Total error = %f",total_error);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,311 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT.mq5 |
|
||||
//| Copyright 2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
|
||||
#include <Math/Stat/Math.mqh>
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
|
||||
#resource "Kernels/fft.cl" as string cl_program
|
||||
#define kernel_init "fft_init"
|
||||
#define kernel_stage "fft_stage"
|
||||
#define kernel_scale "fft_scale"
|
||||
#define NUM_POINTS 16384
|
||||
#define FFT_DIRECTION 1
|
||||
//+------------------------------------------------------------------+
|
||||
//| Fast Fourier transform and its inverse (both recursively) |
|
||||
//| Copyright (C) 2004, Jerome R. Breitenbach. All rights reserved. |
|
||||
//| Reference: |
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Recursive direct FFT transform |
|
||||
//+------------------------------------------------------------------+
|
||||
void fft(const int N,double &x_real[],double &x_imag[],double &X_real[],double &X_imag[])
|
||||
{
|
||||
//--- prepare temporary arrays
|
||||
double XX_real[],XX_imag[];
|
||||
ArrayResize(XX_real,N);
|
||||
ArrayResize(XX_imag,N);
|
||||
//--- calculate FFT by a recursion
|
||||
fft_rec(N,0,1,x_real,x_imag,X_real,X_imag,XX_real,XX_imag);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Recursive inverse FFT transform |
|
||||
//+------------------------------------------------------------------+
|
||||
void ifft(const int N,double &x_real[],double &x_imag[],double &X_real[],double &X_imag[])
|
||||
{
|
||||
int N2=N/2; // half the number of points in IFFT
|
||||
//--- calculate IFFT via reciprocity property of DFT
|
||||
fft(N,X_real,X_imag,x_real,x_imag);
|
||||
x_real[0]=x_real[0]/N;
|
||||
x_imag[0]=x_imag[0]/N;
|
||||
x_real[N2]=x_real[N2]/N;
|
||||
x_imag[N2]=x_imag[N2]/N;
|
||||
for(int i=1; i<N2; i++)
|
||||
{
|
||||
double tmp0=x_real[i]/N;
|
||||
double tmp1=x_imag[i]/N;
|
||||
x_real[i]=x_real[N-i]/N;
|
||||
x_imag[i]=x_imag[N-i]/N;
|
||||
x_real[N-i]=tmp0;
|
||||
x_imag[N-i]=tmp1;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT recursion |
|
||||
//+------------------------------------------------------------------+
|
||||
void fft_rec(const int N,const int offset,const int delta,double &x_real[],double &x_imag[],double &X_real[],double &X_imag[],double &XX_real[],double &XX_imag[])
|
||||
{
|
||||
static const double TWO_PI=(double)(2*M_PI);
|
||||
int N2=N/2; // half the number of points in FFT
|
||||
int k00,k01,k10,k11; // indices for butterflies
|
||||
if(N!=2)
|
||||
{
|
||||
//--- perform recursive step
|
||||
//--- calculate two (N/2)-point DFT's
|
||||
fft_rec(N2,offset,2*delta,x_real,x_imag,XX_real,XX_imag,X_real,X_imag);
|
||||
fft_rec(N2,offset+delta,2*delta,x_real,x_imag,XX_real,XX_imag,X_real,X_imag);
|
||||
//--- combine the two (N/2)-point DFT's into one N-point DFT
|
||||
for(int k=0; k<N2; k++)
|
||||
{
|
||||
k00 = offset + k*delta;
|
||||
k01 = k00 + N2*delta;
|
||||
k10 = offset + 2*k*delta;
|
||||
k11 = k10 + delta;
|
||||
double cs=(double)MathCos(TWO_PI*k/(double)N);
|
||||
double sn=(double)MathSin(TWO_PI*k/(double)N);
|
||||
double tmp0 = cs*XX_real[k11] + sn*XX_imag[k11];
|
||||
double tmp1 = cs*XX_imag[k11] - sn*XX_real[k11];
|
||||
X_real[k01] = XX_real[k10] - tmp0;
|
||||
X_imag[k01] = XX_imag[k10] - tmp1;
|
||||
X_real[k00] = XX_real[k10] + tmp0;
|
||||
X_imag[k00] = XX_imag[k10] + tmp1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
//--- perform 2-point DFT
|
||||
k00=offset;
|
||||
k01=k00+delta;
|
||||
X_real[k01] = x_real[k00] - x_real[k01];
|
||||
X_imag[k01] = x_imag[k00] - x_imag[k01];
|
||||
X_real[k00] = x_real[k00] + x_real[k01];
|
||||
X_imag[k00] = x_imag[k00] + x_imag[k01];
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool FFT_CPU(int direction,int power,double &data_real[],double &data_imag[],ulong &time_cpu)
|
||||
{
|
||||
//--- calculate the number of points
|
||||
int N=1;
|
||||
for(int i=0;i<power;i++)
|
||||
N*=2;
|
||||
//---prepare temporary arrays
|
||||
double XX_real[],XX_imag[];
|
||||
ArrayResize(XX_real,N);
|
||||
ArrayResize(XX_imag,N);
|
||||
//--- CPU calculation start
|
||||
time_cpu=GetMicrosecondCount();
|
||||
if(direction>0)
|
||||
fft(N,data_real,data_imag,XX_real,XX_imag);
|
||||
else
|
||||
ifft(N,XX_real,XX_imag,data_real,data_imag);
|
||||
//--- CPU calculation finished
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu));
|
||||
//--- copy calculated data
|
||||
ArrayCopy(data_real,XX_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(data_imag,XX_imag,0,0,WHOLE_ARRAY);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool FFT_GPU(int direction,int power,double &data_real[],double &data_imag[],ulong &time_gpu)
|
||||
{
|
||||
//--- calculate the number of points
|
||||
int num_points=1;
|
||||
for(int i=0;i<power;i++)
|
||||
num_points*=2;
|
||||
//--- prepare data array for GPU calculation
|
||||
double data[];
|
||||
ArrayResize(data,2*num_points);
|
||||
for(int i=0; i<num_points; i++)
|
||||
{
|
||||
data[2*i]=data_real[i];
|
||||
data[2*i+1]=data_imag[i];
|
||||
}
|
||||
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- check support working with double
|
||||
if(!OpenCL.SupportDouble())
|
||||
{
|
||||
PrintFormat("Working with double (cl_khr_fp64) is not supported on the device.");
|
||||
return(false);
|
||||
}
|
||||
//--- create kernels
|
||||
OpenCL.SetKernelsCount(3);
|
||||
OpenCL.KernelCreate(0,kernel_init);
|
||||
OpenCL.KernelCreate(1,kernel_stage);
|
||||
OpenCL.KernelCreate(2,kernel_scale);
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(2);
|
||||
if(!OpenCL.BufferFromArray(0,data,0,2*num_points,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for input buffer. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferCreate(1,2*num_points*sizeof(double),CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for data buffer. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- determine maximum work-group size
|
||||
int local_size=(int)CLGetInfoInteger(OpenCL.GetKernel(0),CL_KERNEL_WORK_GROUP_SIZE);
|
||||
//--- determine local memory size
|
||||
uint local_mem_size=(uint)CLGetInfoInteger(OpenCL.GetContext(),CL_DEVICE_LOCAL_MEM_SIZE);
|
||||
int points_per_group=(int)local_mem_size/(2*sizeof(double));
|
||||
if(points_per_group>num_points)
|
||||
points_per_group=num_points;
|
||||
//--- set kernel arguments
|
||||
OpenCL.SetArgumentBuffer(0,0,0);
|
||||
OpenCL.SetArgumentBuffer(0,1,1);
|
||||
OpenCL.SetArgumentLocalMemory(0,2,local_mem_size);
|
||||
OpenCL.SetArgument(0,3,points_per_group);
|
||||
OpenCL.SetArgument(0,4,num_points);
|
||||
OpenCL.SetArgument(0,5,direction);
|
||||
//--- OpenCL execute settings
|
||||
int task_dimension=1;
|
||||
uint global_size=(uint)((num_points/points_per_group)*local_size);
|
||||
uint global_work_offset[1]={0};
|
||||
uint global_work_size[1];
|
||||
global_work_size[0]=global_size;
|
||||
uint local_work_size[1];
|
||||
local_work_size[0]=local_size;
|
||||
//--- GPU calculation start
|
||||
time_gpu=GetMicrosecondCount();
|
||||
//-- execute kernel fft_init
|
||||
if(!OpenCL.Execute(0,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_init: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//-- further stages of the FFT
|
||||
if(num_points>points_per_group)
|
||||
{
|
||||
//--- set arguments for kernel 1
|
||||
OpenCL.SetArgumentBuffer(1,0,1);
|
||||
OpenCL.SetArgument(1,2,points_per_group);
|
||||
OpenCL.SetArgument(1,3,direction);
|
||||
for(int stage=2; stage<=num_points/points_per_group; stage<<=1)
|
||||
{
|
||||
OpenCL.SetArgument(1,1,stage);
|
||||
//-- execute kernel fft_stage
|
||||
if(!OpenCL.Execute(1,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_stage: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
}
|
||||
}
|
||||
//--- scale values if performing the inverse FFT
|
||||
if(direction<0)
|
||||
{
|
||||
OpenCL.SetArgumentBuffer(2,0,1);
|
||||
OpenCL.SetArgument(2,1,points_per_group);
|
||||
OpenCL.SetArgument(2,2,num_points);
|
||||
//-- execute kernel fft_scale
|
||||
if(!OpenCL.Execute(2,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_scale: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
}
|
||||
//--- read the results from GPU memory
|
||||
if(!OpenCL.BufferRead(1,data,0,0,2*num_points))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for data_buffer2. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time_gpu=ulong((GetMicrosecondCount()-time_gpu));
|
||||
//--- copy calculated data and release OpenCL handles
|
||||
for(int i=0; i<num_points; i++)
|
||||
{
|
||||
data_real[i]=data[2*i];
|
||||
data_imag[i]=data[2*i+1];
|
||||
}
|
||||
OpenCL.Shutdown();
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
int datacount=NUM_POINTS;
|
||||
int power=(int)(MathLog(NUM_POINTS)/M_LN2);
|
||||
if(MathPow(2,power)!=datacount)
|
||||
{
|
||||
PrintFormat("Number of elements must be power of 2. Elements: %d",datacount);
|
||||
return;
|
||||
}
|
||||
//--- prepare data for FFT calculation
|
||||
double data_real[],data_imag[];
|
||||
ArrayResize(data_real,datacount);
|
||||
ArrayResize(data_imag,datacount);
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
data_real[i]=(double)i;
|
||||
data_imag[i]=0;
|
||||
}
|
||||
int direction=FFT_DIRECTION;
|
||||
//--- data arrays for CPU calculation
|
||||
double CPU_real[],CPU_imag[];
|
||||
ArrayCopy(CPU_real,data_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(CPU_imag,data_imag,0,0,WHOLE_ARRAY);
|
||||
ulong time_cpu=0;
|
||||
//--- calculate FFT using CPU
|
||||
FFT_CPU(direction,power,CPU_real,CPU_imag,time_cpu);
|
||||
|
||||
//--- data arrays for GPU calculation
|
||||
double GPU_real[],GPU_imag[];
|
||||
ArrayCopy(GPU_real,data_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(GPU_imag,data_imag,0,0,WHOLE_ARRAY);
|
||||
ulong time_gpu=0;
|
||||
//--- calculate FFT using GPU
|
||||
if(!FFT_GPU(direction,power,GPU_real,GPU_imag,time_gpu))
|
||||
{
|
||||
PrintFormat("Error in calculation FFT on GPU.");
|
||||
return;
|
||||
}
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
PrintFormat("FFT calculation for %d points.",datacount);
|
||||
PrintFormat("time CPU=%d microseconds, time GPU =%d microseconds, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- determine average error
|
||||
double average_error=0.0;
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
average_error += MathAbs(CPU_real[i]-GPU_real[i]);
|
||||
average_error += MathAbs(CPU_imag[i]-GPU_imag[i]);
|
||||
}
|
||||
average_error=average_error/(datacount*2);
|
||||
PrintFormat("Average error = %f",average_error);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,31 @@
|
||||
//--- by default some GPU doesn't support doubles
|
||||
//--- cl_khr_fp64 directive is used to enable work with doubles
|
||||
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel |
|
||||
//| The bitonic sort kernel does an ascending sort. |
|
||||
//+-----------------------------------------------------------+
|
||||
//| R. Banger,K. Bhattacharyya, OpenCL Programming by Example:|
|
||||
//| A comprehensive guide on OpenCL programming with examples |
|
||||
//| PACKT Publishing, 2013. |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void BitonicSort_GPU(__global double *data,const uint stage,const uint pass)
|
||||
{
|
||||
uint id=get_global_id(0);
|
||||
uint distance = 1<<(stage-pass);
|
||||
uint left_id =(id &(distance-1));
|
||||
left_id+=(id>>(stage-pass))*(distance<<1);
|
||||
uint right_id=left_id+distance;
|
||||
double left_value=data[left_id];
|
||||
double right_value=data[right_id];
|
||||
uint same_direction=(id>>stage)&0x1;
|
||||
uint temp = same_direction?right_id:temp;
|
||||
right_id = same_direction?left_id:right_id;
|
||||
left_id = same_direction?temp:left_id;
|
||||
int compare_res=(left_value<right_value);
|
||||
double greater = compare_res?right_value:left_value;
|
||||
double lesser = compare_res?left_value:right_value;
|
||||
data[left_id] = lesser;
|
||||
data[right_id]= greater;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,156 @@
|
||||
//--- by default some GPU doesn't support doubles
|
||||
//--- cl_khr_fp64 directive is used to enable work with doubles
|
||||
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_init OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_init(__global double2 *in_data,
|
||||
__global double2 *out_data,
|
||||
__local double2 *l_data,
|
||||
uint points_per_group,uint size,int dir)
|
||||
{
|
||||
uint4 br,index;
|
||||
uint points_per_item,g_addr,l_addr,i,fft_index,stage,N2;
|
||||
double2 x1,x2,x3,x4,sum12,diff12,sum34,diff34;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
l_addr = get_local_id(0)*points_per_item;
|
||||
g_addr = get_group_id(0)*points_per_group + l_addr;
|
||||
//--- load data from bit-reversed addresses and perform 4-point FFTs
|
||||
for(i=0; i<points_per_item; i+=4)
|
||||
{
|
||||
index=(uint4)(g_addr,g_addr+1,g_addr+2,g_addr+3);
|
||||
fft_index=size/2;
|
||||
stage=1;
|
||||
N2 =(uint)log2((double)size)-1;
|
||||
br =(index<< N2) & fft_index;
|
||||
br|=(index>> N2) & stage;
|
||||
//--- bit-reverse addresses
|
||||
while(N2>1)
|
||||
{
|
||||
N2-=2;
|
||||
fft_index>>=1;
|
||||
stage<<=1;
|
||||
br |= (index << N2) & fft_index;
|
||||
br |= (index >> N2) & stage;
|
||||
}
|
||||
//--- load global data
|
||||
x1 = in_data[br.s0];
|
||||
x2 = in_data[br.s1];
|
||||
x3 = in_data[br.s2];
|
||||
x4 = in_data[br.s3];
|
||||
|
||||
sum12=x1+x2;
|
||||
diff12= x1-x2;
|
||||
sum34 = x3+x4;
|
||||
diff34=(double2)(x3.s1-x4.s1,x4.s0-x3.s0)*dir;
|
||||
l_data[l_addr]=sum12+sum34;
|
||||
l_data[l_addr+1] = diff12 + diff34;
|
||||
l_data[l_addr+2] = sum12 - sum34;
|
||||
l_data[l_addr+3] = diff12 - diff34;
|
||||
l_addr += 4;
|
||||
g_addr += 4;
|
||||
}
|
||||
//--- perform initial stages of the FFT - each of length N2*2
|
||||
for(N2=4; N2<points_per_item; N2<<=1)
|
||||
{
|
||||
l_addr=get_local_id(0)*points_per_item;
|
||||
for(fft_index=0; fft_index<points_per_item; fft_index+=2*N2)
|
||||
{
|
||||
x1=l_data[l_addr];
|
||||
l_data[l_addr]+=l_data[l_addr+N2];
|
||||
l_data[l_addr+N2]=x1-l_data[l_addr+N2];
|
||||
for(i=1; i<N2; i++)
|
||||
{
|
||||
x3.s0=cos(M_PI_F*i/N2);
|
||||
x3.s1=dir*sin(M_PI_F*i/N2);
|
||||
x2=(double2)(l_data[l_addr+N2+i].s0*x3.s0+l_data[l_addr+N2+i].s1*x3.s1,l_data[l_addr+N2+i].s1*x3.s0-l_data[l_addr+N2+i].s0*x3.s1);
|
||||
l_data[l_addr+N2+i]=l_data[l_addr+i]-x2;
|
||||
l_data[l_addr+i]+=x2;
|
||||
}
|
||||
l_addr+=2*N2;
|
||||
}
|
||||
}
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
//--- perform FFT with other items in group - each of length N2*2
|
||||
stage=2;
|
||||
for(N2=points_per_item; N2<points_per_group; N2<<=1)
|
||||
{
|
||||
br.s0=(get_local_id(0)+(get_local_id(0)/stage)*stage) *(points_per_item/2);
|
||||
size = br.s0 % (N2*2);
|
||||
for(i=br.s0; i<br.s0+points_per_item/2; i++)
|
||||
{
|
||||
x3.s0=cos(M_PI_F*size/N2);
|
||||
x3.s1=dir*sin(M_PI_F*size/N2);
|
||||
x2=(double2)(l_data[N2+i].s0*x3.s0+l_data[N2+i].s1*x3.s1,l_data[N2+i].s1*x3.s0-l_data[N2+i].s0*x3.s1);
|
||||
l_data[N2+i]=l_data[i]-x2;
|
||||
l_data[i]+=x2;
|
||||
size++;
|
||||
}
|
||||
stage<<=1;
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
}
|
||||
//--- store results in global memory
|
||||
l_addr = get_local_id(0)*points_per_item;
|
||||
g_addr = get_group_id(0)*points_per_group + l_addr;
|
||||
for(i=0; i<points_per_item; i+=4)
|
||||
{
|
||||
out_data[g_addr]=l_data[l_addr];
|
||||
out_data[g_addr+1] = l_data[l_addr+1];
|
||||
out_data[g_addr+2] = l_data[l_addr+2];
|
||||
out_data[g_addr+3] = l_data[l_addr+3];
|
||||
g_addr += 4;
|
||||
l_addr += 4;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_stage OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_stage(__global double2 *g_data,uint stage,uint points_per_group,int dir)
|
||||
{
|
||||
uint points_per_item,addr,N,ang,i;
|
||||
double c,s;
|
||||
double2 input1,input2,w;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
addr=(get_group_id(0)+(get_group_id(0)/stage)*stage)*(points_per_group/2)+get_local_id(0)*(points_per_item/2);
|
||||
N=points_per_group*(stage/2);
|
||||
ang=addr%(N*2);
|
||||
|
||||
for(i=addr; i<addr+points_per_item/2; i++)
|
||||
{
|
||||
c = cos(M_PI_F*ang/N);
|
||||
s = dir*sin(M_PI_F*ang/N);
|
||||
input1 = g_data[i];
|
||||
input2 = g_data[i+N];
|
||||
w=(double2)(input2.s0*c+input2.s1*s,input2.s1*c-input2.s0*s);
|
||||
g_data[i]=input1+w;
|
||||
g_data[i+N]=input1-w;
|
||||
ang++;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_scale OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_scale(__global double2 *g_data,uint points_per_group,uint scale)
|
||||
{
|
||||
uint points_per_item,addr,i;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
addr=get_group_id(0)*points_per_group+get_local_id(0)*points_per_item;
|
||||
|
||||
for(i=addr; i<addr+points_per_item; i++)
|
||||
{
|
||||
g_data[i]/=scale;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,69 @@
|
||||
//--- by default some GPU doesn't support doubles
|
||||
//--- cl_khr_fp64 directive is used to enable work with doubles
|
||||
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel for matrix multiplication |
|
||||
//| using global work groups |
|
||||
//+-----------------------------------------------------------+
|
||||
//| http://gpgpu-computing4.blogspot.ru/2009/09/ |
|
||||
//| /matrix-multiplication-2-opencl.html |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void MatrixMult_GPU1(__global double *matrix_a,
|
||||
__global double *matrix_b,
|
||||
__global double *matrix_c,
|
||||
int rows_a,int cols_a,int cols_b)
|
||||
{
|
||||
int i=get_global_id(0);
|
||||
int j=get_global_id(1);
|
||||
double sum=0.0;
|
||||
for(int k=0; k<cols_a; k++)
|
||||
{
|
||||
sum+=matrix_a[cols_a*i+k]*matrix_b[cols_b*k+j];
|
||||
}
|
||||
matrix_c[cols_b*i+j]=sum;
|
||||
}
|
||||
#define BLOCK_SIZE 10
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel for matrix multiplication |
|
||||
//| using local groups with common local memory |
|
||||
//+-----------------------------------------------------------+
|
||||
//| http://gpgpu-computing4.blogspot.ru/2009/10/ |
|
||||
//| /matrix-multiplication-3-opencl.html |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void MatrixMult_GPU2(__global double *matrix_a,
|
||||
__global double *matrix_b,
|
||||
__global double *matrix_c,
|
||||
int rows_a,int cols_a,int cols_b)
|
||||
{
|
||||
int group_i=get_group_id(0);
|
||||
int group_j=get_group_id(1);
|
||||
|
||||
int i=get_local_id(0);
|
||||
int j=get_local_id(1);
|
||||
|
||||
int offset_b=BLOCK_SIZE*group_i;
|
||||
int offset_a_start=cols_a*BLOCK_SIZE*group_j;
|
||||
double sum=(float)0.0;
|
||||
|
||||
for(int offset_a=offset_a_start;
|
||||
offset_a<offset_a_start+cols_a;
|
||||
offset_a+=BLOCK_SIZE,
|
||||
offset_b+=BLOCK_SIZE*cols_b)
|
||||
{
|
||||
__local double submatrix_a[BLOCK_SIZE][BLOCK_SIZE];
|
||||
__local double submatrix_b[BLOCK_SIZE][BLOCK_SIZE];
|
||||
|
||||
submatrix_a[i][j]=matrix_a[offset_a+cols_a*i+j];
|
||||
submatrix_b[i][j]=matrix_b[offset_b+cols_b*i+j];
|
||||
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
|
||||
for(int k=0; k<BLOCK_SIZE; k++)
|
||||
sum+=submatrix_a[i][k]*submatrix_b[k][j];
|
||||
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
}
|
||||
int offset_c=BLOCK_SIZE*(cols_b*group_j+group_i);
|
||||
matrix_c[offset_c+cols_b*i+j]=sum;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,51 @@
|
||||
//--- by default some GPU doesn't support doubles
|
||||
//--- cl_khr_fp64 directive is used to enable work with doubles
|
||||
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
||||
//+------------------------------------------------------------------+
|
||||
//| Morlet wavelet function |
|
||||
//+------------------------------------------------------------------+
|
||||
double Morlet(const double t)
|
||||
{
|
||||
return exp(-t*t*0.5)*cos(M_2_PI*t);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| OpenCL kernel function |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void Wavelet_GPU(__global double *data,int datacount,int x_size,int y_size,__global double *result)
|
||||
{
|
||||
size_t i = get_global_id(0);
|
||||
size_t j = get_global_id(1);
|
||||
double a1=(double)10e-10;
|
||||
double a2=(double)15.0;
|
||||
double da=(a2-a1)/(double)y_size;
|
||||
double db=((double)datacount-(double)0.0)/x_size;
|
||||
double a=a1+j*da;
|
||||
double b=0+i*db;
|
||||
uint norm=1;
|
||||
double B=(double)1.0; //Morlet
|
||||
double B_inv=(double)1.0/B;
|
||||
double a_inv=(double)1.0/a;
|
||||
double dt=(double)1.0;
|
||||
double coef=(double)0.0;
|
||||
if(norm==0)
|
||||
coef=sqrt(a_inv);
|
||||
else
|
||||
{
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
double arg=(dt*k-b)*a_inv;
|
||||
arg=-B_inv*arg*arg;
|
||||
coef=coef+exp(arg);
|
||||
}
|
||||
}
|
||||
double sum=(float)0.0;
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
double arg=(dt*k-b)*a_inv;
|
||||
sum+=data[k]*Morlet(arg);
|
||||
}
|
||||
sum=sum/coef;
|
||||
uint pos=(int)(j*x_size+i);
|
||||
result[pos]=sum;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,225 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult.mq5 |
|
||||
//| Copyright 2016-2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2016-2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
//--- OpenCL kernels
|
||||
#resource "Kernels/matrixmult.cl" as string cl_program
|
||||
#define BLOCK_SIZE 10
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool MatrixMult_CPU(const double &matrix_a[],const double &matrix_b[],double &matrix_c[],
|
||||
const int rows_a,const int cols_a,const int cols_b,ulong &time_cpu)
|
||||
{
|
||||
int size=rows_a*cols_b;
|
||||
if(ArrayResize(matrix_c,size)!=size)
|
||||
return(false);
|
||||
//--- CPU calculation started
|
||||
time_cpu=GetMicrosecondCount();
|
||||
for(int i=0; i<rows_a; i++)
|
||||
{
|
||||
for(int j=0; j<cols_b; j++)
|
||||
{
|
||||
double sum=0.0;
|
||||
for(int k=0; k<cols_a; k++)
|
||||
{
|
||||
sum+=matrix_a[cols_a*i+k]*matrix_b[cols_b*k+j];
|
||||
}
|
||||
matrix_c[cols_b*i+j]=sum;
|
||||
}
|
||||
}
|
||||
//--- CPU calculation finished
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu)/1000);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool MatrixMult_GPU(const double &matrix_a[],const double &matrix_b[],double &matrix1_c[],double &matrix2_c[],
|
||||
const int rows_a,const int cols_a,const int cols_b,const int size_a,const int size_b,
|
||||
const int size_c,ulong &time1_gpu,ulong &time2_gpu)
|
||||
{
|
||||
const int task_dimension=2;
|
||||
//--- prepare matrices for result
|
||||
if(ArrayResize(matrix1_c,size_c)!=size_c || ArrayResize(matrix2_c,size_c)!=size_c)
|
||||
return(false);
|
||||
ArrayFill(matrix1_c,0,size_c,(double)0.0);
|
||||
ArrayFill(matrix2_c,0,size_c,(double)0.0);
|
||||
//--- OpenCL
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- check support working with double
|
||||
if(!OpenCL.SupportDouble())
|
||||
{
|
||||
PrintFormat("Working with double (cl_khr_fp64) is not supported on the device.");
|
||||
return(false);
|
||||
}
|
||||
//--- create kernels
|
||||
OpenCL.SetKernelsCount(2);
|
||||
OpenCL.KernelCreate(0,"MatrixMult_GPU1");
|
||||
OpenCL.KernelCreate(1,"MatrixMult_GPU2");
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(3);
|
||||
//---
|
||||
if(!OpenCL.BufferFromArray(0,matrix_a,0,size_a,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for matrix A. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferFromArray(1,matrix_b,0,size_b,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for matrix B. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferCreate(2,size_c*sizeof(double),CL_MEM_WRITE_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for matrix C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- prepare arguments for kernel 0
|
||||
int kernel_index=0;
|
||||
OpenCL.SetArgumentBuffer(kernel_index,0,0);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,1,1);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,2,2);
|
||||
OpenCL.SetArgument(kernel_index,3,rows_a);
|
||||
OpenCL.SetArgument(kernel_index,4,cols_a);
|
||||
OpenCL.SetArgument(kernel_index,5,cols_b);
|
||||
//--- set task dimension a_rows x b_cols
|
||||
uint global_work_size[2];
|
||||
//--- set dimensions
|
||||
global_work_size[0]=rows_a;
|
||||
global_work_size[1]=cols_b;
|
||||
uint global_work_offset[2]={0,0};
|
||||
//--- GPU calculation start kernel 0
|
||||
time1_gpu=GetMicrosecondCount();
|
||||
if(!OpenCL.Execute(kernel_index,task_dimension,global_work_offset,global_work_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferRead(2,matrix1_c,0,0,size_c))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for matrix1 C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time1_gpu=ulong((GetMicrosecondCount()-time1_gpu)/1000);
|
||||
|
||||
//--- prepare arguments for kernel 1
|
||||
kernel_index=1;
|
||||
//--- set arguments
|
||||
OpenCL.SetArgumentBuffer(kernel_index,0,0);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,1,1);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,2,2);
|
||||
OpenCL.SetArgument(kernel_index,3,rows_a);
|
||||
OpenCL.SetArgument(kernel_index,4,cols_a);
|
||||
OpenCL.SetArgument(kernel_index,5,cols_b);
|
||||
uint local_work_size[2];
|
||||
local_work_size[0]=BLOCK_SIZE;
|
||||
local_work_size[1]=BLOCK_SIZE;
|
||||
//--- GPU calculation start, kernel1
|
||||
time2_gpu=GetMicrosecondCount();
|
||||
if(!OpenCL.Execute(kernel_index,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferRead(2,matrix2_c,0,0,size_c))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for matrix2 C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time2_gpu=ulong((GetMicrosecondCount()-time2_gpu)/1000);
|
||||
//--- remove OpenCL objects
|
||||
OpenCL.Shutdown();
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//--- matrix A 1000x2000
|
||||
int rows_a=1000;
|
||||
int cols_a=2000;
|
||||
//--- matrix B 2000x1000
|
||||
int rows_b=cols_a;
|
||||
int cols_b=1000;
|
||||
//--- matrix C 1000x1000
|
||||
int rows_c=rows_a;
|
||||
int cols_c=cols_b;
|
||||
//--- matrix A: size=rows_a*cols_a
|
||||
int size_a=rows_a*cols_a;
|
||||
int size_b=rows_b*cols_b;
|
||||
int size_c=rows_c*cols_c;
|
||||
//--- prepare matrix A
|
||||
double matrix_a[];
|
||||
ArrayResize(matrix_a,rows_a*cols_a);
|
||||
for(int i=0; i<rows_a; i++)
|
||||
for(int j=0; j<cols_a; j++)
|
||||
{
|
||||
matrix_a[i*cols_a+j]=(double)(10*MathRand()/32767);
|
||||
}
|
||||
//--- prepare matrix B
|
||||
double matrix_b[];
|
||||
ArrayResize(matrix_b,rows_b*cols_b);
|
||||
for(int i=0; i<rows_b; i++)
|
||||
for(int j=0; j<cols_b; j++)
|
||||
{
|
||||
matrix_b[i*cols_b+j]=(double)(10*MathRand()/32767);
|
||||
}
|
||||
//--- CPU: calculate matrix product matrix_a*matrix_b
|
||||
double matrix_c_cpu[];
|
||||
ulong time_cpu=0;
|
||||
if(!MatrixMult_CPU(matrix_a,matrix_b,matrix_c_cpu,rows_a,cols_a,cols_b,time_cpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on CPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//--- calculate matrix product using GPU
|
||||
double matrix_c_gpu_method1[];
|
||||
double matrix_c_gpu_method2[];
|
||||
ulong time_gpu_method1=0;
|
||||
ulong time_gpu_method2=0;
|
||||
if(!MatrixMult_GPU(matrix_a,matrix_b,matrix_c_gpu_method1,matrix_c_gpu_method2,rows_a,cols_a,cols_b,size_a,size_b,size_c,time_gpu_method1,time_gpu_method2))
|
||||
{
|
||||
PrintFormat("Error in calculation on GPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio1=0;
|
||||
double CPU_GPU_ratio2=0;
|
||||
if(time_gpu_method1!=0)
|
||||
CPU_GPU_ratio1=1.0*time_cpu/time_gpu_method1;
|
||||
if(time_gpu_method2!=0)
|
||||
CPU_GPU_ratio2=1.0*time_cpu/time_gpu_method2;
|
||||
PrintFormat("time CPU=%d ms, time GPU global work groups =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu_method1,CPU_GPU_ratio1);
|
||||
PrintFormat("time CPU=%d ms, time GPU local work groups =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu_method2,CPU_GPU_ratio2);
|
||||
//--- check calculations
|
||||
double total_error1=0;
|
||||
double total_error2=0;
|
||||
for(int i=0; i<rows_c; i++)
|
||||
{
|
||||
for(int j=0; j<cols_c; j++)
|
||||
{
|
||||
int pos=cols_c*i+j;
|
||||
total_error1+=MathAbs(matrix_c_gpu_method1[pos]-matrix_c_cpu[pos]);
|
||||
total_error2+=MathAbs(matrix_c_gpu_method2[pos]-matrix_c_cpu[pos]);
|
||||
}
|
||||
}
|
||||
PrintFormat("Total error for method 1 = %f",total_error1);
|
||||
PrintFormat("Total error for method 2 = %f",total_error2);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,419 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Wavelet.mq5 |
|
||||
//| Copyright 2016-2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2016-2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
|
||||
#include <Math/Stat/Math.mqh>
|
||||
#include <Graphics/Graphic.mqh>
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
|
||||
#define CPU_DATA 1
|
||||
#define GPU_DATA 2
|
||||
|
||||
#define SIZE_X 600
|
||||
#define SIZE_Y 200
|
||||
|
||||
#resource "Kernels/wavelet.cl" as string cl_program
|
||||
//+------------------------------------------------------------------+
|
||||
//| CWavelet |
|
||||
//+------------------------------------------------------------------+
|
||||
class CWavelet
|
||||
{
|
||||
protected:
|
||||
int m_xsize;
|
||||
int m_ysize;
|
||||
int m_maxcolor;
|
||||
string m_res_name;
|
||||
string m_label_name;
|
||||
uchar m_palette[3*256];
|
||||
//---
|
||||
double m_data[];
|
||||
double m_wavelet_data_CPU[];
|
||||
double m_wavelet_data_GPU[];
|
||||
uint m_bmp_buffer[];
|
||||
|
||||
COpenCL m_OpenCL;
|
||||
|
||||
double Morlet(const double t);
|
||||
void ShowWaveletData(const double &m_wavelet_data[]);
|
||||
int GetPalColor(const int index);
|
||||
void Blend(const uint c1,const uint c2,const uint r1,const uint g1,const uint b1,const uint r2,const uint g2,const uint b2);
|
||||
bool WaveletCPU(const double &data[],const int datacount,const int x_size,const int y_size,const int i,const int j,const bool norm,double &result[]);
|
||||
public:
|
||||
//---
|
||||
void Create(const string name,const int x0,const int y0,const int x_size,const int y_size);
|
||||
bool CalculateWavelet_CPU(const double &data[],uint &time);
|
||||
bool CalculateWavelet_GPU(double &data[],uint &time);
|
||||
void ShowWavelet(const int mode);
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
//| Morlet wavelet function |
|
||||
//+------------------------------------------------------------------+
|
||||
double CWavelet::Morlet(const double t)
|
||||
{
|
||||
double v=t;
|
||||
double res=MathExp(-v*v*0.5)*MathCos(M_2_PI*v);
|
||||
return ((double)res);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| GetPalColor |
|
||||
//+------------------------------------------------------------------+
|
||||
int CWavelet::GetPalColor(const int index)
|
||||
{
|
||||
int ind=index;
|
||||
if(ind<=0)
|
||||
ind=0;
|
||||
if(ind>255)
|
||||
ind=255;
|
||||
int idx=3*(ind);
|
||||
uchar r=m_palette[idx];
|
||||
uchar g=m_palette[idx+1];
|
||||
uchar b=m_palette[idx+2];
|
||||
//---
|
||||
return(b+256*g+65536*r);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Gradient palette |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::Blend(const uint c1,const uint c2,const uint r1,const uint g1,const uint b1,const uint r2,const uint g2,const uint b2)
|
||||
{
|
||||
int n=int(c2-c1);
|
||||
for(int i=0; i<=n; i++)
|
||||
{
|
||||
if((c1+i+2)<ArraySize(m_palette))
|
||||
{
|
||||
m_palette[3*(c1+i)]=uchar(MathRound(1*(r1*(n-i)+r2*i)*1.0/n));
|
||||
m_palette[3*(c1+i)+1]=uchar(MathRound(1*(g1*(n-i)+g2*i)*1.0/n));
|
||||
m_palette[3*(c1+i)+2]=uchar(MathRound(1*(b1*(n-i)+b2*i)*1.0/n));
|
||||
}
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Create |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::Create(const string name,const int x0,const int y0,const int x_size,const int y_size)
|
||||
{
|
||||
//---
|
||||
m_xsize=x_size;
|
||||
m_ysize=y_size;
|
||||
int size=m_xsize*m_ysize;
|
||||
ArrayResize(m_bmp_buffer,size);
|
||||
ArrayFill(m_bmp_buffer,0,size,0);
|
||||
ArrayResize(m_wavelet_data_CPU,size);
|
||||
ArrayResize(m_wavelet_data_GPU,size);
|
||||
ArrayFill(m_wavelet_data_CPU,0,size,0);
|
||||
ArrayFill(m_wavelet_data_GPU,0,size,0);
|
||||
m_res_name=name;
|
||||
m_label_name=m_res_name;
|
||||
StringToUpper(m_label_name);
|
||||
ResourceCreate(m_res_name,m_bmp_buffer,m_xsize,m_ysize,0,0,0,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
ObjectCreate(0,m_label_name,OBJ_BITMAP_LABEL,0,0,0);
|
||||
ObjectSetInteger(0,m_label_name,OBJPROP_XDISTANCE,x0);
|
||||
ObjectSetInteger(0,m_label_name,OBJPROP_YDISTANCE,y0);
|
||||
ObjectSetString(0,m_label_name,OBJPROP_BMPFILE,NULL);
|
||||
ObjectSetString(0,m_label_name,OBJPROP_BMPFILE,"::"+m_label_name);
|
||||
//---
|
||||
m_maxcolor=100;
|
||||
Blend(0,20,0,0,95,0,0,246);
|
||||
Blend(21,40,0,0,246,0,236,226);
|
||||
Blend(41,60,0,236,226,226,246,0);
|
||||
Blend(61,80,226,246,0,226,0,0);
|
||||
Blend(81,100,226,0,0,123,0,0);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| WaveletCPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::WaveletCPU(const double &data[],const int datacount,const int x_size,const int y_size,const int i,const int j,const bool norm,double &result[])
|
||||
{
|
||||
double a1=(double)10e-10;
|
||||
double a2=(double)15.0;
|
||||
double da=(double)(a2-a1)/y_size;
|
||||
double db=(double)(datacount-0)/x_size;
|
||||
int pos=j*x_size+i;
|
||||
//---
|
||||
double a=a1+j*da;
|
||||
double b=i*db;
|
||||
double B=(double)1.0; //Morlet
|
||||
double B_inv=(double)1.0/B;
|
||||
double a_inv=(double)1/a;
|
||||
double dt=(double)1.0;
|
||||
double coef=(double)0.0;
|
||||
if(!norm)
|
||||
coef=(double)MathSqrt(a_inv);
|
||||
else
|
||||
{
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
double arg=(dt*k-b)*a_inv;
|
||||
arg=-B_inv*arg*arg;
|
||||
coef+=(double)MathExp(arg);
|
||||
}
|
||||
}
|
||||
double sum=0.0;
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
double arg=(dt*k-b)*a_inv;
|
||||
sum+=data[k]*Morlet(arg);
|
||||
}
|
||||
sum/=coef;
|
||||
result[pos]=sum;
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| CalculateWavelet_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::CalculateWavelet_CPU(const double &data[],uint &time)
|
||||
{
|
||||
time=GetTickCount();
|
||||
int datacount=ArraySize(data);
|
||||
ArrayCopy(m_data,data,0,0,WHOLE_ARRAY);
|
||||
for(int i=0; i<m_xsize; i++)
|
||||
{
|
||||
for(int j=0; j<m_ysize; j++)
|
||||
{
|
||||
WaveletCPU(m_data,datacount,m_xsize,m_ysize,i,j,true,m_wavelet_data_CPU);
|
||||
}
|
||||
}
|
||||
time=GetTickCount()-time;
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| CalculateWavelet_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::CalculateWavelet_GPU(double &data[],uint &time)
|
||||
{
|
||||
int datacount=ArraySize(data);
|
||||
if(!m_OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- check support working with double
|
||||
if(!m_OpenCL.SupportDouble())
|
||||
{
|
||||
PrintFormat("Working with double (cl_khr_fp64) is not supported on the device.");
|
||||
return(false);
|
||||
}
|
||||
//---
|
||||
m_OpenCL.SetKernelsCount(1);
|
||||
m_OpenCL.KernelCreate(0,"Wavelet_GPU");
|
||||
//---
|
||||
m_OpenCL.SetBuffersCount(2);
|
||||
if(!m_OpenCL.BufferFromArray(0,data,0,datacount,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!m_OpenCL.BufferCreate(1,m_xsize*m_ysize*sizeof(double),CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
m_OpenCL.SetArgumentBuffer(0,0,0);
|
||||
m_OpenCL.SetArgumentBuffer(0,4,1);
|
||||
//---
|
||||
ArrayResize(m_wavelet_data_GPU,m_xsize*m_ysize);
|
||||
uint work[2];
|
||||
uint offset[2]={0,0};
|
||||
//--- set dimensions
|
||||
work[0]=m_xsize;
|
||||
work[1]=m_ysize;
|
||||
//--- set parameters and write data to buffer
|
||||
m_OpenCL.SetArgument(0,1,datacount);
|
||||
m_OpenCL.SetArgument(0,2,m_xsize);
|
||||
m_OpenCL.SetArgument(0,3,m_ysize);
|
||||
time=GetTickCount();
|
||||
//--- GPU calculation start
|
||||
if(!m_OpenCL.Execute(0,2,offset,work))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!m_OpenCL.BufferRead(1,m_wavelet_data_GPU,0,0,m_xsize*m_ysize))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for m_wavelet_data_GPU array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finish
|
||||
time=GetTickCount()-time;
|
||||
//---
|
||||
m_OpenCL.Shutdown();
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| ShowWavelet |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::ShowWavelet(const int mode)
|
||||
{
|
||||
if(mode==CPU_DATA)
|
||||
ShowWaveletData(m_wavelet_data_CPU);
|
||||
else
|
||||
if(mode==GPU_DATA)
|
||||
ShowWaveletData(m_wavelet_data_GPU);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| ShowWaveletData |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::ShowWaveletData(const double &m_wavelet_data[])
|
||||
{
|
||||
//--- calculate min/max and range
|
||||
int count=ArraySize(m_wavelet_data);
|
||||
double min_value=m_wavelet_data[0];
|
||||
double max_value=m_wavelet_data[0];
|
||||
for(int i=1; i<count; i++)
|
||||
{
|
||||
min_value=MathMin(min_value,m_wavelet_data[i]);
|
||||
max_value=MathMax(max_value,m_wavelet_data[i]);
|
||||
}
|
||||
double range=max_value-min_value;
|
||||
if(range>0)
|
||||
{
|
||||
for(int j=0; j<m_ysize; j++)
|
||||
{
|
||||
for(int i=0; i<m_xsize; i++)
|
||||
{
|
||||
int pos=j*m_xsize+i;
|
||||
int colindex=int(m_maxcolor*(m_wavelet_data[pos]-min_value)/range);
|
||||
m_bmp_buffer[pos]=GetPalColor(colindex);
|
||||
}
|
||||
}
|
||||
//--- show image
|
||||
ResourceCreate(m_res_name,m_bmp_buffer,m_xsize,m_ysize,0,0,0,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
ChartRedraw();
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Weirstrass function |
|
||||
//+------------------------------------------------------------------+
|
||||
double Weirstrass(double x,double a,double b)
|
||||
{
|
||||
double sum=0.0;
|
||||
double b0=b;
|
||||
double a0=a;
|
||||
for(int n=0; n<35; n++)
|
||||
{
|
||||
double v=b0*(double)MathCos(a0*M_PI*x);
|
||||
sum=sum+v;
|
||||
a0=a0*a;
|
||||
b0=b0*b;
|
||||
}
|
||||
return(sum);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareModelData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PrepareModelData(double &price_data[],const int datacount)
|
||||
{
|
||||
ArrayResize(price_data,datacount);
|
||||
//--- Weirstrass function
|
||||
double x1=0;
|
||||
double x2=2;
|
||||
double dx=(x2-x1)/datacount;
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
price_data[i]=Weirstrass(x1+dx*i,(double)3,(double)0.62);
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PreparePriceData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PreparePriceData(const string symbol,ENUM_TIMEFRAMES timeframe,double &price_data[],const int datacount)
|
||||
{
|
||||
ArrayResize(price_data,datacount);
|
||||
CopyClose(symbol,timeframe,0,datacount,price_data);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareMomentumData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PrepareMomentumData(double &price_data[],double &momentum_data[],const int momentum_period)
|
||||
{
|
||||
int size=ArraySize(price_data);
|
||||
int datacount=size-momentum_period;
|
||||
//---
|
||||
ArrayResize(momentum_data,datacount);
|
||||
for(int i=0; i<datacount; i+=1)
|
||||
{
|
||||
momentum_data[i]=price_data[i+momentum_period]-price_data[i];
|
||||
}
|
||||
ArrayCopy(price_data,price_data,momentum_period,0,datacount);
|
||||
ArrayResize(price_data,datacount);
|
||||
//--- rescale momentum data
|
||||
double min_value=momentum_data[0];
|
||||
double max_value=momentum_data[0];
|
||||
for(int i=1; i<datacount; i++)
|
||||
{
|
||||
double value=momentum_data[i];
|
||||
if(momentum_data[i]>max_value)
|
||||
max_value=value;
|
||||
|
||||
if(momentum_data[i]<min_value)
|
||||
min_value=value;
|
||||
}
|
||||
double range=max_value-min_value;
|
||||
for(int i=0; i<datacount; i+=1)
|
||||
momentum_data[i]=-1+2*(momentum_data[i]-min_value)/range;
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//---
|
||||
int momentum_period=8;
|
||||
double price_data[];
|
||||
double momentum_data[];
|
||||
PrepareModelData(price_data,SIZE_X+momentum_period);
|
||||
//PreparePriceData("EURUSD",PERIOD_M1,price_data,SIZE_X+momentum_period);
|
||||
PrepareMomentumData(price_data,momentum_data,momentum_period);
|
||||
//---
|
||||
CGraphic graph_price;
|
||||
CGraphic graph_momentum;
|
||||
graph_price.Create(0,"price",0,0,0,SIZE_X+130+6,SIZE_Y);
|
||||
graph_price.XAxis().MaxGrace(0);
|
||||
graph_price.HistorySymbolSize(10);
|
||||
graph_price.CurveAdd(price_data,ColorToARGB(clrRed,255),CURVE_LINES,"Price");
|
||||
graph_price.CurvePlotAll();
|
||||
graph_price.Redraw(true);
|
||||
graph_price.Update();
|
||||
//---
|
||||
graph_momentum.Create(0,"momentum",0,0,SIZE_Y,SIZE_X+130+6,SIZE_Y+SIZE_Y);
|
||||
graph_momentum.XAxis().MaxGrace(0);
|
||||
graph_momentum.HistorySymbolSize(10);
|
||||
graph_momentum.CurveAdd(momentum_data,ColorToARGB(clrBlue,255),CURVE_LINES,"Momentum");
|
||||
graph_momentum.CurvePlotAll();
|
||||
graph_momentum.Redraw(true);
|
||||
graph_momentum.Update();
|
||||
//---
|
||||
CWavelet wavelet;
|
||||
//---
|
||||
uint time_cpu=0;
|
||||
wavelet.Create("Wavelet",50,2*SIZE_Y,SIZE_X,SIZE_Y);
|
||||
if(!wavelet.CalculateWavelet_CPU(momentum_data,time_cpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on CPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//wavelet.ShowWavelet(CPU_DATA);
|
||||
uint time_gpu=0;
|
||||
if(!wavelet.CalculateWavelet_GPU(momentum_data,time_gpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on GPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
wavelet.ShowWavelet(GPU_DATA);
|
||||
//---
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
//---
|
||||
PrintFormat("time CPU=%d ms, time GPU=%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- Sleep 10 seconds
|
||||
Sleep(10000);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,213 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort.mq5 |
|
||||
//| Copyright 2016-2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2016-2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
//--- COpenCL class
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
#resource "Kernels/bitonicsort.cl" as string cl_program
|
||||
//+------------------------------------------------------------------+
|
||||
//| QuickSortAscending |
|
||||
//+------------------------------------------------------------------+
|
||||
//| The function sorts array[] QuickSort algorithm. |
|
||||
//| |
|
||||
//| Arguments: |
|
||||
//| array : Array with values to sort |
|
||||
//| first : First element index |
|
||||
//| last : Last element index |
|
||||
//| |
|
||||
//| Return value: None |
|
||||
//+------------------------------------------------------------------+
|
||||
void QuickSortAscending(float &array[],int first,int last)
|
||||
{
|
||||
int i,j;
|
||||
float p_float,t_float;
|
||||
if(first<0 || last<0)
|
||||
return;
|
||||
i=first;
|
||||
j=last;
|
||||
while(i<last)
|
||||
{
|
||||
p_float=array[(first+last)>>1];
|
||||
while(i<j)
|
||||
{
|
||||
while(array[i]<p_float)
|
||||
{
|
||||
if(i==ArraySize(array)-1)
|
||||
break;
|
||||
i++;
|
||||
}
|
||||
while(array[j]>p_float)
|
||||
{
|
||||
if(j==0)
|
||||
break;
|
||||
j--;
|
||||
}
|
||||
if(i<=j)
|
||||
{
|
||||
//-- swap elements i and j
|
||||
t_float=array[i];
|
||||
array[i]=array[j];
|
||||
array[j]=t_float;
|
||||
i++;
|
||||
if(j==0)
|
||||
break;
|
||||
j--;
|
||||
}
|
||||
}
|
||||
if(first<j)
|
||||
QuickSortAscending(array,first,j);
|
||||
first=i;
|
||||
j=last;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| QuickSort_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool QuickSort_CPU(float &data_array[],ulong &time_cpu)
|
||||
{
|
||||
int data_count=ArraySize(data_array);
|
||||
if(data_count<=1)
|
||||
return(false);
|
||||
//--- sort values on CPU
|
||||
time_cpu=GetMicrosecondCount();
|
||||
QuickSortAscending(data_array,0,data_count-1);
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu)/1000);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool BitonicSort_GPU(COpenCL &OpenCL,float &data_array[],ulong &time_gpu)
|
||||
{
|
||||
int data_count=ArraySize(data_array);
|
||||
if(data_count<=1)
|
||||
return(false);
|
||||
|
||||
OpenCL.SetKernelsCount(1);
|
||||
OpenCL.KernelCreate(0,"BitonicSort_GPU");
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(1);
|
||||
if(!OpenCL.BufferFromArray(0,data_array,0,data_count,CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
OpenCL.SetArgumentBuffer(0,0,0);
|
||||
//---
|
||||
uint work_offset[1]={0};
|
||||
uint global_size[1];
|
||||
global_size[0]=data_count>>1;
|
||||
//---
|
||||
uint passes_total=0;
|
||||
uint stages_total=0;
|
||||
//---
|
||||
for(uint temp=data_count; temp>1; temp>>=1)
|
||||
stages_total++;
|
||||
//--- GPU calculation start
|
||||
time_gpu=GetMicrosecondCount();
|
||||
for(uint stage=0; stage<stages_total; stage++)
|
||||
{
|
||||
//--- set stage of the algorithm
|
||||
OpenCL.SetArgument(0,1,stage);
|
||||
for(uint pass=0; pass<stage+1; pass++)
|
||||
{
|
||||
//--- set pass of the current stage
|
||||
OpenCL.SetArgument(0,2,pass);
|
||||
//--- execute kernel
|
||||
if(!OpenCL.Execute(0,1,work_offset,global_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
else
|
||||
passes_total++;
|
||||
}
|
||||
}
|
||||
//---
|
||||
if(!OpenCL.BufferRead(0,data_array,0,0,data_count))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for data array C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finish
|
||||
time_gpu=ulong((GetMicrosecondCount()-time_gpu)/1000);
|
||||
PrintFormat("Bitonic sort finished. Total stages=%d, total passes=%d",stages_total,passes_total);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareDataArray |
|
||||
//+------------------------------------------------------------------+
|
||||
bool PrepareDataArray(ulong global_memory_size,float &data[],int &data_count)
|
||||
{
|
||||
int pwr_max=(int)(MathLog(global_memory_size/sizeof(float))/MathLog(2));
|
||||
int pwr=(int)MathMax(15,pwr_max-4);
|
||||
data_count=(int)MathPow(2,pwr);
|
||||
//--- prepare array and generate random data
|
||||
if(ArrayResize(data,data_count)<data_count)
|
||||
{
|
||||
Print("Error in ArrayResize. Error code=",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
for(int i=0; i<data_count; i++)
|
||||
data[i]=(float)(100000000*MathRand()/32767.0);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//--- OpenCL
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
long global_memory_size=0;
|
||||
if(!OpenCL.GetGlobalMemorySize(global_memory_size))
|
||||
{
|
||||
Print("Error in request of global memory size. Error code=",GetLastError());
|
||||
return;
|
||||
}
|
||||
float data_cpu[];
|
||||
int data_count;
|
||||
//--- prepare array with random values
|
||||
if(PrepareDataArray(global_memory_size,data_cpu,data_count)==false)
|
||||
return;
|
||||
//--- copy array values for sorting on GPU
|
||||
float data_gpu[];
|
||||
if(ArrayCopy(data_gpu,data_cpu,0,0,data_count)!=data_count)
|
||||
return;
|
||||
//--- Quick sort values using CPU
|
||||
ulong time_cpu=0;
|
||||
if(!QuickSort_CPU(data_cpu,time_cpu))
|
||||
return;
|
||||
//--- Bitonic sort values using GPU
|
||||
ulong time_gpu=0;
|
||||
if(!BitonicSort_GPU(OpenCL,data_gpu,time_gpu))
|
||||
return;
|
||||
//--- remove OpenCL objects
|
||||
OpenCL.Shutdown();
|
||||
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
PrintFormat("time CPU=%d ms, time GPU =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- check calculations
|
||||
float total_error=0;
|
||||
for(int i=0; i<data_count; i++)
|
||||
{
|
||||
total_error+=MathAbs(data_gpu[i]-data_cpu[i]);
|
||||
}
|
||||
PrintFormat("Total error = %f",total_error);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,306 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT.mq5 |
|
||||
//| Copyright 2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
|
||||
#include <Math/Stat/Math.mqh>
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
|
||||
#resource "Kernels/fft.cl" as string cl_program
|
||||
#define kernel_init "fft_init"
|
||||
#define kernel_stage "fft_stage"
|
||||
#define kernel_scale "fft_scale"
|
||||
#define NUM_POINTS 16384
|
||||
#define FFT_DIRECTION 1
|
||||
//+------------------------------------------------------------------+
|
||||
//| Fast Fourier transform and its inverse (both recursively) |
|
||||
//| Copyright (C) 2004, Jerome R. Breitenbach. All rights reserved. |
|
||||
//| Reference: |
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Recursive direct FFT transform |
|
||||
//+------------------------------------------------------------------+
|
||||
void fft(const int N,float &x_real[],float &x_imag[],float &X_real[],float &X_imag[])
|
||||
{
|
||||
//--- prepare temporary arrays
|
||||
float XX_real[],XX_imag[];
|
||||
ArrayResize(XX_real,N);
|
||||
ArrayResize(XX_imag,N);
|
||||
//--- calculate FFT by a recursion
|
||||
fft_rec(N,0,1,x_real,x_imag,X_real,X_imag,XX_real,XX_imag);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Recursive inverse FFT transform |
|
||||
//+------------------------------------------------------------------+
|
||||
void ifft(const int N,float &x_real[],float &x_imag[],float &X_real[],float &X_imag[])
|
||||
{
|
||||
int N2=N/2; // half the number of points in IFFT
|
||||
//--- calculate IFFT via reciprocity property of DFT
|
||||
fft(N,X_real,X_imag,x_real,x_imag);
|
||||
x_real[0]=x_real[0]/N;
|
||||
x_imag[0]=x_imag[0]/N;
|
||||
x_real[N2]=x_real[N2]/N;
|
||||
x_imag[N2]=x_imag[N2]/N;
|
||||
for(int i=1; i<N2; i++)
|
||||
{
|
||||
float tmp0=x_real[i]/N;
|
||||
float tmp1=x_imag[i]/N;
|
||||
x_real[i]=x_real[N-i]/N;
|
||||
x_imag[i]=x_imag[N-i]/N;
|
||||
x_real[N-i]=tmp0;
|
||||
x_imag[N-i]=tmp1;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT recursion |
|
||||
//+------------------------------------------------------------------+
|
||||
void fft_rec(const int N,const int offset,const int delta,float &x_real[],float &x_imag[],float &X_real[],float &X_imag[],float &XX_real[],float &XX_imag[])
|
||||
{
|
||||
static const float TWO_PI=(float)(2*M_PI);
|
||||
int N2=N/2; // half the number of points in FFT
|
||||
int k00,k01,k10,k11; // indices for butterflies
|
||||
if(N!=2)
|
||||
{
|
||||
//--- perform recursive step
|
||||
//--- calculate two (N/2)-point DFT's
|
||||
fft_rec(N2,offset,2*delta,x_real,x_imag,XX_real,XX_imag,X_real,X_imag);
|
||||
fft_rec(N2,offset+delta,2*delta,x_real,x_imag,XX_real,XX_imag,X_real,X_imag);
|
||||
//--- combine the two (N/2)-point DFT's into one N-point DFT
|
||||
for(int k=0; k<N2; k++)
|
||||
{
|
||||
k00 = offset + k*delta;
|
||||
k01 = k00 + N2*delta;
|
||||
k10 = offset + 2*k*delta;
|
||||
k11 = k10 + delta;
|
||||
float cs=(float)MathCos(TWO_PI*k/(float)N);
|
||||
float sn=(float)MathSin(TWO_PI*k/(float)N);
|
||||
float tmp0 = cs*XX_real[k11] + sn*XX_imag[k11];
|
||||
float tmp1 = cs*XX_imag[k11] - sn*XX_real[k11];
|
||||
X_real[k01] = XX_real[k10] - tmp0;
|
||||
X_imag[k01] = XX_imag[k10] - tmp1;
|
||||
X_real[k00] = XX_real[k10] + tmp0;
|
||||
X_imag[k00] = XX_imag[k10] + tmp1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
//--- perform 2-point DFT
|
||||
k00=offset;
|
||||
k01=k00+delta;
|
||||
X_real[k01] = x_real[k00] - x_real[k01];
|
||||
X_imag[k01] = x_imag[k00] - x_imag[k01];
|
||||
X_real[k00] = x_real[k00] + x_real[k01];
|
||||
X_imag[k00] = x_imag[k00] + x_imag[k01];
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool FFT_CPU(int direction,int power,float &data_real[],float &data_imag[],ulong &time_cpu)
|
||||
{
|
||||
//--- calculate the number of points
|
||||
int N=1;
|
||||
for(int i=0;i<power;i++)
|
||||
N*=2;
|
||||
//---prepare temporary arrays
|
||||
float XX_real[],XX_imag[];
|
||||
ArrayResize(XX_real,N);
|
||||
ArrayResize(XX_imag,N);
|
||||
//--- CPU calculation start
|
||||
time_cpu=GetMicrosecondCount();
|
||||
if(direction>0)
|
||||
fft(N,data_real,data_imag,XX_real,XX_imag);
|
||||
else
|
||||
ifft(N,XX_real,XX_imag,data_real,data_imag);
|
||||
//--- CPU calculation finished
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu));
|
||||
//--- copy calculated data
|
||||
ArrayCopy(data_real,XX_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(data_imag,XX_imag,0,0,WHOLE_ARRAY);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool FFT_GPU(int direction,int power,float &data_real[],float &data_imag[],ulong &time_gpu)
|
||||
{
|
||||
//--- calculate the number of points
|
||||
int num_points=1;
|
||||
for(int i=0;i<power;i++)
|
||||
num_points*=2;
|
||||
//--- prepare data array for GPU calculation
|
||||
float data[];
|
||||
ArrayResize(data,2*num_points);
|
||||
for(int i=0; i<num_points; i++)
|
||||
{
|
||||
data[2*i]=data_real[i];
|
||||
data[2*i+1]=data_imag[i];
|
||||
}
|
||||
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- create kernels
|
||||
OpenCL.SetKernelsCount(3);
|
||||
OpenCL.KernelCreate(0,kernel_init);
|
||||
OpenCL.KernelCreate(1,kernel_stage);
|
||||
OpenCL.KernelCreate(2,kernel_scale);
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(2);
|
||||
if(!OpenCL.BufferFromArray(0,data,0,2*num_points,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for input buffer. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferCreate(1,2*num_points*sizeof(float),CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for data buffer. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- determine maximum work-group size
|
||||
int local_size=(int)CLGetInfoInteger(OpenCL.GetKernel(0),CL_KERNEL_WORK_GROUP_SIZE);
|
||||
//--- determine local memory size
|
||||
uint local_mem_size=(uint)CLGetInfoInteger(OpenCL.GetContext(),CL_DEVICE_LOCAL_MEM_SIZE);
|
||||
int points_per_group=(int)local_mem_size/(2*sizeof(float));
|
||||
if(points_per_group>num_points)
|
||||
points_per_group=num_points;
|
||||
//--- set kernel arguments
|
||||
OpenCL.SetArgumentBuffer(0,0,0);
|
||||
OpenCL.SetArgumentBuffer(0,1,1);
|
||||
OpenCL.SetArgumentLocalMemory(0,2,local_mem_size);
|
||||
OpenCL.SetArgument(0,3,points_per_group);
|
||||
OpenCL.SetArgument(0,4,num_points);
|
||||
OpenCL.SetArgument(0,5,direction);
|
||||
//--- OpenCL execute settings
|
||||
int task_dimension=1;
|
||||
uint global_size=(uint)((num_points/points_per_group)*local_size);
|
||||
uint global_work_offset[1]={0};
|
||||
uint global_work_size[1];
|
||||
global_work_size[0]=global_size;
|
||||
uint local_work_size[1];
|
||||
local_work_size[0]=local_size;
|
||||
//--- GPU calculation start
|
||||
time_gpu=GetMicrosecondCount();
|
||||
//-- execute kernel fft_init
|
||||
if(!OpenCL.Execute(0,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_init: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//-- further stages of the FFT
|
||||
if(num_points>points_per_group)
|
||||
{
|
||||
//--- set arguments for kernel 1
|
||||
OpenCL.SetArgumentBuffer(1,0,1);
|
||||
OpenCL.SetArgument(1,2,points_per_group);
|
||||
OpenCL.SetArgument(1,3,direction);
|
||||
for(int stage=2; stage<=num_points/points_per_group; stage<<=1)
|
||||
{
|
||||
OpenCL.SetArgument(1,1,stage);
|
||||
//-- execute kernel fft_stage
|
||||
if(!OpenCL.Execute(1,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_stage: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
}
|
||||
}
|
||||
//--- scale values if performing the inverse FFT
|
||||
if(direction<0)
|
||||
{
|
||||
OpenCL.SetArgumentBuffer(2,0,1);
|
||||
OpenCL.SetArgument(2,1,points_per_group);
|
||||
OpenCL.SetArgument(2,2,num_points);
|
||||
//-- execute kernel fft_scale
|
||||
if(!OpenCL.Execute(2,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("fft_scale: Error in CLExecute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
}
|
||||
//--- read the results from GPU memory
|
||||
if(!OpenCL.BufferRead(1,data,0,0,2*num_points))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for data_buffer2. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time_gpu=ulong((GetMicrosecondCount()-time_gpu));
|
||||
//--- copy calculated data and release OpenCL handles
|
||||
for(int i=0; i<num_points; i++)
|
||||
{
|
||||
data_real[i]=data[2*i];
|
||||
data_imag[i]=data[2*i+1];
|
||||
}
|
||||
OpenCL.Shutdown();
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
int datacount=NUM_POINTS;
|
||||
int power=(int)(MathLog(NUM_POINTS)/M_LN2);
|
||||
if(MathPow(2,power)!=datacount)
|
||||
{
|
||||
PrintFormat("Number of elements must be power of 2. Elements: %d",datacount);
|
||||
return;
|
||||
}
|
||||
//--- prepare data for FFT calculation
|
||||
float data_real[],data_imag[];
|
||||
ArrayResize(data_real,datacount);
|
||||
ArrayResize(data_imag,datacount);
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
data_real[i]=(float)i;
|
||||
data_imag[i]=0;
|
||||
}
|
||||
|
||||
int direction=FFT_DIRECTION;
|
||||
//--- data arrays for CPU calculation
|
||||
float CPU_real[],CPU_imag[];
|
||||
ArrayCopy(CPU_real,data_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(CPU_imag,data_imag,0,0,WHOLE_ARRAY);
|
||||
ulong time_cpu=0;
|
||||
//--- calculate FFT using CPU
|
||||
FFT_CPU(direction,power,CPU_real,CPU_imag,time_cpu);
|
||||
|
||||
//--- data arrays for GPU calculation
|
||||
float GPU_real[],GPU_imag[];
|
||||
ArrayCopy(GPU_real,data_real,0,0,WHOLE_ARRAY);
|
||||
ArrayCopy(GPU_imag,data_imag,0,0,WHOLE_ARRAY);
|
||||
ulong time_gpu=0;
|
||||
//--- calculate FFT using GPU
|
||||
if(!FFT_GPU(direction,power,GPU_real,GPU_imag,time_gpu))
|
||||
{
|
||||
PrintFormat("Error in calculation FFT on GPU.");
|
||||
return;
|
||||
}
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
PrintFormat("FFT calculation for %d points.",datacount);
|
||||
PrintFormat("time CPU=%d microseconds, time GPU =%d microseconds, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- determine average error
|
||||
float average_error=0.0;
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
average_error += (float)MathAbs(CPU_real[i]-GPU_real[i]);
|
||||
average_error += (float)MathAbs(CPU_imag[i]-GPU_imag[i]);
|
||||
}
|
||||
average_error=average_error/(datacount*2);
|
||||
PrintFormat("Average error = %f",average_error);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,28 @@
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel |
|
||||
//| The bitonic sort kernel does an ascending sort. |
|
||||
//+-----------------------------------------------------------+
|
||||
//| R. Banger,K. Bhattacharyya, OpenCL Programming by Example:|
|
||||
//| A comprehensive guide on OpenCL programming with examples |
|
||||
//| PACKT Publishing, 2013. |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void BitonicSort_GPU(__global float *data,const uint stage,const uint pass)
|
||||
{
|
||||
uint id=get_global_id(0);
|
||||
uint distance = 1<<(stage-pass);
|
||||
uint left_id =(id &(distance-1));
|
||||
left_id+=(id>>(stage-pass))*(distance<<1);
|
||||
uint right_id=left_id+distance;
|
||||
float left_value=data[left_id];
|
||||
float right_value=data[right_id];
|
||||
uint same_direction=(id>>stage)&0x1;
|
||||
uint temp = same_direction?right_id:temp;
|
||||
right_id = same_direction?left_id:right_id;
|
||||
left_id = same_direction?temp:left_id;
|
||||
int compare_res=(left_value<right_value);
|
||||
float greater = compare_res?right_value:left_value;
|
||||
float lesser = compare_res?left_value:right_value;
|
||||
data[left_id] = lesser;
|
||||
data[right_id]= greater;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,153 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_init OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_init(__global float2 *in_data,
|
||||
__global float2 *out_data,
|
||||
__local float2 *l_data,
|
||||
uint points_per_group,uint size,int dir)
|
||||
{
|
||||
uint4 br,index;
|
||||
uint points_per_item,g_addr,l_addr,i,fft_index,stage,N2;
|
||||
float2 x1,x2,x3,x4,sum12,diff12,sum34,diff34;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
l_addr = get_local_id(0)*points_per_item;
|
||||
g_addr = get_group_id(0)*points_per_group + l_addr;
|
||||
//--- load data from bit-reversed addresses and perform 4-point FFTs
|
||||
for(i=0; i<points_per_item; i+=4)
|
||||
{
|
||||
index=(uint4)(g_addr,g_addr+1,g_addr+2,g_addr+3);
|
||||
fft_index=size/2;
|
||||
stage=1;
|
||||
N2 =(uint)log2((float)size)-1;
|
||||
br =(index<< N2) & fft_index;
|
||||
br|=(index>> N2) & stage;
|
||||
//--- bit-reverse addresses
|
||||
while(N2>1)
|
||||
{
|
||||
N2-=2;
|
||||
fft_index>>=1;
|
||||
stage<<=1;
|
||||
br |= (index << N2) & fft_index;
|
||||
br |= (index >> N2) & stage;
|
||||
}
|
||||
//--- load global data
|
||||
x1 = in_data[br.s0];
|
||||
x2 = in_data[br.s1];
|
||||
x3 = in_data[br.s2];
|
||||
x4 = in_data[br.s3];
|
||||
|
||||
sum12=x1+x2;
|
||||
diff12= x1-x2;
|
||||
sum34 = x3+x4;
|
||||
diff34=(float2)(x3.s1-x4.s1,x4.s0-x3.s0)*dir;
|
||||
l_data[l_addr]=sum12+sum34;
|
||||
l_data[l_addr+1] = diff12 + diff34;
|
||||
l_data[l_addr+2] = sum12 - sum34;
|
||||
l_data[l_addr+3] = diff12 - diff34;
|
||||
l_addr += 4;
|
||||
g_addr += 4;
|
||||
}
|
||||
//--- perform initial stages of the FFT - each of length N2*2
|
||||
for(N2=4; N2<points_per_item; N2<<=1)
|
||||
{
|
||||
l_addr=get_local_id(0)*points_per_item;
|
||||
for(fft_index=0; fft_index<points_per_item; fft_index+=2*N2)
|
||||
{
|
||||
x1=l_data[l_addr];
|
||||
l_data[l_addr]+=l_data[l_addr+N2];
|
||||
l_data[l_addr+N2]=x1-l_data[l_addr+N2];
|
||||
for(i=1; i<N2; i++)
|
||||
{
|
||||
x3.s0=cos(M_PI_F*i/N2);
|
||||
x3.s1=dir*sin(M_PI_F*i/N2);
|
||||
x2=(float2)(l_data[l_addr+N2+i].s0*x3.s0+l_data[l_addr+N2+i].s1*x3.s1,l_data[l_addr+N2+i].s1*x3.s0-l_data[l_addr+N2+i].s0*x3.s1);
|
||||
l_data[l_addr+N2+i]=l_data[l_addr+i]-x2;
|
||||
l_data[l_addr+i]+=x2;
|
||||
}
|
||||
l_addr+=2*N2;
|
||||
}
|
||||
}
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
//--- perform FFT with other items in group - each of length N2*2
|
||||
stage=2;
|
||||
for(N2=points_per_item; N2<points_per_group; N2<<=1)
|
||||
{
|
||||
br.s0=(get_local_id(0)+(get_local_id(0)/stage)*stage) *(points_per_item/2);
|
||||
size = br.s0 % (N2*2);
|
||||
for(i=br.s0; i<br.s0+points_per_item/2; i++)
|
||||
{
|
||||
x3.s0=cos(M_PI_F*size/N2);
|
||||
x3.s1=dir*sin(M_PI_F*size/N2);
|
||||
x2=(float2)(l_data[N2+i].s0*x3.s0+l_data[N2+i].s1*x3.s1,l_data[N2+i].s1*x3.s0-l_data[N2+i].s0*x3.s1);
|
||||
l_data[N2+i]=l_data[i]-x2;
|
||||
l_data[i]+=x2;
|
||||
size++;
|
||||
}
|
||||
stage<<=1;
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
}
|
||||
//--- store results in global memory
|
||||
l_addr = get_local_id(0)*points_per_item;
|
||||
g_addr = get_group_id(0)*points_per_group + l_addr;
|
||||
for(i=0; i<points_per_item; i+=4)
|
||||
{
|
||||
out_data[g_addr]=l_data[l_addr];
|
||||
out_data[g_addr+1] = l_data[l_addr+1];
|
||||
out_data[g_addr+2] = l_data[l_addr+2];
|
||||
out_data[g_addr+3] = l_data[l_addr+3];
|
||||
g_addr += 4;
|
||||
l_addr += 4;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_stage OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_stage(__global float2 *g_data,uint stage,uint points_per_group,int dir)
|
||||
{
|
||||
uint points_per_item,addr,N,ang,i;
|
||||
float c,s;
|
||||
float2 input1,input2,w;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
addr=(get_group_id(0)+(get_group_id(0)/stage)*stage)*(points_per_group/2)+get_local_id(0)*(points_per_item/2);
|
||||
N=points_per_group*(stage/2);
|
||||
ang=addr%(N*2);
|
||||
|
||||
for(i=addr; i<addr+points_per_item/2; i++)
|
||||
{
|
||||
c = cos(M_PI_F*ang/N);
|
||||
s = dir*sin(M_PI_F*ang/N);
|
||||
input1 = g_data[i];
|
||||
input2 = g_data[i+N];
|
||||
w=(float2)(input2.s0*c+input2.s1*s,input2.s1*c-input2.s0*s);
|
||||
g_data[i]=input1+w;
|
||||
g_data[i+N]=input1-w;
|
||||
ang++;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| fft_scale OpenCL kernel for Fast Fourier Transfrom |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Matthew Scarpino, "OpenCL in Action: How to accelerate graphics |
|
||||
//| and computations", Manning, 2012, Chapter 14. |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void fft_scale(__global float2 *g_data,uint points_per_group,uint scale)
|
||||
{
|
||||
uint points_per_item,addr,i;
|
||||
|
||||
points_per_item=points_per_group/get_local_size(0);
|
||||
addr=get_group_id(0)*points_per_group+get_local_id(0)*points_per_item;
|
||||
|
||||
for(i=addr; i<addr+points_per_item; i++)
|
||||
{
|
||||
g_data[i]/=scale;
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,66 @@
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel for matrix multiplication |
|
||||
//| using global work groups |
|
||||
//+-----------------------------------------------------------+
|
||||
//| http://gpgpu-computing4.blogspot.ru/2009/09/ |
|
||||
//| /matrix-multiplication-2-opencl.html |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void MatrixMult_GPU1(__global float *matrix_a,
|
||||
__global float *matrix_b,
|
||||
__global float *matrix_c,
|
||||
int rows_a,int cols_a,int cols_b)
|
||||
{
|
||||
int i=get_global_id(0);
|
||||
int j=get_global_id(1);
|
||||
float sum=0.0;
|
||||
for(int k=0; k<cols_a; k++)
|
||||
{
|
||||
sum+=matrix_a[cols_a*i+k]*matrix_b[cols_b*k+j];
|
||||
}
|
||||
matrix_c[cols_b*i+j]=sum;
|
||||
}
|
||||
#define BLOCK_SIZE 10
|
||||
//+-----------------------------------------------------------+
|
||||
//| OpenCL kernel for matrix multiplication |
|
||||
//| using local groups with common local memory |
|
||||
//+-----------------------------------------------------------+
|
||||
//| http://gpgpu-computing4.blogspot.ru/2009/10/ |
|
||||
//| /matrix-multiplication-3-opencl.html |
|
||||
//+-----------------------------------------------------------+
|
||||
__kernel void MatrixMult_GPU2(__global float *matrix_a,
|
||||
__global float *matrix_b,
|
||||
__global float *matrix_c,
|
||||
int rows_a,int cols_a,int cols_b)
|
||||
{
|
||||
int group_i=get_group_id(0);
|
||||
int group_j=get_group_id(1);
|
||||
|
||||
int i=get_local_id(0);
|
||||
int j=get_local_id(1);
|
||||
|
||||
int offset_b=BLOCK_SIZE*group_i;
|
||||
int offset_a_start=cols_a*BLOCK_SIZE*group_j;
|
||||
float sum=(float)0.0;
|
||||
|
||||
for(int offset_a=offset_a_start;
|
||||
offset_a<offset_a_start+cols_a;
|
||||
offset_a+=BLOCK_SIZE,
|
||||
offset_b+=BLOCK_SIZE*cols_b)
|
||||
{
|
||||
__local float submatrix_a[BLOCK_SIZE][BLOCK_SIZE];
|
||||
__local float submatrix_b[BLOCK_SIZE][BLOCK_SIZE];
|
||||
|
||||
submatrix_a[i][j]=matrix_a[offset_a+cols_a*i+j];
|
||||
submatrix_b[i][j]=matrix_b[offset_b+cols_b*i+j];
|
||||
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
|
||||
for(int k=0; k<BLOCK_SIZE; k++)
|
||||
sum+=submatrix_a[i][k]*submatrix_b[k][j];
|
||||
|
||||
barrier(CLK_LOCAL_MEM_FENCE);
|
||||
}
|
||||
int offset_c=BLOCK_SIZE*(cols_b*group_j+group_i);
|
||||
matrix_c[offset_c+cols_b*i+j]=sum;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,48 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Morlet wavelet function |
|
||||
//+------------------------------------------------------------------+
|
||||
float Morlet(const float t)
|
||||
{
|
||||
return exp(-t*t*0.5)*cos(M_2_PI*t);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| OpenCL kernel function |
|
||||
//+------------------------------------------------------------------+
|
||||
__kernel void Wavelet_GPU(__global float *data,int datacount,int x_size,int y_size,__global float *result)
|
||||
{
|
||||
size_t i = get_global_id(0);
|
||||
size_t j = get_global_id(1);
|
||||
float a1=(float)10e-10;
|
||||
float a2=(float)15.0;
|
||||
float da=(a2-a1)/(float)y_size;
|
||||
float db=((float)datacount-(float)0.0)/x_size;
|
||||
float a=a1+j*da;
|
||||
float b=0+i*db;
|
||||
uint norm=1;
|
||||
float B=(float)1.0; //Morlet
|
||||
float B_inv=(float)1.0/B;
|
||||
float a_inv=(float)1.0/a;
|
||||
float dt=(float)1.0;
|
||||
float coef=(float)0.0;
|
||||
if(norm==0)
|
||||
coef=sqrt(a_inv);
|
||||
else
|
||||
{
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
float arg=(dt*k-b)*a_inv;
|
||||
arg=-B_inv*arg*arg;
|
||||
coef=coef+exp(arg);
|
||||
}
|
||||
}
|
||||
float sum=(float)0.0;
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
float arg=(dt*k-b)*a_inv;
|
||||
sum+=data[k]*Morlet(arg);
|
||||
}
|
||||
sum=sum/coef;
|
||||
uint pos=(int)(j*x_size+i);
|
||||
result[pos]=sum;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,219 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult.mq5 |
|
||||
//| Copyright 2016-2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2016-2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
//--- OpenCL kernels
|
||||
#resource "Kernels/matrixmult.cl" as string cl_program
|
||||
#define BLOCK_SIZE 10
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool MatrixMult_CPU(const float &matrix_a[],const float &matrix_b[],float &matrix_c[],
|
||||
const int rows_a,const int cols_a,const int cols_b,ulong &time_cpu)
|
||||
{
|
||||
int size=rows_a*cols_b;
|
||||
if(ArrayResize(matrix_c,size)!=size)
|
||||
return(false);
|
||||
//--- CPU calculation started
|
||||
time_cpu=GetMicrosecondCount();
|
||||
for(int i=0; i<rows_a; i++)
|
||||
{
|
||||
for(int j=0; j<cols_b; j++)
|
||||
{
|
||||
float sum=0.0;
|
||||
for(int k=0; k<cols_a; k++)
|
||||
{
|
||||
sum+=matrix_a[cols_a*i+k]*matrix_b[cols_b*k+j];
|
||||
}
|
||||
matrix_c[cols_b*i+j]=sum;
|
||||
}
|
||||
}
|
||||
//--- CPU calculation finished
|
||||
time_cpu=ulong((GetMicrosecondCount()-time_cpu)/1000);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool MatrixMult_GPU(const float &matrix_a[],const float &matrix_b[],float &matrix1_c[],float &matrix2_c[],
|
||||
const int rows_a,const int cols_a,const int cols_b,const int size_a,const int size_b,
|
||||
const int size_c,ulong &time1_gpu,ulong &time2_gpu)
|
||||
{
|
||||
const int task_dimension=2;
|
||||
//--- prepare matrices for result
|
||||
if(ArrayResize(matrix1_c,size_c)!=size_c || ArrayResize(matrix2_c,size_c)!=size_c)
|
||||
return(false);
|
||||
ArrayFill(matrix1_c,0,size_c,(float)0.0);
|
||||
ArrayFill(matrix2_c,0,size_c,(float)0.0);
|
||||
//--- OpenCL
|
||||
COpenCL OpenCL;
|
||||
if(!OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- create kernels
|
||||
OpenCL.SetKernelsCount(2);
|
||||
OpenCL.KernelCreate(0,"MatrixMult_GPU1");
|
||||
OpenCL.KernelCreate(1,"MatrixMult_GPU2");
|
||||
//--- create buffers
|
||||
OpenCL.SetBuffersCount(3);
|
||||
//---
|
||||
if(!OpenCL.BufferFromArray(0,matrix_a,0,size_a,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for matrix A. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferFromArray(1,matrix_b,0,size_b,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for matrix B. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferCreate(2,size_c*sizeof(float),CL_MEM_WRITE_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for matrix C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- prepare arguments for kernel 0
|
||||
int kernel_index=0;
|
||||
OpenCL.SetArgumentBuffer(kernel_index,0,0);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,1,1);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,2,2);
|
||||
OpenCL.SetArgument(kernel_index,3,rows_a);
|
||||
OpenCL.SetArgument(kernel_index,4,cols_a);
|
||||
OpenCL.SetArgument(kernel_index,5,cols_b);
|
||||
//--- set task dimension a_rows x b_cols
|
||||
uint global_work_size[2];
|
||||
//--- set dimensions
|
||||
global_work_size[0]=rows_a;
|
||||
global_work_size[1]=cols_b;
|
||||
uint global_work_offset[2]={0,0};
|
||||
//--- GPU calculation start kernel 0
|
||||
time1_gpu=GetMicrosecondCount();
|
||||
if(!OpenCL.Execute(kernel_index,task_dimension,global_work_offset,global_work_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferRead(2,matrix1_c,0,0,size_c))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for matrix1 C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time1_gpu=ulong((GetMicrosecondCount()-time1_gpu)/1000);
|
||||
|
||||
//--- prepare arguments for kernel 1
|
||||
kernel_index=1;
|
||||
//--- set arguments
|
||||
OpenCL.SetArgumentBuffer(kernel_index,0,0);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,1,1);
|
||||
OpenCL.SetArgumentBuffer(kernel_index,2,2);
|
||||
OpenCL.SetArgument(kernel_index,3,rows_a);
|
||||
OpenCL.SetArgument(kernel_index,4,cols_a);
|
||||
OpenCL.SetArgument(kernel_index,5,cols_b);
|
||||
uint local_work_size[2];
|
||||
local_work_size[0]=BLOCK_SIZE;
|
||||
local_work_size[1]=BLOCK_SIZE;
|
||||
//--- GPU calculation start, kernel1
|
||||
time2_gpu=GetMicrosecondCount();
|
||||
if(!OpenCL.Execute(kernel_index,task_dimension,global_work_offset,global_work_size,local_work_size))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!OpenCL.BufferRead(2,matrix2_c,0,0,size_c))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for matrix2 C. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finished
|
||||
time2_gpu=ulong((GetMicrosecondCount()-time2_gpu)/1000);
|
||||
//--- remove OpenCL objects
|
||||
OpenCL.Shutdown();
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//--- matrix A 1000x2000
|
||||
int rows_a=1000;
|
||||
int cols_a=2000;
|
||||
//--- matrix B 2000x1000
|
||||
int rows_b=cols_a;
|
||||
int cols_b=1000;
|
||||
//--- matrix C 1000x1000
|
||||
int rows_c=rows_a;
|
||||
int cols_c=cols_b;
|
||||
//--- matrix A: size=rows_a*cols_a
|
||||
int size_a=rows_a*cols_a;
|
||||
int size_b=rows_b*cols_b;
|
||||
int size_c=rows_c*cols_c;
|
||||
//--- prepare matrix A
|
||||
float matrix_a[];
|
||||
ArrayResize(matrix_a,rows_a*cols_a);
|
||||
for(int i=0; i<rows_a; i++)
|
||||
for(int j=0; j<cols_a; j++)
|
||||
{
|
||||
matrix_a[i*cols_a+j]=(float)(10*MathRand()/32767);
|
||||
}
|
||||
//--- prepare matrix B
|
||||
float matrix_b[];
|
||||
ArrayResize(matrix_b,rows_b*cols_b);
|
||||
for(int i=0; i<rows_b; i++)
|
||||
for(int j=0; j<cols_b; j++)
|
||||
{
|
||||
matrix_b[i*cols_b+j]=(float)(10*MathRand()/32767);
|
||||
}
|
||||
//--- CPU: calculate matrix product matrix_a*matrix_b
|
||||
float matrix_c_cpu[];
|
||||
ulong time_cpu=0;
|
||||
if(!MatrixMult_CPU(matrix_a,matrix_b,matrix_c_cpu,rows_a,cols_a,cols_b,time_cpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on CPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//--- calculate matrix product using GPU
|
||||
float matrix_c_gpu_method1[];
|
||||
float matrix_c_gpu_method2[];
|
||||
ulong time_gpu_method1=0;
|
||||
ulong time_gpu_method2=0;
|
||||
if(!MatrixMult_GPU(matrix_a,matrix_b,matrix_c_gpu_method1,matrix_c_gpu_method2,rows_a,cols_a,cols_b,size_a,size_b,size_c,time_gpu_method1,time_gpu_method2))
|
||||
{
|
||||
PrintFormat("Error in calculation on GPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//--- calculate CPU/GPU ratio
|
||||
double CPU_GPU_ratio1=0;
|
||||
double CPU_GPU_ratio2=0;
|
||||
if(time_gpu_method1!=0)
|
||||
CPU_GPU_ratio1=1.0*time_cpu/time_gpu_method1;
|
||||
if(time_gpu_method2!=0)
|
||||
CPU_GPU_ratio2=1.0*time_cpu/time_gpu_method2;
|
||||
PrintFormat("time CPU=%d ms, time GPU global work groups =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu_method1,CPU_GPU_ratio1);
|
||||
PrintFormat("time CPU=%d ms, time GPU local work groups =%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu_method2,CPU_GPU_ratio2);
|
||||
//--- check calculations
|
||||
float total_error1=0;
|
||||
float total_error2=0;
|
||||
for(int i=0; i<rows_c; i++)
|
||||
{
|
||||
for(int j=0; j<cols_c; j++)
|
||||
{
|
||||
int pos=cols_c*i+j;
|
||||
total_error1+=MathAbs(matrix_c_gpu_method1[pos]-matrix_c_cpu[pos]);
|
||||
total_error2+=MathAbs(matrix_c_gpu_method2[pos]-matrix_c_cpu[pos]);
|
||||
}
|
||||
}
|
||||
PrintFormat("Total error for method 1 = %f",total_error1);
|
||||
PrintFormat("Total error for method 2 = %f",total_error2);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,436 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Wavelet.mq5 |
|
||||
//| Copyright 2016-2017, MetaQuotes Software Corp. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2016-2017, MetaQuotes Software Corp."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
|
||||
#include <Math/Stat/Math.mqh>
|
||||
#include <Graphics/Graphic.mqh>
|
||||
#include <OpenCL/OpenCL.mqh>
|
||||
|
||||
#define CPU_DATA 1
|
||||
#define GPU_DATA 2
|
||||
|
||||
#define SIZE_X 600
|
||||
#define SIZE_Y 200
|
||||
|
||||
#resource "Kernels/wavelet.cl" as string cl_program
|
||||
//+------------------------------------------------------------------+
|
||||
//| CWavelet |
|
||||
//+------------------------------------------------------------------+
|
||||
class CWavelet
|
||||
{
|
||||
protected:
|
||||
int m_xsize;
|
||||
int m_ysize;
|
||||
int m_maxcolor;
|
||||
string m_res_name;
|
||||
string m_label_name;
|
||||
uchar m_palette[3*256];
|
||||
//---
|
||||
float m_data[];
|
||||
float m_wavelet_data_CPU[];
|
||||
float m_wavelet_data_GPU[];
|
||||
uint m_bmp_buffer[];
|
||||
|
||||
COpenCL m_OpenCL;
|
||||
|
||||
float Morlet(const float t);
|
||||
void ShowWaveletData(const float &m_wavelet_data[]);
|
||||
int GetPalColor(const int index);
|
||||
void Blend(const uint c1,const uint c2,const uint r1,const uint g1,const uint b1,const uint r2,const uint g2,const uint b2);
|
||||
bool WaveletCPU(const float &data[],const int datacount,const int x_size,const int y_size,const int i,const int j,const bool norm,float &result[]);
|
||||
public:
|
||||
//---
|
||||
void Create(const string name,const int x0,const int y0,const int x_size,const int y_size);
|
||||
bool CalculateWavelet_CPU(const float &data[],uint &time);
|
||||
bool CalculateWavelet_GPU(float &data[],uint &time);
|
||||
void ShowWavelet(const int mode);
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
//| Morlet wavelet function |
|
||||
//+------------------------------------------------------------------+
|
||||
float CWavelet::Morlet(const float t)
|
||||
{
|
||||
double v=t;
|
||||
double res=MathExp(-v*v*0.5)*MathCos(M_2_PI*v);
|
||||
return ((float)res);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| GetPalColor |
|
||||
//+------------------------------------------------------------------+
|
||||
int CWavelet::GetPalColor(const int index)
|
||||
{
|
||||
int ind=index;
|
||||
if(ind<=0)
|
||||
ind=0;
|
||||
if(ind>255)
|
||||
ind=255;
|
||||
int idx=3*(ind);
|
||||
uchar r=m_palette[idx];
|
||||
uchar g=m_palette[idx+1];
|
||||
uchar b=m_palette[idx+2];
|
||||
//---
|
||||
return(b+256*g+65536*r);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Gradient palette |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::Blend(const uint c1,const uint c2,const uint r1,const uint g1,const uint b1,const uint r2,const uint g2,const uint b2)
|
||||
{
|
||||
int n=int(c2-c1);
|
||||
for(int i=0; i<=n; i++)
|
||||
{
|
||||
if((c1+i+2)<ArraySize(m_palette))
|
||||
{
|
||||
m_palette[3*(c1+i)]=uchar(MathRound(1*(r1*(n-i)+r2*i)*1.0/n));
|
||||
m_palette[3*(c1+i)+1]=uchar(MathRound(1*(g1*(n-i)+g2*i)*1.0/n));
|
||||
m_palette[3*(c1+i)+2]=uchar(MathRound(1*(b1*(n-i)+b2*i)*1.0/n));
|
||||
}
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Create |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::Create(const string name,const int x0,const int y0,const int x_size,const int y_size)
|
||||
{
|
||||
//---
|
||||
m_xsize=x_size;
|
||||
m_ysize=y_size;
|
||||
int size=m_xsize*m_ysize;
|
||||
ArrayResize(m_bmp_buffer,size);
|
||||
ArrayFill(m_bmp_buffer,0,size,0);
|
||||
ArrayResize(m_wavelet_data_CPU,size);
|
||||
ArrayResize(m_wavelet_data_GPU,size);
|
||||
ArrayFill(m_wavelet_data_CPU,0,size,0);
|
||||
ArrayFill(m_wavelet_data_GPU,0,size,0);
|
||||
m_res_name=name;
|
||||
m_label_name=m_res_name;
|
||||
StringToUpper(m_label_name);
|
||||
ResourceCreate(m_res_name,m_bmp_buffer,m_xsize,m_ysize,0,0,0,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
ObjectCreate(0,m_label_name,OBJ_BITMAP_LABEL,0,0,0);
|
||||
ObjectSetInteger(0,m_label_name,OBJPROP_XDISTANCE,x0);
|
||||
ObjectSetInteger(0,m_label_name,OBJPROP_YDISTANCE,y0);
|
||||
ObjectSetString(0,m_label_name,OBJPROP_BMPFILE,NULL);
|
||||
ObjectSetString(0,m_label_name,OBJPROP_BMPFILE,"::"+m_label_name);
|
||||
m_maxcolor=100;
|
||||
//Blend(0,100,0,0,0,255,255,255);
|
||||
Blend(0,20,0,0,95,0,0,246);
|
||||
Blend(21,40,0,0,246,0,236,226);
|
||||
Blend(41,60,0,236,226,226,246,0);
|
||||
Blend(61,80,226,246,0,226,0,0);
|
||||
Blend(81,100,226,0,0,123,0,0);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| WaveletCPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::WaveletCPU(const float &data[],const int datacount,const int x_size,const int y_size,const int i,const int j,const bool norm,float &result[])
|
||||
{
|
||||
float a1=(float)10e-10;
|
||||
float a2=(float)15.0;
|
||||
float da=(float)(a2-a1)/y_size;
|
||||
float db=(float)(datacount-0)/x_size;
|
||||
int pos=j*x_size+i;
|
||||
//---
|
||||
float a=a1+j*da;
|
||||
float b=i*db;
|
||||
float B=(float)1.0; //Morlet
|
||||
float B_inv=(float)1.0/B;
|
||||
float a_inv=(float)1/a;
|
||||
float dt=(float)1.0;
|
||||
float coef=(float)0.0;
|
||||
if(!norm)
|
||||
coef=(float)MathSqrt(a_inv);
|
||||
else
|
||||
{
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
float arg=(dt*k-b)*a_inv;
|
||||
arg=-B_inv*arg*arg;
|
||||
coef+=(float)MathExp(arg);
|
||||
}
|
||||
}
|
||||
float sum=0.0;
|
||||
for(int k=0; k<datacount; k++)
|
||||
{
|
||||
float arg=(dt*k-b)*a_inv;
|
||||
sum+=data[k]*Morlet(arg);
|
||||
}
|
||||
sum/=coef;
|
||||
result[pos]=sum;
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| CalculateWavelet_CPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::CalculateWavelet_CPU(const float &data[],uint &time)
|
||||
{
|
||||
time=GetTickCount();
|
||||
int datacount=ArraySize(data);
|
||||
ArrayCopy(m_data,data,0,0,WHOLE_ARRAY);
|
||||
for(int i=0; i<m_xsize; i++)
|
||||
{
|
||||
for(int j=0; j<m_ysize; j++)
|
||||
{
|
||||
WaveletCPU(m_data,datacount,m_xsize,m_ysize,i,j,true,m_wavelet_data_CPU);
|
||||
}
|
||||
}
|
||||
time=GetTickCount()-time;
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| CalculateWavelet_GPU |
|
||||
//+------------------------------------------------------------------+
|
||||
bool CWavelet::CalculateWavelet_GPU(float &data[],uint &time)
|
||||
{
|
||||
int datacount=ArraySize(data);
|
||||
if(!m_OpenCL.Initialize(cl_program,true))
|
||||
{
|
||||
PrintFormat("Error in OpenCL initialization. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//---
|
||||
m_OpenCL.SetKernelsCount(1);
|
||||
m_OpenCL.KernelCreate(0,"Wavelet_GPU");
|
||||
//---
|
||||
m_OpenCL.SetBuffersCount(2);
|
||||
if(!m_OpenCL.BufferFromArray(0,data,0,datacount,CL_MEM_READ_ONLY))
|
||||
{
|
||||
PrintFormat("Error in BufferFromArray for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!m_OpenCL.BufferCreate(1,m_xsize*m_ysize*sizeof(float),CL_MEM_READ_WRITE))
|
||||
{
|
||||
PrintFormat("Error in BufferCreate for data array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
m_OpenCL.SetArgumentBuffer(0,0,0);
|
||||
m_OpenCL.SetArgumentBuffer(0,4,1);
|
||||
//---
|
||||
ArrayResize(m_wavelet_data_GPU,m_xsize*m_ysize);
|
||||
uint work[2];
|
||||
uint offset[2]={0,0};
|
||||
//--- set dimensions
|
||||
work[0]=m_xsize;
|
||||
work[1]=m_ysize;
|
||||
//--- set parameters and write data to buffer
|
||||
m_OpenCL.SetArgument(0,1,datacount);
|
||||
m_OpenCL.SetArgument(0,2,m_xsize);
|
||||
m_OpenCL.SetArgument(0,3,m_ysize);
|
||||
time=GetTickCount();
|
||||
//--- GPU calculation start
|
||||
if(!m_OpenCL.Execute(0,2,offset,work))
|
||||
{
|
||||
PrintFormat("Error in Execute. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
if(!m_OpenCL.BufferRead(1,m_wavelet_data_GPU,0,0,m_xsize*m_ysize))
|
||||
{
|
||||
PrintFormat("Error in BufferRead for m_wavelet_data_GPU array. Error code=%d",GetLastError());
|
||||
return(false);
|
||||
}
|
||||
//--- GPU calculation finish
|
||||
time=GetTickCount()-time;
|
||||
//---
|
||||
m_OpenCL.Shutdown();
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| ShowWavelet |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::ShowWavelet(const int mode)
|
||||
{
|
||||
if(mode==CPU_DATA)
|
||||
ShowWaveletData(m_wavelet_data_CPU);
|
||||
else
|
||||
if(mode==GPU_DATA)
|
||||
ShowWaveletData(m_wavelet_data_GPU);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| ShowWaveletData |
|
||||
//+------------------------------------------------------------------+
|
||||
void CWavelet::ShowWaveletData(const float &m_wavelet_data[])
|
||||
{
|
||||
//--- calculate min/max and range
|
||||
int count=ArraySize(m_wavelet_data);
|
||||
float min_value=m_wavelet_data[0];
|
||||
float max_value=m_wavelet_data[0];
|
||||
for(int i=1; i<count; i++)
|
||||
{
|
||||
min_value=MathMin(min_value,m_wavelet_data[i]);
|
||||
max_value=MathMax(max_value,m_wavelet_data[i]);
|
||||
}
|
||||
float range=max_value-min_value;
|
||||
if(range>0)
|
||||
{
|
||||
for(int j=0; j<m_ysize; j++)
|
||||
{
|
||||
for(int i=0; i<m_xsize; i++)
|
||||
{
|
||||
int pos=j*m_xsize+i;
|
||||
int colindex=int(m_maxcolor*(m_wavelet_data[pos]-min_value)/range);
|
||||
m_bmp_buffer[pos]=GetPalColor(colindex);
|
||||
}
|
||||
}
|
||||
//--- show image
|
||||
ResourceCreate(m_res_name,m_bmp_buffer,m_xsize,m_ysize,0,0,0,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
ChartRedraw();
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Weirstrass function |
|
||||
//+------------------------------------------------------------------+
|
||||
float Weirstrass(float x,float a,float b)
|
||||
{
|
||||
float sum=0.0;
|
||||
float b0=b;
|
||||
float a0=a;
|
||||
for(int n=0; n<35; n++)
|
||||
{
|
||||
float v=b0*(float)MathCos(a0*M_PI*x);
|
||||
sum=sum+v;
|
||||
a0=a0*a;
|
||||
b0=b0*b;
|
||||
}
|
||||
return(sum);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareModelData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PrepareModelData(float &price_data[],const int datacount)
|
||||
{
|
||||
ArrayResize(price_data,datacount);
|
||||
//--- Weirstrass function
|
||||
float x1=0;
|
||||
float x2=2;
|
||||
float dx=(x2-x1)/datacount;
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
price_data[i]=Weirstrass(x1+dx*i,(float)3,(float)0.62);
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PreparePriceData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PreparePriceData(const string symbol,ENUM_TIMEFRAMES timeframe,float &price_data[],const int datacount)
|
||||
{
|
||||
ArrayResize(price_data,datacount);
|
||||
|
||||
double price_data_double[];
|
||||
CopyClose(symbol,timeframe,0,datacount,price_data_double);
|
||||
|
||||
int size=ArraySize(price_data_double);
|
||||
for(int i=0; i<size; i++)
|
||||
{
|
||||
price_data[i]=(float)price_data_double[i];
|
||||
}
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| PrepareMomentumData |
|
||||
//+------------------------------------------------------------------+
|
||||
void PrepareMomentumData(float &price_data[],float &momentum_data[],const int momentum_period)
|
||||
{
|
||||
int size=ArraySize(price_data);
|
||||
int datacount=size-momentum_period;
|
||||
//---
|
||||
ArrayResize(momentum_data,datacount);
|
||||
for(int i=0; i<datacount; i+=1)
|
||||
{
|
||||
momentum_data[i]=price_data[i+momentum_period]-price_data[i];
|
||||
}
|
||||
ArrayCopy(price_data,price_data,momentum_period,0,datacount);
|
||||
ArrayResize(price_data,datacount);
|
||||
//--- rescale momentum data
|
||||
float min_value=momentum_data[0];
|
||||
float max_value=momentum_data[0];
|
||||
for(int i=1; i<datacount; i++)
|
||||
{
|
||||
float value=momentum_data[i];
|
||||
if(momentum_data[i]>max_value)
|
||||
max_value=value;
|
||||
|
||||
if(momentum_data[i]<min_value)
|
||||
min_value=value;
|
||||
}
|
||||
float range=max_value-min_value;
|
||||
for(int i=0; i<datacount; i+=1)
|
||||
momentum_data[i]=-1+2*(momentum_data[i]-min_value)/range;
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
//---
|
||||
int momentum_period=8;
|
||||
float price_data[];
|
||||
float momentum_data[];
|
||||
PrepareModelData(price_data,SIZE_X+momentum_period);
|
||||
//PreparePriceData("EURUSD",PERIOD_M1,price_data,SIZE_X+momentum_period);
|
||||
PrepareMomentumData(price_data,momentum_data,momentum_period);
|
||||
|
||||
double price_data_double[];
|
||||
double momentum_data_double[];
|
||||
|
||||
int datacount=ArraySize(price_data);
|
||||
ArrayResize(price_data_double,datacount);
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
price_data_double[i]=(double)price_data[i];
|
||||
}
|
||||
datacount=ArraySize(momentum_data);
|
||||
ArrayResize(momentum_data_double,datacount);
|
||||
for(int i=0; i<datacount; i++)
|
||||
{
|
||||
momentum_data_double[i]=(double)momentum_data[i];
|
||||
}
|
||||
|
||||
CGraphic graph_price;
|
||||
CGraphic graph_momentum;
|
||||
graph_price.Create(0,"price",0,0,0,SIZE_X+130+6,SIZE_Y);
|
||||
graph_price.XAxis().MaxGrace(0);
|
||||
graph_price.HistorySymbolSize(10);
|
||||
graph_price.CurveAdd(price_data_double,ColorToARGB(clrRed,255),CURVE_LINES,"Price");
|
||||
graph_price.CurvePlotAll();
|
||||
graph_price.Redraw(true);
|
||||
graph_price.Update();
|
||||
//---
|
||||
graph_momentum.Create(0,"momentum",0,0,SIZE_Y,SIZE_X+130+6,SIZE_Y+SIZE_Y);
|
||||
graph_momentum.XAxis().MaxGrace(0);
|
||||
graph_momentum.HistorySymbolSize(10);
|
||||
graph_momentum.CurveAdd(momentum_data_double,ColorToARGB(clrBlue,255),CURVE_LINES,"Momentum");
|
||||
graph_momentum.CurvePlotAll();
|
||||
graph_momentum.Redraw(true);
|
||||
graph_momentum.Update();
|
||||
//---
|
||||
uint time_cpu=0;
|
||||
CWavelet wavelet;
|
||||
wavelet.Create("Wavelet",50,2*SIZE_Y,SIZE_X,SIZE_Y);
|
||||
if(!wavelet.CalculateWavelet_CPU(momentum_data,time_cpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on CPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
//wavelet.ShowWavelet(CPU_DATA);
|
||||
uint time_gpu=0;
|
||||
if(!wavelet.CalculateWavelet_GPU(momentum_data,time_gpu))
|
||||
{
|
||||
PrintFormat("Error in calculation on GPU. Error code=%d",GetLastError());
|
||||
return;
|
||||
}
|
||||
wavelet.ShowWavelet(GPU_DATA);
|
||||
//---
|
||||
double CPU_GPU_ratio=0;
|
||||
if(time_gpu!=0)
|
||||
CPU_GPU_ratio=1.0*time_cpu/time_gpu;
|
||||
//---
|
||||
PrintFormat("time CPU=%d ms, time GPU=%d ms, CPU/GPU ratio: %f",time_cpu,time_gpu,CPU_GPU_ratio);
|
||||
//--- Sleep 10 seconds
|
||||
Sleep(10000);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Reference in New Issue
Block a user