cuda: host

/*
* 电子科技大学 电子工程学院 张舒 编写,未经允许不得用于商业用途
* BP nerve neworl research
* Host code.
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cutil.h>
#include "cublas.h"
#include <nerve_kernel.cu>

#define size_SamInEx  ((int)(SamNum*(InDim+1)))
#define size_SamOut  ((int)(SamNum*OutDim))
#define size_W1Ex   ((int)(HiddenUnitNum*(InDim+1)))
#define size_W2Ex   ((int)(OutDim*(HiddenUnitNum+1)))


int iDivUp(int, int);
void InitSample(int, int, int, int, float*, float*, float*, float*);
void train(float, int, int, int, int, int, float*, float*, float*, float*);
void logsig(float*, float*, float*, int, int);

extern "C"
////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////

int
main(int argc, char** argv)
{
const float lr = 0.1f;
const int MaxEpochs = 20000;
const int HiddenUnitNum = 10;
const int InDim = 2;
const int OutDim = 3;
const int SamNum = 200;

float* h_SamInEx;
float* h_SamOut;
float* h_W1Ex;
float* h_W2Ex;

CUT_DEVICE_INIT();
cublasInit();

CUDA_SAFE_CALL( cudaMallocHost((void**) &h_SamInEx, sizeof(float)*size_SamInEx));
CUDA_SAFE_CALL( cudaMallocHost((void**) &h_SamOut, sizeof(float)*size_SamOut));
CUDA_SAFE_CALL( cudaMallocHost((void**) &h_W1Ex, sizeof(float)*size_W1Ex));
CUDA_SAFE_CALL( cudaMallocHost((void**) &h_W2Ex, sizeof(float)*size_W2Ex));

InitSample(SamNum, InDim, OutDim, HiddenUnitNum, h_SamInEx, h_SamOut, h_W1Ex, h_W2Ex);
    train(lr, MaxEpochs, HiddenUnitNum, InDim, OutDim, SamNum, h_SamInEx, h_SamOut, h_W1Ex, h_W2Ex);

printf("\n");
for(int i=0; i <size_W1Ex; i++)
{
printf("%3.3f ",h_W1Ex[i]);
}

printf("\n");
for(int i=0; i <size_W2Ex; i++)
{
printf("%3.3f ",h_W2Ex[i]);
}

FILE *p;
p = fopen("D:\\0_W1Ex.dat", "wb");
fwrite(h_W1Ex,sizeof(float),size_W1Ex,p);
fclose(p);

p = fopen("D:\\0_W2Ex.dat", "wb");
fwrite(h_W2Ex,sizeof(float),size_W2Ex,p);
fclose(p);


CUDA_SAFE_CALL( cudaFreeHost((h_SamInEx)));
CUDA_SAFE_CALL( cudaFreeHost((h_SamOut)));
CUDA_SAFE_CALL( cudaFreeHost((h_W1Ex)));
CUDA_SAFE_CALL( cudaFreeHost((h_W2Ex)));

h_SamInEx = NULL;
h_SamOut = NULL;
h_W1Ex=NULL;
h_W2Ex=NULL;

cublasShutdown();
CUT_EXIT(argc, argv);
}

int
iDivUp(int a, int b){
    return ((a % b) != 0) ? (a / b + 1) : (a / b);
}

void InitSample(int SamNum, int InDim, int OutDim, int HiddenUnitNum, float* h_SamInEx, float* h_SamOut, float* h_W1Ex,float* h_W2Ex)
{
srand(clock());
FILE *p;
p = fopen("D:\\samin.dat", "rb");

for(int i=0; i <size_SamInEx; i++)
{
h_SamInEx[i] = 1.0f;
}

for(int i=0; i <SamNum; i++)
{
for(int j=0; j <InDim;j++)
{
fread(&h_SamInEx[j*SamNum+i],sizeof(float),1,p);
}
}
fclose(p);

for(int i=0; i <size_SamInEx; i++)
{
printf("%3.4f\t",h_SamInEx[i]);
}

printf("\n");

p = fopen("D:\\samout.dat", "rb");
for(int i=0; i <SamNum; i++)
{
for(int j=0; j <OutDim;j++)
{
fread(&h_SamOut[j*SamNum+i],sizeof(float),1,p);
}
}
/* for(int i=0; i <size_SamOut; i++)
{
printf("%3.5f\t",h_SamOut[i]);
}*/
fclose(p);

