Init mql5 project
This commit is contained in:
Bell
Binary file not shown.
@@ -0,0 +1,239 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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];
|
||||
uint passes_total=0;
|
||||
uint stages_total=0;
|
||||
|
||||
global_size[0]=data_count>>1;
|
||||
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 with %d size. Error code=%d",data_count,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/2/sizeof(double))/MathLog(2));
|
||||
int pwr=(int)MathMax(15,pwr_max-4);
|
||||
//--- prepare array and generate random data
|
||||
data_count=(int)MathPow(2,pwr);
|
||||
|
||||
if(data_count<4096)
|
||||
data_count=4096;
|
||||
|
||||
if(data_count>4*1024*1024)
|
||||
data_count=4*1024*1024;
|
||||
|
||||
Print(data_count," elements in double array");
|
||||
//---
|
||||
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()
|
||||
{
|
||||
COpenCL OpenCL;
|
||||
//--- 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;
|
||||
}
|
||||
//--- prepare array with random values
|
||||
double data_cpu[];
|
||||
int data_count=0;
|
||||
|
||||
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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,354 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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 1024
|
||||
#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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| |
|
||||
//+------------------------------------------------------------------+
|
||||
bool ExecutionWait(COpenCL& OpenCL,int kernel_index)
|
||||
{
|
||||
for(int i=0;;)
|
||||
{
|
||||
ENUM_OPENCL_EXECUTION_STATUS status=OpenCL.ExecutionStatus(kernel_index);
|
||||
if(status==CL_COMPLETE)
|
||||
{
|
||||
Print("i=",i);
|
||||
return(true);
|
||||
}
|
||||
if(status<0 || IsStopped())
|
||||
break;
|
||||
if(++i%10000==0)
|
||||
{
|
||||
Print("running... i=",i," execute_status=",status);
|
||||
Sleep(1);
|
||||
}
|
||||
}
|
||||
//---
|
||||
return(false);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| 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 workgroup_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);
|
||||
local_mem_size/=3;
|
||||
local_mem_size*=2;
|
||||
local_mem_size/=num_points;
|
||||
local_mem_size*=num_points;
|
||||
Print("local_mem_size=",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)*workgroup_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]=workgroup_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);
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,0))
|
||||
return(false);
|
||||
Print("fft_init passed");
|
||||
//-- further stages of the FFT
|
||||
if(num_points>points_per_group)
|
||||
{
|
||||
Print("num_points=",num_points," points_per_group=",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);
|
||||
}
|
||||
Print("fft_stage ",stage," passed");
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,1))
|
||||
return(false);
|
||||
}
|
||||
//--- scale values if performing the inverse FFT
|
||||
if(direction<0)
|
||||
{
|
||||
Print("direction=",direction);
|
||||
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);
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,2))
|
||||
return(false);
|
||||
Print("fft_scale passed");
|
||||
}
|
||||
//--- 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);
|
||||
}
|
||||
Print("buffer read");
|
||||
//--- 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();
|
||||
Print("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;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,225 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,419 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Wavelet.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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)<m_palette.Size())
|
||||
{
|
||||
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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,232 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| BitonicSort.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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];
|
||||
uint passes_total=0;
|
||||
uint stages_total=0;
|
||||
|
||||
global_size[0]=data_count>>1;
|
||||
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 with %d size. Error code=%d",data_count,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,float &data[],int &data_count)
|
||||
{
|
||||
int pwr_max=(int)(MathLog(global_memory_size/2/sizeof(float))/MathLog(2));
|
||||
int pwr=(int)MathMax(15,pwr_max-4);
|
||||
//--- prepare array and generate random data
|
||||
data_count=(int)MathPow(2,pwr);
|
||||
|
||||
if(data_count<4096)
|
||||
data_count=4096;
|
||||
|
||||
if(data_count>4*1024*1024)
|
||||
data_count=4*1024*1024;
|
||||
|
||||
Print(data_count," elements in float array");
|
||||
//---
|
||||
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()
|
||||
{
|
||||
COpenCL OpenCL;
|
||||
//--- 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;
|
||||
}
|
||||
//--- prepare array with random values
|
||||
float data_cpu[];
|
||||
int data_count;
|
||||
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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,349 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| FFT.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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 1024
|
||||
#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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| |
|
||||
//+------------------------------------------------------------------+
|
||||
bool ExecutionWait(COpenCL& OpenCL,int kernel_index)
|
||||
{
|
||||
for(int i=0;;)
|
||||
{
|
||||
ENUM_OPENCL_EXECUTION_STATUS status=OpenCL.ExecutionStatus(kernel_index);
|
||||
if(status==CL_COMPLETE)
|
||||
{
|
||||
Print("i=",i);
|
||||
return(true);
|
||||
}
|
||||
if(status<0 || IsStopped())
|
||||
break;
|
||||
if(++i%10000==0)
|
||||
{
|
||||
Print("running... i=",i," execute_status=",status);
|
||||
Sleep(1);
|
||||
}
|
||||
}
|
||||
//---
|
||||
return(false);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| 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 workgroup_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);
|
||||
local_mem_size/=3;
|
||||
local_mem_size*=2;
|
||||
local_mem_size/=num_points;
|
||||
local_mem_size*=num_points;
|
||||
Print("local_mem_size=",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)*workgroup_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]=workgroup_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);
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,0))
|
||||
return(false);
|
||||
Print("fft_init passed");
|
||||
//-- further stages of the FFT
|
||||
if(num_points>points_per_group)
|
||||
{
|
||||
Print("num_points=",num_points," points_per_group=",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);
|
||||
}
|
||||
Print("fft_stage ",stage," passed");
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,1))
|
||||
return(false);
|
||||
}
|
||||
//--- scale values if performing the inverse FFT
|
||||
if(direction<0)
|
||||
{
|
||||
Print("direction=",direction);
|
||||
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);
|
||||
}
|
||||
if(!ExecutionWait(OpenCL,2))
|
||||
return(false);
|
||||
Print("fft_scale passed");
|
||||
}
|
||||
//--- 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);
|
||||
}
|
||||
Print("buffer read");
|
||||
//--- 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();
|
||||
Print("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;
|
||||
};
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,219 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| MatrixMult.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
@@ -0,0 +1,436 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Wavelet.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#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)<m_palette.Size())
|
||||
{
|
||||
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);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
@@ -0,0 +1,223 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| OpenCL code |
|
||||
//+------------------------------------------------------------------+
|
||||
//| Based on pixel shader by Alexander Alekseev aka TDM-2014 |
|
||||
//| https://www.shadertoy.com/view/Ms2SD1 |
|
||||
//+------------------------------------------------------------------+
|
||||
#define NUM_STEPS (int)8
|
||||
#define PI (float)3.1415
|
||||
#define EPSILON (float)1e-3
|
||||
#define EPSILON_NRM ((float)0.1/iResolution.x)
|
||||
#define ITER_GEOMETRY (int)3
|
||||
#define ITER_FRAGMENT (int)5
|
||||
#define SEA_HEIGHT (float)0.6
|
||||
#define SEA_CHOPPY (float)4.0
|
||||
#define SEA_SPEED (float)0.8
|
||||
#define SEA_FREQ (float)0.16
|
||||
#define SEA_BASE vec3(0.1,0.19,0.22)
|
||||
#define SEA_WATER_COLOR vec3(0.8,0.9,0.6)
|
||||
#define SEA_TIME ((float)1.0 + iGlobalTime * SEA_SPEED)
|
||||
|
||||
float max1(float v1,float v2) { if (v1>v2) return v1; return v2; }
|
||||
float dot2(float2 v1,float2 v2) { return v1.x*v2.x+v1.y*v2.y; }
|
||||
float dot3(float3 v1,float3 v2) { return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z; }
|
||||
float fract1(float p) { return p-floor(p); }
|
||||
float2 fract2(float2 p) { return p-floor(p); }
|
||||
float mix1(float x1,float x2,float a) { return(x1+(x2-x1)*a); }
|
||||
float2 mix2(float2 x1,float2 x2,float2 a) { float2 r={x1.x+(x2.x-x1.x)*a.x,x1.y+(x2.y-x1.y)*a.y }; return(r); }
|
||||
float3 mix3(float3 x1,float3 x2,float3 a) { float3 r={x1.x+(x2.x-x1.x)*a.x,x1.y+(x2.y-x1.y)*a.y,x1.z+(x2.z-x1.z)*a.z }; return(r); }
|
||||
float3 vec31(float x) { float3 r={x,x,x}; return(r); }
|
||||
float3 vec3(float x,float y,float z) { float3 r={x,y,z}; return(r); }
|
||||
float2 vec2(float x,float y) { float2 r={x,y}; return(r); }
|
||||
float4 vec4(float3 x,float y) { float4 r={x.x,x.y,x.z,y}; return(r); }
|
||||
float3 reflect3(float3 I,float3 N) { return(I-(float)2.0*dot3(N,I)*N); }
|
||||
float diffuse(float3 n,float3 l,float p) { return pow(dot3(n,l) * (float)0.4 + (float)0.6,p); }
|
||||
float2 abs2(float2 f) { if(f.x<0) f.x=-f.x; if(f.y<0) f.y=-f.y; return(f); }
|
||||
float hash(float2 p) { return fract1(sin(dot2(p,vec2((float)127.1,(float)311.7)))*43758.5453123); }
|
||||
float specular(float3 n,float3 l,float3 e,float s) { return pow(max1(dot3(reflect3(e,n),l),0.0),s)*((s+(float)8.0)/(float)(3.1415 * 8.0)); }
|
||||
|
||||
float3 getSkyColor(float3 e)
|
||||
{
|
||||
e.y=max1(e.y,0.0);
|
||||
float3 r2;
|
||||
r2.x = (float)1.0-e.y;
|
||||
r2.y = r2.x;
|
||||
r2.z = (float)0.6+r2.x *(float)0.4;
|
||||
r2.x = r2.x * r2.x;
|
||||
return r2;
|
||||
}
|
||||
|
||||
void fromEuler(float3 ang,float3 *r1,float3 *r2,float3 *r3)
|
||||
{
|
||||
float2 a1 = vec2(sin(ang.x),cos(ang.x));
|
||||
float2 a2 = vec2(sin(ang.y),cos(ang.y));
|
||||
float2 a3 = vec2(sin(ang.z),cos(ang.z));
|
||||
*r1 = vec3(a1.y*a3.y+a1.x*a2.x*a3.x,a1.y*a2.x*a3.x+a3.y*a1.x,-a2.y*a3.x);
|
||||
*r2 = vec3(-a2.y*a1.x,a1.y*a2.y,a2.x);
|
||||
*r3 = vec3(a3.y*a1.x*a2.x+a1.y*a3.x,a1.x*a3.x-a1.y*a3.y*a2.x,a2.y*a3.y);
|
||||
}
|
||||
|
||||
float noise(float2 p)
|
||||
{
|
||||
float2 i = floor(p);
|
||||
float2 f = fract2(p);
|
||||
float2 u = f * f * ((float)3.0 - (float)2.0 * f);
|
||||
float mx1=mix1(hash(i+vec2(0.0,0.0)),hash(i+vec2(1.0,0.0)),u.x);
|
||||
float mx2=mix1(hash(i+vec2(0.0,1.0)),hash(i+vec2(1.0,1.0)),u.x);
|
||||
return (float)-1.0 + (float)2.0*mix(mx1,mx2,u.y);
|
||||
}
|
||||
|
||||
float sea_octave(float2 uv, float choppy)
|
||||
{
|
||||
uv += noise(uv);
|
||||
float2 wv = (float)1.0-abs2(sin(uv));
|
||||
float2 swv = abs2(cos(uv));
|
||||
wv = mix2(wv,swv,wv);
|
||||
return pow((float)1.0-pow(wv.x * wv.y,(float)0.65),choppy);
|
||||
}
|
||||
|
||||
float map(float3 p,float iGlobalTime)
|
||||
{
|
||||
float freq = SEA_FREQ;
|
||||
float amp = SEA_HEIGHT;
|
||||
float choppy = SEA_CHOPPY;
|
||||
float2 uv = p.xz; uv.x *= (float)0.75;
|
||||
float d, h = 0.0;
|
||||
for(int i = 0; i < ITER_GEOMETRY; i++)
|
||||
{
|
||||
d = sea_octave((uv+SEA_TIME)*freq,choppy);
|
||||
d += sea_octave((uv-SEA_TIME)*freq,choppy);
|
||||
h += d * amp;
|
||||
float2 uvt={(float)1.6*uv.x+(float)1.2*uv.y, (float)-1.2*uv.x+(float)1.6*uv.y };
|
||||
uv = uvt;
|
||||
freq *= (float)1.9; amp *= (float)0.22;
|
||||
choppy = mix1(choppy,(float)1.0,(float)0.2);
|
||||
}
|
||||
return p.y - h;
|
||||
}
|
||||
|
||||
float map_detailed(float3 p,float iGlobalTime)
|
||||
{
|
||||
float freq = SEA_FREQ;
|
||||
float amp = SEA_HEIGHT;
|
||||
float choppy = SEA_CHOPPY;
|
||||
float2 uv = p.xz; uv.x *= (float)0.75;
|
||||
float d, h = 0.0;
|
||||
for(int i = 0; i < ITER_FRAGMENT; i++)
|
||||
{
|
||||
d =sea_octave((uv+SEA_TIME)*freq,choppy);
|
||||
d+=sea_octave((uv-SEA_TIME)*freq,choppy);
|
||||
h+=d * amp;
|
||||
float2 uvt={ (float)1.6*uv.x+(float)1.2*uv.y, (float)-1.2*uv.x+(float)1.6*uv.y };
|
||||
uv = uvt;
|
||||
freq *= (float)1.9; amp *= (float)0.22;
|
||||
choppy = mix1(choppy,(float)1.0,(float)0.2);
|
||||
}
|
||||
return p.y - h;
|
||||
}
|
||||
|
||||
float3 getSeaColor(float3 p, float3 n, float3 l, float3 eye, float3 dist)
|
||||
{
|
||||
float fresnel = clamp((float)1.0 - dot3(n,-eye), (float)0.0, (float)1.0);
|
||||
fresnel = pow(fresnel,(float)3.0) * (float)0.65;
|
||||
float3 reflected = getSkyColor(reflect3(eye,n));
|
||||
float3 refracted = SEA_BASE + diffuse(n,l,(float)80.0) * SEA_WATER_COLOR * (float)0.12;
|
||||
float3 color = mix3(refracted,reflected,fresnel);
|
||||
float atten = max1((float)1.0 - dot(dist,dist) * (float)0.001, (float)0.0);
|
||||
color += SEA_WATER_COLOR * (p.y - SEA_HEIGHT) * (float)0.18 * atten;
|
||||
color += vec31(specular(n,l,eye,(float)60.0));
|
||||
|
||||
if(isnan(color.x))
|
||||
color.x=0.0;
|
||||
if(isnan(color.y))
|
||||
color.y=0.0;
|
||||
if(isnan(color.z))
|
||||
color.z=0.0;
|
||||
|
||||
return color;
|
||||
}
|
||||
|
||||
float3 getNormal(float3 p, float eps,float iGlobalTime)
|
||||
{
|
||||
float3 n;
|
||||
n.y = map_detailed(p,iGlobalTime);
|
||||
n.x = map_detailed(vec3(p.x+eps,p.y,p.z),iGlobalTime) - n.y;
|
||||
n.z = map_detailed(vec3(p.x,p.y,p.z+eps),iGlobalTime) - n.y;
|
||||
n.y = eps;
|
||||
return normalize(n);
|
||||
}
|
||||
|
||||
float heightMapTracing(float3 ori, float3 dir,float3 *p,float iGlobalTime)
|
||||
{
|
||||
float tm = (float)0.0;
|
||||
float tx = (float)1000.0;
|
||||
float hx = map(ori + dir * tx,iGlobalTime);
|
||||
if(hx > (float)0.0) return tx;
|
||||
float hm = map(ori + dir * tm,iGlobalTime);
|
||||
float tmid = (float)0.0;
|
||||
for(int i = 0; i < NUM_STEPS; i++)
|
||||
{
|
||||
tmid = mix(tm,tx, hm/(hm-hx));
|
||||
*p = ori + dir * tmid;
|
||||
float hmid = map(*p,iGlobalTime);
|
||||
if(hmid < (float)0.0)
|
||||
{
|
||||
tx = tmid;
|
||||
hx = hmid;
|
||||
}
|
||||
else
|
||||
{
|
||||
tm = tmid;
|
||||
hm = hmid;
|
||||
}
|
||||
}
|
||||
return tmid;
|
||||
}
|
||||
|
||||
float4 mainImage(float2 fragCoord,float iGlobalTime,float2 iResolution)
|
||||
{
|
||||
float4 fragColor;
|
||||
float2 iMouse= {0,0};
|
||||
float2 uv = fragCoord.xy / iResolution;
|
||||
uv = uv * (float)2.0 - (float)1.0;
|
||||
uv.x *= iResolution.x / iResolution.y;
|
||||
float time = iGlobalTime * (float)0.3 + iMouse.x*(float)0.01;
|
||||
float3 ang = vec3(sin(time*(float)3.0)*(float)0.1,sin(time)*(float)0.2+(float)0.3,time);
|
||||
float3 ori = vec3(0.0,3.5,time*(float)5.0);
|
||||
float3 dir = normalize(vec3(uv.x,uv.y,(float)-2.0)); dir.z += length(uv) * (float)0.15;
|
||||
dir = normalize(dir);
|
||||
float3 r1,r2,r3;
|
||||
fromEuler(ang,&r1,&r2,&r3);
|
||||
float3 r;
|
||||
r.x=r1.x*dir.x+r1.y*dir.y+r1.z*dir.z;
|
||||
r.y=r2.x*dir.x+r2.y*dir.y+r2.z*dir.z;
|
||||
r.z=r3.x*dir.x+r3.y*dir.y+r3.z*dir.z;
|
||||
dir=r;
|
||||
float3 p;
|
||||
heightMapTracing(ori,dir,&p,iGlobalTime);
|
||||
float3 dist = p - ori;
|
||||
float3 n = getNormal(p, dot(dist,dist) * EPSILON_NRM,iGlobalTime);
|
||||
float3 light = normalize(vec3((float)0.0,(float)1.0,(float)0.8));
|
||||
float3 seacol=getSeaColor(p,n,light,dir,dist);
|
||||
float3 color = mix(getSkyColor(dir),seacol,pow(smoothstep((float)0.0,(float)-0.05,dir.y),(float)0.3));
|
||||
fragColor = vec4(pow(color,vec31((float)0.75)),(float)1.0);
|
||||
return(fragColor);
|
||||
}
|
||||
|
||||
__kernel void Seascape(float iGlobalTime,__global uint *out)
|
||||
{
|
||||
size_t w = get_global_size(0);
|
||||
size_t h = get_global_size(1);
|
||||
float2 iRes = {(float)w,(float)h};
|
||||
size_t gx = get_global_id(0);
|
||||
size_t gy = get_global_id(1);
|
||||
float2 coord={gx,gy};
|
||||
float4 res=mainImage(coord,iGlobalTime,iRes);
|
||||
uint b=(uint)(res.z*255);
|
||||
uint g=(uint)(res.y*255);
|
||||
uint r=(uint)(res.x*255);
|
||||
if(b>255) b=255;
|
||||
if(r>255) r=255;
|
||||
if(g>255) g=255;
|
||||
out[w*((uint)(iRes.y-1)-gy)+gx] = (r<<16)|(g<<8)|b;
|
||||
};
|
||||
Binary file not shown.
@@ -0,0 +1,193 @@
|
||||
//+------------------------------------------------------------------+
|
||||
//| Seascape.mq5 |
|
||||
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
||||
//| https://www.mql5.com |
|
||||
//+------------------------------------------------------------------+
|
||||
#property copyright "Copyright 2000-2025, MetaQuotes Ltd."
|
||||
#property link "https://www.mql5.com"
|
||||
#property version "1.00"
|
||||
|
||||
#resource "seascape.cl" as string ExtCL
|
||||
|
||||
uint ExtSizeX=640; // Width
|
||||
uint ExtSizeY=360; // Height
|
||||
//+------------------------------------------------------------------+
|
||||
//| Initialize OpenCL engine |
|
||||
//+------------------------------------------------------------------+
|
||||
bool ModelInitialize(int &cl_ctx,int &cl_prg,int &cl_krn,int &cl_mem)
|
||||
{
|
||||
//--- initializing OpenCL objects
|
||||
if((cl_ctx=CLContextCreate())==INVALID_HANDLE)
|
||||
{
|
||||
Print("OpenCL not found");
|
||||
return(false);
|
||||
}
|
||||
//--- compile the program
|
||||
string build_log;
|
||||
if((cl_prg=CLProgramCreate(cl_ctx,ExtCL,build_log))==INVALID_HANDLE)
|
||||
{
|
||||
//--- free
|
||||
CLContextFree(cl_ctx);
|
||||
cl_ctx=INVALID_HANDLE;
|
||||
//---
|
||||
Print("OpenCL program create failed: ",build_log);
|
||||
return(false);
|
||||
}
|
||||
//--- kernel
|
||||
if((cl_krn=CLKernelCreate(cl_prg,"Seascape"))==INVALID_HANDLE)
|
||||
{
|
||||
//--- free
|
||||
CLProgramFree(cl_prg);
|
||||
cl_prg=INVALID_HANDLE;
|
||||
|
||||
CLContextFree(cl_ctx);
|
||||
cl_ctx=INVALID_HANDLE;
|
||||
//---
|
||||
Print("OpenCL kernel create failed");
|
||||
return(false);
|
||||
}
|
||||
//--- create buffer
|
||||
ExtSizeX=(uint)ChartGetInteger(0,CHART_WIDTH_IN_PIXELS);
|
||||
ExtSizeY=(uint)ChartGetInteger(0,CHART_HEIGHT_IN_PIXELS);
|
||||
|
||||
if((cl_mem=CLBufferCreate(cl_ctx,ExtSizeX*ExtSizeY*sizeof(uint),CL_MEM_READ_WRITE))==INVALID_HANDLE)
|
||||
{
|
||||
//--- free objects
|
||||
CLKernelFree(cl_krn);
|
||||
cl_krn=INVALID_HANDLE;
|
||||
|
||||
CLProgramFree(cl_prg);
|
||||
cl_prg=INVALID_HANDLE;
|
||||
|
||||
CLContextFree(cl_ctx);
|
||||
cl_ctx=INVALID_HANDLE;
|
||||
//---
|
||||
Print("OpenCL buffer create failed");
|
||||
return(false);
|
||||
}
|
||||
|
||||
CLSetKernelArgMem(cl_krn,1,cl_mem);
|
||||
//---
|
||||
return(true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Check and resize the window |
|
||||
//+------------------------------------------------------------------+
|
||||
bool ModelResize(const int cl_ctx,const int cl_krn,int &cl_mem)
|
||||
{
|
||||
uint width =(uint)ChartGetInteger(0,CHART_WIDTH_IN_PIXELS);
|
||||
uint height=(uint)ChartGetInteger(0,CHART_HEIGHT_IN_PIXELS);
|
||||
//--- correction
|
||||
width=(width + 7) & ~7; // align to 8
|
||||
if(width<8)
|
||||
width=8;
|
||||
|
||||
if(height<8)
|
||||
height=8;
|
||||
//--- the same size?
|
||||
if(ExtSizeX!=width || ExtSizeY!=height)
|
||||
{
|
||||
ExtSizeX=width;
|
||||
ExtSizeY=height;
|
||||
//--- create buffer
|
||||
if(cl_mem!=INVALID_HANDLE)
|
||||
CLBufferFree(cl_mem);
|
||||
|
||||
if((cl_mem=CLBufferCreate(cl_ctx,ExtSizeX*ExtSizeY*sizeof(uint),CL_MEM_READ_WRITE))!=INVALID_HANDLE)
|
||||
{
|
||||
CLSetKernelArgMem(cl_krn,1,cl_mem);
|
||||
return(true);
|
||||
}
|
||||
}
|
||||
//--- no changes
|
||||
return(false);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Shutdown OpenCL engine |
|
||||
//+------------------------------------------------------------------+
|
||||
void ModelShutdown(int cl_ctx,int cl_prg,int cl_krn,int cl_mem)
|
||||
{
|
||||
//--- remove OpenCL objects
|
||||
if(cl_mem!=INVALID_HANDLE)
|
||||
{
|
||||
CLBufferFree(cl_mem);
|
||||
cl_mem=INVALID_HANDLE;
|
||||
}
|
||||
|
||||
if(cl_krn!=INVALID_HANDLE)
|
||||
{
|
||||
CLKernelFree(cl_krn);
|
||||
cl_krn=INVALID_HANDLE;
|
||||
}
|
||||
|
||||
if(cl_prg!=INVALID_HANDLE)
|
||||
{
|
||||
CLProgramFree(cl_prg);
|
||||
cl_prg=INVALID_HANDLE;
|
||||
}
|
||||
|
||||
if(cl_ctx!=INVALID_HANDLE)
|
||||
{
|
||||
CLContextFree(cl_ctx);
|
||||
cl_ctx=INVALID_HANDLE;
|
||||
}
|
||||
//---
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
//| Script program start function |
|
||||
//+------------------------------------------------------------------+
|
||||
void OnStart()
|
||||
{
|
||||
int cl_ctx=INVALID_HANDLE,cl_prg=INVALID_HANDLE;
|
||||
int cl_krn=INVALID_HANDLE,cl_mem=INVALID_HANDLE;
|
||||
//--- prepare the chart
|
||||
ChartSetInteger(0,CHART_SHOW,false);
|
||||
ChartRedraw();
|
||||
//--- initializing OpenCL objects
|
||||
if(ModelInitialize(cl_ctx,cl_prg,cl_krn,cl_mem))
|
||||
{
|
||||
uint buf[];
|
||||
uint work[2];
|
||||
uint offset[2]={0,0};
|
||||
string objname="OpenCL_"+IntegerToString(ChartID());
|
||||
string resname="::Seascape_"+IntegerToString(ChartID());
|
||||
//--- create object and empty picture
|
||||
ObjectCreate(0,objname,OBJ_BITMAP_LABEL,0,0,0);
|
||||
ObjectSetInteger(0,objname,OBJPROP_XDISTANCE,0);
|
||||
ObjectSetInteger(0,objname,OBJPROP_YDISTANCE,0);
|
||||
|
||||
ArrayResize(buf,ExtSizeX*ExtSizeY);
|
||||
ResourceCreate(resname,buf,ExtSizeX,ExtSizeY,0,0,ExtSizeX,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
ObjectSetString(0,objname,OBJPROP_BMPFILE,resname);
|
||||
//--- render
|
||||
work[0]=ExtSizeX;
|
||||
work[1]=ExtSizeY;
|
||||
|
||||
for(float time=0;!IsStopped();time+=0.04f)
|
||||
{
|
||||
//--- check the resolution
|
||||
if(ModelResize(cl_ctx,cl_krn,cl_mem))
|
||||
{
|
||||
work[0]=ExtSizeX;
|
||||
work[1]=ExtSizeY;
|
||||
ArrayResize(buf,ExtSizeX*ExtSizeY);
|
||||
}
|
||||
//--- rendering the frame
|
||||
CLSetKernelArg(cl_krn,0,time);
|
||||
CLExecute(cl_krn,2,offset,work);
|
||||
//--- take the frame data, save in memory and draw it
|
||||
CLBufferRead(cl_mem,buf);
|
||||
ResourceCreate(resname,buf,ExtSizeX,ExtSizeY,0,0,ExtSizeX,COLOR_FORMAT_XRGB_NOALPHA);
|
||||
|
||||
ChartRedraw();
|
||||
Sleep(0);
|
||||
}
|
||||
//--- remove OpenCL objects
|
||||
ModelShutdown(cl_ctx,cl_prg,cl_krn,cl_mem);
|
||||
//--- remove object
|
||||
ObjectDelete(0,objname);
|
||||
}
|
||||
//---
|
||||
ChartSetInteger(0,CHART_SHOW,true);
|
||||
}
|
||||
//+------------------------------------------------------------------+
|
||||
Binary file not shown.
Reference in New Issue
Block a user