printf("\n");

for(int i=0; i <size_W1Ex; i++)
{
h_W1Ex[i]=0.2f*rand()/(float)RAND_MAX - 0.1f;
}
for(int i=0; i <size_W2Ex; i++)
{
h_W2Ex[i]=0.2f*rand()/(float)RAND_MAX - 0.1f;
}
/*
p = fopen("D:\\W1Ex.dat", "rb");
for(int i=0; i <(InDim+1); i++)
{
for(int j=0; j <HiddenUnitNum;j++)
{
fread(&h_W1Ex[j*(InDim+1)+i],sizeof(float),1,p);
}
}

for(int i=0; i <size_W1Ex; i++)
{
printf("%3.3f ",h_W1Ex[i]);
}
fclose(p);

printf("\n");

p = fopen("D:\\W2Ex.dat", "rb");
for(int i=0; i <(HiddenUnitNum+1); i++)
{
for(int j=0; j <OutDim;j++)
{
fread(&h_W2Ex[j*(HiddenUnitNum+1)+i],sizeof(float),1,p);
}
}

for(int i=0; i <size_W2Ex; i++)
{
printf("%3.3f ",h_W2Ex[i]);
}
fclose(p);
*/
printf("\n");

p=NULL;
}


void
train(float lr, int MaxEpochs, int HiddenUnitNum, int InDim,int OutDim,int SamNum, float* h_SamInEx, float* h_SamOut, float* h_W1Ex, float* h_W2Ex)
{
    float* d_W1Ex;
float* d_W2Ex;
float* d_W2;
float* d_SamInEx;
float* d_SamOut;
float* d_error;
float* HiddenOutEx;
float* NetworkOut;
float* Delta1;
float* Delta2;

unsigned int timer = 0;

    dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
    dim3 HiddenOutExgrid(iDivUp(SamNum, BLOCK_SIZE), iDivUp((HiddenUnitNum+1), BLOCK_SIZE));
dim3 HiddenOutgrid(iDivUp(SamNum, BLOCK_SIZE), iDivUp((HiddenUnitNum), BLOCK_SIZE));
dim3 SamOutgrid(iDivUp(SamNum, BLOCK_SIZE), iDivUp(OutDim, BLOCK_SIZE));

CUDA_SAFE_CALL(cudaMalloc((void**) &d_error, sizeof(float)*size_SamOut));
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_W1Ex, sizeof(float)*size_W1Ex));
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_W2Ex, sizeof(float)*size_W2Ex));
CUDA_SAFE_CALL(cudaMalloc((void**) &d_W2, sizeof(float)*OutDim*HiddenUnitNum));
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_SamInEx, sizeof(float)*size_SamInEx));
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_SamOut, sizeof(float)*size_SamOut));
CUDA_SAFE_CALL(cudaMalloc((void**) &HiddenOutEx, sizeof(float)*(HiddenUnitNum+1)*SamNum));
CUDA_SAFE_CALL(cudaMalloc((void**) &NetworkOut, sizeof(float)*size_SamOut));
CUDA_SAFE_CALL(cudaMalloc((void**) &Delta1, sizeof(float)*HiddenUnitNum*SamNum));
CUDA_SAFE_CALL(cudaMalloc((void**) &Delta2, sizeof(float)*size_SamOut)); 
    CUDA_SAFE_CALL(cudaMemcpy(d_W1Ex, h_W1Ex, sizeof(float)*size_W1Ex,cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(d_W2Ex, h_W2Ex, sizeof(float)*size_W2Ex,cudaMemcpyHostToDevice));
CUDA_SAFE_CALL(cudaMemcpy(d_SamInEx, h_SamInEx, sizeof(float)*size_SamInEx,cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(d_SamOut, h_SamOut, sizeof(float)*size_SamOut,cudaMemcpyHostToDevice));

CUT_SAFE_CALL( cutCreateTimer( &timer));
    CUT_SAFE_CALL( cutStartTimer( timer));

for(int l=0; l < MaxEpochs; l++)
{

for(int i=0; i < OutDim; i++)
{
CUDA_SAFE_CALL(cudaMemcpy(d_W2 + HiddenUnitNum * i, d_W2Ex + (HiddenUnitNum + 1) * i, sizeof(float)*HiddenUnitNum,cudaMemcpyDeviceToDevice));
}

cublasSgemm('n','n',SamNum, HiddenUnitNum, (InDim+1), 1.0f,  d_SamInEx, SamNum, d_W1Ex, (InDim+1), 0.0f, HiddenOutEx, SamNum );
logsig1 < < <HiddenOutExgrid, threads>>>(HiddenOutEx, SamNum, HiddenUnitNum);
cublasSgemm('n','n',SamNum, OutDim, (HiddenUnitNum+1), 1.0f,  HiddenOutEx, SamNum, d_W2Ex, (HiddenUnitNum+1), 0.0f, NetworkOut, SamNum );
logsig2 < < <SamOutgrid,threads>>>(NetworkOut, SamNum, OutDim);
dotsub < < <SamOutgrid,threads>>>(Delta2, d_SamOut, NetworkOut, SamNum, OutDim);
getdelta < < <SamOutgrid,threads>>>(Delta2, NetworkOut, SamNum, OutDim);
cublasSgemm('t','n',(HiddenUnitNum+1), OutDim, SamNum, lr,  HiddenOutEx, SamNum, Delta2, SamNum, 1.0f, d_W2Ex, (HiddenUnitNum+1) );
cublasSgemm('n','t', SamNum, HiddenUnitNum, OutDim, 1.0f,  Delta2, SamNum, d_W2, HiddenUnitNum, 0.0f, Delta1, SamNum );
getdelta < < <HiddenOutgrid, threads>>>(Delta1, HiddenOutEx, SamNum, HiddenUnitNum);
cublasSgemm('t','n', (InDim+1), HiddenUnitNum, SamNum, lr, d_SamInEx, SamNum, Delta1, SamNum, 1.0f, d_W1Ex, (InDim+1));
}
CUT_SAFE_CALL( cutStopTimer( timer));
    printf( "Processing time: %f (ms)\n", cutGetTimerValue( timer));
    CUT_SAFE_CALL( cutDeleteTimer( timer));
CUDA_SAFE_CALL(cudaMemcpy(h_W1Ex, d_W1Ex, sizeof(float)*size_W1Ex,cudaMemcpyDeviceToHost));
CUDA_SAFE_CALL(cudaMemcpy(h_W2Ex, d_W2Ex, sizeof(float)*size_W2Ex,cudaMemcpyDeviceToHost));

    CUDA_SAFE_CALL(cudaFree(d_SamInEx));
    CUDA_SAFE_CALL(cudaFree(d_SamOut));
    CUDA_SAFE_CALL(cudaFree(d_W1Ex));
CUDA_SAFE_CALL(cudaFree(d_W2Ex));
CUDA_SAFE_CALL(cudaFree(d_SamOut));
CUDA_SAFE_CALL(cudaFree(HiddenOutEx));
CUDA_SAFE_CALL(cudaFree(NetworkOut));
CUDA_SAFE_CALL(cudaFree(Delta1));
CUDA_SAFE_CALL(cudaFree(Delta2));
} 

/*
* BP nerve neworl research
* Device code.
*/

#ifndef _MATRIXMUL_KERNEL_H_
#define _MATRIXMUL_KERNEL_H_

#include <stdio.h>

#define BLOCK_SIZE 16

__global__ void
logsig1( float* A, int wA, int hA)
{
    const int posy = blockIdx.y * blockDim.y + threadIdx.y;
const int posx = blockIdx.x * blockDim.x + threadIdx.x;
const int pos = posy * wA + posx;

if(posx < wA )
{
if(posy < hA )
    A[pos] = 1/(1+expf(-A[pos]));
else if(posy == hA)
A[pos] = 1.0f;
}
}

__global__ void
logsig2( float* A, int wA, int hA)
{
    const int posy = blockIdx.y * blockDim.y + threadIdx.y;
const int posx = blockIdx.x * blockDim.x + threadIdx.x;
const int pos = posy * wA + posx;

if((posx < wA)&&(posy < hA))
{
    A[pos] = 1/(1+expf(-A[pos]));
}
}

__global__ void
dotsub( float* C, float* A, float* B, int wA, int hA)
{
    const int posy = blockIdx.y * blockDim.y + threadIdx.y;
const int posx = blockIdx.x * blockDim.x + threadIdx.x;
const int pos = posy * wA + posx;

if((posx < wA)&&(posy < hA))
{
    C[pos] = A[pos]-B[pos];
}
}

__global__ void
getdelta( float* C, float* A, int wA, int hA)
{
    const int tidx = threadIdx.x;
const int tidy = threadIdx.y;
const int posy = blockIdx.y * blockDim.y + threadIdx.y;
const int posx = blockIdx.x * blockDim.x + threadIdx.x;
const int pos = posy * wA + posx;

//__shared__ s_A[BLOCK_SIZE][BLOCK_SIZE];

if((posx < wA)&&(posy < hA))
{
//s_A[tidy][tidx] = A[pos];
    C[pos] *= ((1.0f - A[pos])*A[pos]);
}
__syncthreads();
}

#endif 

运行后用于评估结果的matlab程序：
注意编译时工程要加入依赖项cublas.lib，不能用emu模式编译


fid=fopen('d:\0_W1Ex.dat','rb')
W1Ex=fread(fid,[3,11],'float32')
fid=fopen('d:\0_W2Ex.dat','rb')
W2Ex=fread(fid,[11,3],'float32')
fclose(fid)
W1Ex=W1Ex'
W2Ex=W2Ex'


InDim = 2;
OutDim = 3;
HiddenUnitNum=10;
W1=W1Ex(:,1:InDim);
B1=W1Ex(:,InDim+1);
W2=W2Ex(:,1:HiddenUnitNum);
B2=W2Ex(:,1+HiddenUnitNum);
TestSamNum=5000;
TestSamIn=(rand(2,TestSamNum)-0.5)*4;

TestHiddenOut=logsig(W1*TestSamIn+repmat(B1,1,TestSamNum));
TestNetworkOut=logsig(W2*TestHiddenOut + repmat(B2,1,TestSamNum));
[Val,NNClass]=max(TestNetworkOut);

TestTargetOut=[];
for i=1:TestSamNum;
Sam=TestSamIn(:,i);
x=Sam(1,1);
y=Sam(2,1);
if((x>-1)&(x <1))==1
if((y>x/2+0.5)&(y <1))==1
TestTargetOut=[TestTargetOut 2];
elseif((y <-0.5)&(y>-1.5))==1;
TestTargetOut=[TestTargetOut 3];
else
TestTargetOut=[TestTargetOut 1];
end
else
TestTargetOut=[TestTargetOut 1];
end
end

NNC1Flag=abs(NNClass-1) <0.1;
NNC2Flag=abs(NNClass-2) <0.1;
NNC3Flag=abs(NNClass-3) <0.1;

TargetC1Flag=abs(TestTargetOut-1) <0.1;
TargetC2Flag=abs(TestTargetOut-2) <0.1;
TargetC3Flag=abs(TestTargetOut-3) <0.1;

Test_C1_num=sum(NNC1Flag)
Test_C2_num=sum(NNC2Flag)
Test_C3_num=sum(NNC3Flag)

Test_C1_C1=1.0*NNC1Flag*TargetC1Flag'
Test_C1_C2=1.0*NNC1Flag*TargetC2Flag'
Test_C1_C3=1.0*NNC1Flag*TargetC3Flag'

Test_C2_c1=1.0*NNC2Flag*TargetC1Flag'
Test_C2_C2=1.0*NNC2Flag*TargetC2Flag'
Test_C2_C3=1.0*NNC2Flag*TargetC3Flag'

Test_C3_C1=1.0*NNC3Flag*TargetC1Flag'
Test_C3_C2=1.0*NNC3Flag*TargetC2Flag'
Test_C3_C3=1.0*NNC3Flag*TargetC3Flag'

Test_Correct=(Test_C1_C1+Test_C2_C2+Test_C3_C3)/TestSamNum