Pet整理

1. UDP通信

1.1 UdpClient

#include <WinSock2.h>
#include <Ws2tcpip.h>
#include <stdio.h>
#include <windows.h>

void udpClient(){
	//加载套接字库（必须）
	WORD wVersionRequested;
	WSADATA wsaData;
	int err;
	wVersionRequested = MAKEWORD(1, 1);  //版本1.1
	err = WSAStartup(wVersionRequested, &wsaData);  // 加载套接字库需要指定版本
	
	/*
	AF_UNIX（本机通信）
	AF_INET（TCP/IP – IPv4）
	AF_INET6（TCP/IP – IPv6）
	
	SOCK_STREAM（TCP流）
	SOCK_DGRAM（UDP数据报）
	SOCK_RAW（原始套接字）
	*/
	
	if (err != 0)
	{
		printf("WSAStartup failed with error: %d\n", err);
		return 1;
	}
	//低字节，高字节
	if (LOBYTE(wsaData.wVersion) != 1 || HIBYTE(wsaData.wVersion) != 1)
	{
		printf("Could not find a usable version of Winsock.dll\n");
		WSACleanup();
	}
	else
		printf("The Winsock 1.1 dll was found okay\n");

	//创建套接字
	SOCKET sockSrv = socket(AF_INET, SOCK_DGRAM, 0);

	SOCKADDR_IN addrSrv;
	addrSrv.sin_addr.S_un.S_addr = inet_addr("127.0.0.1");
	//inet_pton(AF_INET, "115.156.129.14", &addrClient[0].sin_addr.S_un.S_addr);
	addrSrv.sin_family = AF_INET;
	addrSrv.sin_port = htons(6000);

	//发送
	sendto(sockSrv, "client", strlen("client") + 1, 0, (SOCKADDR*)&addrSrv, sizeof(SOCKADDR));

	closesocket(sockSrv);
	WSACleanup();
}

1.2 UdpServer

void udpServer(){

	//加载套接字库（必须）
	WORD wVersionRequested;
	WSADATA wsaData;
	int err;
	wVersionRequested = MAKEWORD(1, 1);  //版本1.1
	err = WSAStartup(wVersionRequested, &wsaData);  // 加载套接字库需要指定版本
	
	/*
	AF_UNIX（本机通信）
	AF_INET（TCP/IP – IPv4）
	AF_INET6（TCP/IP – IPv6）
	
	SOCK_STREAM（TCP流）
	SOCK_DGRAM（UDP数据报）
	SOCK_RAW（原始套接字）
	*/
	
	if (err != 0)
	{
		printf("WSAStartup failed with error: %d\n", err);
		return 1;
	}
	//低字节，高字节
	if (LOBYTE(wsaData.wVersion) != 1 || HIBYTE(wsaData.wVersion) != 1)
	{
		printf("Could not find a usable version of Winsock.dll\n");
		WSACleanup();
	}
	else
		printf("The Winsock 1.1 dll was found okay\n");

	SOCKADDR_IN addrSrv;
	addrSrv.sin_addr.S_un.S_addr = inet_addr("115.156.131.120");    //inet_addr("115.156.131.120");          // 表示32位的IP地址
	addrSrv.sin_family = AF_INET;           // 表示地址类型
	addrSrv.sin_port = htons(888);         // 表示端口号
	
	//套接字，地址
	bind(sockSrv, (SOCKADDR*)&addrSrv, sizeof(SOCKADDR));           // 将套接字与地址绑定

	SOCKADDR_IN addrClient;
	int len = sizeof(SOCKADDR);
	char buff[100];
	
	//接收
	recvfrom(sockSrv, buff, 100, 0, (SOCKADDR*)&addrClient, &len);
	printf("%s\n", buff);
	
	closesocket(sockSrv);
	WSACleanup();
	
}

1.3 UdpSocket.h

#ifndef __UUDPSOCKET_H
#define __UUDPSOCKET_H

#include "TypeDef.h"
#include <cstddef>

class UdpSocket
{
public:
    UdpSocket();
    virtual ~UdpSocket();

    bool connect(const char *pHost, uint16 nPort);
    void close();
    int32 recvFrom(char *pBuf, uint32 nLen);
    int32 sendTo(const char * pBuf, uint32 nLen, const char *pIp, uint16 nPort);

private:
    ptrdiff_t           m_nSocket;

};

#endif

1.4 UdpSocket.cpp

#include "UdpSocket.h"
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cerrno>
#if defined I_OS_LINUX
#include <sys/socket.h>
#include <netinet/in.h>
#include <arpa/inet.h>
#include <unistd.h>
#elif defined I_OS_WIN32
#include <winsock2.h>
#pragma comment(lib,"ws2_32.lib")
#endif

/*********************** 
 Constrctor Function 
***********************/
UdpSocket::UdpSocket()
{
}

/***********************
 Destructor Function 
 **********************/
UdpSocket::~UdpSocket()
{
    close();
}

/**********************************************************
Function:         Connect(const char *pHost, uint16 nPort)
Description:      To connect appointed host Ip and Port
Return:           bool
Others:           None
**********************************************************/
bool UdpSocket::connect(const char* pHost, uint16 nPort)
{
    m_nSocket = socket(AF_INET, SOCK_DGRAM, IPPROTO_UDP);          // 使用UDP协议进行传输

    /* The receive buffer is set to 2MB */
    uint32 recv_buf = 8 * 1024 * 1024;
    setsockopt(m_nSocket, SOL_SOCKET, SO_RCVBUF, (const char*)&recv_buf, sizeof(uint32));    // 设置缓冲区的大小为

    /* The timeout is set to 3s */
#if defined I_OS_LINUX
    struct timeval TimeOut;
    TimeOut.tv_sec = 3;
    TimeOut.tv_usec = 0;
    setsockopt(m_nSocket, SOL_SOCKET, SO_RCVTIMEO, &TimeOut, sizeof(TimeOut));
#elif defined I_OS_WIN32
    uint32 TimeOut;
    TimeOut = 3;
    setsockopt(m_nSocket, SOL_SOCKET, SO_RCVTIMEO, (const char*)&TimeOut, sizeof(TimeOut));  // 设置接受时限为3ms
#endif

    sockaddr_in addr_sock;                           // 设置地址
    addr_sock.sin_family = AF_INET;                  // IPv4
    addr_sock.sin_port = htons(nPort);               // 绑定端口号
    addr_sock.sin_addr.s_addr = inet_addr(pHost);    // 绑定主机的IP地址

    return 0 == bind(m_nSocket, (sockaddr *)&addr_sock, sizeof(sockaddr_in));
}

/**********************************************************
Function:         Close
Description:      To close this socket
Return:           None
Others:           None
**********************************************************/
void UdpSocket::close()
{
    /* To shutdown for exit block */
    uint32 tmp = 1;
    
#if defined I_OS_LINUX
    shutdown(m_nSocket, SHUT_RDWR);
    setsockopt(m_nSocket, SOL_SOCKET, SO_REUSEADDR, &tmp, sizeof(uint32));
    ::close(m_nSocket);
#elif defined I_OS_WIN32
    shutdown(m_nSocket, SD_BOTH);
    setsockopt(m_nSocket, SOL_SOCKET, SO_REUSEADDR, (const char*)(&tmp), sizeof(uint32));   // 打开或关闭地址复用功能，这里是打开地址复用功能
    closesocket(m_nSocket);
#endif
}

/**********************************************************
Function:         RecvFrom(char * pBuf, uint32 nLen)
Description:      To receive data
Return:           int32
Others:           None
**********************************************************/
int32 UdpSocket::recvFrom(char * pBuf, uint32 nLen)
{
    sockaddr_in addr_sock;
    int32 len = sizeof(sockaddr);

#if defined I_OS_LINUX
    return recvfrom(m_nSocket, pBuf, nLen, 0, (sockaddr*)&addr_sock, reinterpret_cast<socklen_t *>(&len));
#elif defined I_OS_WIN32
    return recvfrom(m_nSocket, pBuf, nLen, 0, (sockaddr*)&addr_sock, &len);
#endif
}

/**********************************************************
Function:         Send(const char * pBuf, const uint32 nLen, const char *pIp, uint16 nPort)
Description:      To send data
Return:           int32
Others:           None
**********************************************************/
int32 UdpSocket::sendTo(const char * pBuf, const uint32 nLen, const char *pIp, uint16 nPort)
{
    sockaddr_in addr_sock;
    uint32 len = sizeof(sockaddr); 

    addr_sock.sin_family = AF_INET;
    addr_sock.sin_port = htons(nPort);
    addr_sock.sin_addr.s_addr = inet_addr(pIp);

    return sendto(m_nSocket, pBuf, nLen, 0, (sockaddr*)&addr_sock, len);
}

1.5 示例

#include "UdpSocket.h"

UdpSocket m_socket;

if (!m_socket.connect("192.168.11.124", 888)) {
	std::cout << "m_socket connect error." << std::endl;
}
else {
	std::cout << "m_socket is connected." << std::endl;
}

uint32 temp = 0;
if (m_socket.recvFrom((char*)temp, size)) {
	std::cout << "receive succeed." << std::endl;
} 
else {
	std::cout << "receive failed." << std::endl;
}

2. 文件读写

2.1 UFileReadSave.h

#ifndef __FILE_READ_SAVE_H
#define __FILE_READ_SAVE_H

#include <string>

class UFileReadSave
{
public:
	UFileReadSave(void) {};//inline
	~UFileReadSave(void) {};
public:
	/**********************************************************
	Description:		Get the bytes of file.
	Arguments:
	{
	strFilePath:	File path.
	fileByteNum:	The Bytes will be read.
	}
	Return:				Status.
	**********************************************************/
	static bool GetFileBytes(std::string& strFilePath, size_t& fileByteNum); // 20200810Updated
	/**********************************************************
	Description:		Read the binary file to an memory.
	Arguments:
	{
	pDst:			The destination memory.
	strFilePath:	File path.
	fileByteNum:	The Bytes will be read.
	}
	Return:				Status.
	**********************************************************/
	static bool ReadFile(char* pSrc, std::string& strFilePath, size_t fileByteNum);
	/**********************************************************
	Description:		Read the binary file to an memory with multi-processor.
	Arguments:
	{
	pDst:			The destination memory.
	strFilePath:	File path.
	fileByteNum:	The Bytes will be read.
	}
	Return:				Status.
	**********************************************************/
	static bool ReadFileMultiPrcs(void* pSrc, std::string& strFilePath, size_t fileByteNum, int threadNum = 0);
	/**********************************************************
	Description:		Write the binary file to the disk.
	Arguments:
	{
	pSrc:			The source data.
	strFilePath:	File path.
	fileByteNum:	The Bytes will be saved.
	}
	Return:				Status.
	**********************************************************/
	static bool SaveFile(char* pSrc, std::string& strFilePath, size_t fileByteNum);
};


#endif

2.2 UFileReadSave.cpp

#include "UFileReadSave.h"
#include <fstream>
#include <iostream>
#include <string>
#include <stdio.h>
#include <thread>
#include <vector>

using namespace std;


bool UFileReadSave::GetFileBytes(std::string& strFilePath, size_t& fileByteNum)
{
	ifstream projFile(strFilePath.c_str(), ios::binary);
	if (!projFile.is_open())
	{
		cout << strFilePath << " is losted!" << endl;        // file losted
		return false;
	}
	projFile.seekg(0, ios::end);
	fileByteNum = projFile.tellg();
	projFile.close();
	return true;
}

bool UFileReadSave::ReadFile(char* pDst, std::string& strFilePath, size_t fileByteNum)
{
	ifstream projFile(strFilePath.c_str(), ios::binary);
	if (!projFile.is_open())
	{
		cout << strFilePath << " is losted!" << endl;       // file losted
		return false;
	}
	projFile.seekg(0, ios::end);
	size_t fileLen = projFile.tellg();
	if (fileLen < fileByteNum)//file too small
	{
		cout << strFilePath << " only has " << fileLen << " bytes which is less than required " << fileByteNum << " bytes." << endl;
		return false;                                       // 20200908Updated
	}
	projFile.seekg(0, ios::beg);
	projFile.read(pDst, fileByteNum);
	projFile.close();
	return true;
}

bool UFileReadSave::SaveFile(char* pSrc, std::string& strFilePath, size_t fileByteNum)
{
	ofstream projFile(strFilePath.c_str(), ios::binary);
	if (!projFile.is_open())
	{
		cout << strFilePath << " open failed! Please select another path" << endl;//file losted
		return false;
	}
	projFile.write(pSrc, fileByteNum);
	projFile.close();
	return true;
}

2.3 使用方法

#include "UFileReadSave.h"

float* mich = new float[LORNum];
UFileReadSave::ReadFile((char*)mich, m_pPara->m_strMichPath, LORNum * sizeof(float))

2.3 多线程读取文件

• 多线程读取文件到一个数组

void test() {
	// 1. 读取文件大小，得到总的数据量
	for (unsigned i = 0; i < lChannelNum; ++i) {
		// 读取样本文件
		std::string lStrSampleFilePath = mCoinPetPara.m_strPETSavePath +
			"/thread" + std::to_string(mCoinPetPara.m_nThreadNum) +
			"-channel" + std::to_string(i);
		UFileReadSave::GetFileBytes(lStrSampleFilePath, lFileSize[i]);     // 给一个路径和一个存放位置即可
		//printf("FileSize[%d] = %f\n", i, double(lFileSize[i]));
		// 每个样本文件的数据量
		lSinglesNum[i] = lFileSize[i] / sizeof(DataFrameV2);
		printf("FileNum[%d] = %f\n", i, double(lSinglesNum[i]));
		lTotalSinglesNum += lSinglesNum[i];
		if (i >= 1) {
			lSingleReadBegin[i] += lSingleReadBegin[i - 1] + lSinglesNum[i - 1];
		}
	}

	DataFrameV2* lUdpFrame = new DataFrameV2[lTotalSinglesNum];
	
	std::vector<std::thread> lReadTasks;
	for (unsigned i = 0; i < lChannelNum; ++i) {
		lReadTasks.emplace_back([&, i] {
			// 读取样本文件
			std::string lStrSampleFilePath = mCoinPetPara.m_strPETSavePath +
			"/thread" + std::to_string(mCoinPetPara.m_nThreadNum) +
			"-channel" + std::to_string(i);
			UFileReadSave::ReadFile((char*)(lUdpFrame + lSingleReadBegin[i]), lStrSampleFilePath, lFileSize[i])
			});
	}
	for (auto& item : lReadTasks) {
		item.join();
	}


}

3. LM双指数拟合算法

3.1 LM拟合算法

#ifndef LMFIT_H
#define LMFIT_H

#include <math.h>
#include <string>
#include "maqr.h"
#include <iostream>
#include "lmstruct.h"

#define DERIV_STEP 1e-5
#define num_params 4  //参数个数
#define total_data 8
#define MAX_ITER 500
#define MIN_SETP 1e-12
#define GTOL 1e-12

double Func(double* x, double* par, int n) {

	return par[0] * exp(par[1] * x[n]) + par[2] * exp(par[3] * x[n]);
}

double Deriv(double(*Func)(double* x, double* par, int n), double* x, double* par, int j, int k) {
	double par1[num_params];
	double par2[num_params];
	for (unsigned i = 0; i < num_params; ++i) {
		par1[i] = par2[i] = par[i];
	}
	par1[k] -= DERIV_STEP;
	par2[k] += DERIV_STEP;

	double p1 = Func(x, par1, j);
	double p2 = Func(x, par2, j);

	return (p2 - p1) / (2 * DERIV_STEP);

}

double calBottom(double* hlm, float u, double* c) {
	double b[num_params] = { 0 };

	double sum = 0;
	for (unsigned i = 0; i < num_params; ++i) {
		b[i] = u * hlm[i];
		sum += hlm[i] * (b[i] - c[i]);
	}

	return sum;
}

//求取方程左边方程
void CalLeft(double* Jf, float u, double* a) {
	for (unsigned i = 0; i < num_params; ++i) {
		for (unsigned j = 0; j < num_params; ++j) {
			for (unsigned k = 0; k < total_data; ++k) {
				a[i * num_params + j] += Jf[k * num_params + i] * Jf[k * num_params + j];
			}
			a[i * num_params + j] += u;
		}
	}
}
void CalRight(double* q, double* Jf, double* r, double* c, double* d) {
	//double c[num_params] = { 0 };

	for (unsigned i = 0; i < num_params; ++i) {
		for (unsigned k = 0; k < total_data; ++k) {
			c[i] += Jf[k * num_params + i] * r[k];
		}
	}

	for (unsigned i = 0; i < num_params; ++i) {
		for (unsigned j = 0; j < num_params; ++j) {
			d[i] += q[j * num_params + i] * c[j];
		}
	}
}
double norm(double* r, unsigned m) {
	double sum = 0;
	for (unsigned i = 0; i < m; ++i) {
		sum += r[i] * r[i];
	}
	return sqrt(sum);
}

/*利用L - M方法求解 原方程（J^(T)*J + u * I） * hlm = J ^ (T)*r
	对（J^(T)*J + u * I）做QR分解，并作运算变换得到
	R * hlm = Q^(T) * J^(T) * r  根据此方程求解hlm
	利用实际下降量和预测下降量的比值来判断拟合效果
*/
void LM(double(*Func)(double* x, double* par, int n), double* x, double* y, double* par, lm_status_struct* s) {
	int m = 8; //数据个数
	double* r = new double[m];   //存放阈值与拟合值的误差
	double* r_temp = new double[m];
	double* par_temp = new double[num_params];
	double* Jf = new double[m * num_params];  // 雅各比矩阵

	double last_mse = 0;

	float u = 1, v = 2;
	double mse = 0;

	s->outcome = 1;              // 迭代100次退出

	int i = 0;
	for (i = 0; i < MAX_ITER; ++i) {
		mse = 0;
		double mse_temp = 0;
		//计算雅各比行列式 同时得到误差
		for (unsigned j = 0; j < m; ++j) {
			//计算y与拟合的值之间的误差
			r[j] = y[j] - Func(x, par, j);
			mse += r[j] * r[j];

			for (unsigned k = 0; k < num_params; ++k) {
				Jf[j * num_params + k] = Deriv(Func, x, par, j, k);
			}
		}
		mse /= m;  //求取平均误差
		//std::cout << "mse: " << mse << std::endl;

		//for (unsigned j = 0; j < total_data; ++j) {
		//	for (unsigned k = 0; k < num_params; ++k) {
		//		std::cout << Jf[j * num_params + k] << " ";
		//	}
		//	std::cout << std::endl;
		//}
		//std::cout << std::endl;


		double hlm[num_params] = { 0 };
		double a[num_params * num_params] = { 0 };  // 用来存放（J^(T) * J + u * I）的结果
		double q[num_params * num_params] = { 0 };  //a矩阵QR分解后的Q矩阵
		double d[num_params] = { 0 };

		CalLeft(Jf, u, a);   //计算转换方程的等式左边的值 结果存放在a中


		//QR分解  分解完以后a = R， q = Q
		maqr(a, num_params, num_params, q, s);

		//计算右边Q^(T) * J^(T) * r
		double c[num_params] = { 0 };
		CalRight(q, Jf, r, c, d);  //计算方程右边的值，用d来存放

		//检测退化
		if (norm(c, num_params) < GTOL) {
			//std::cout << "norm(c, num_params) < GTOL" << std::endl;
			s->outcome = 3;           // 矩阵退化退出
			break;
		}

		for (int i = num_params - 1; i >= 0; --i) {
			for (int j = i + 1; j < num_params; ++j) {
				d[i] -= a[i * num_params + j] * hlm[j];
			}
			hlm[i] = d[i] / a[i * num_params + i];  //求解步长
			par_temp[i] = par[i] + hlm[i];  //更新参数
		}

		for (unsigned j = 0; j < total_data; ++j) {

			r_temp[j] = y[j] - Func(x, par_temp, j);
			mse_temp += r_temp[j] * r_temp[j];  //计算更新参数后的总误差
		}
		mse_temp /= total_data;

		//std::cout << "mse_temp: " << mse_temp << std::endl;

		double bottom = calBottom(hlm, u, c);   //计算判断拟合效果式子的分母
		double qqq = (mse - mse_temp) / (0.5 * bottom);        // 这里是一个矩阵运算，判断拟合效果的依据


		//std::cout << "qqq: " << qqq << std::endl;

		if (qqq > 0) {

			double s = 1.0 / 3.0;
			v = 2;

			mse = mse_temp;            // 更新误差
			for (unsigned j = 0; j < num_params; ++j) {
				par[j] = par_temp[j];
			}
			double temp = 1 - pow(2 * qqq - 1, 3);

			if (s > temp)
			{
				u = u * s;
			}
			else
			{
				u = u * temp;
			}
		}
		else {
			u = u * v;
			v = 2 * v;
		}


		double setp = norm(hlm, num_params);


		// The difference in mse is very small, so quit
		// 退出条件：步长收敛
		if (setp <= MIN_SETP) {
			//std::cout << "setp <= MIN_SETP" << std::endl;
			s->outcome = 4;          // 步长很小退出
			break;
		}
		//退出条件:两次更新后的误差差值很小，已经收敛
		if (qqq > 0 && fabs(mse - last_mse) < 1e-3) {
			//std::cout << "qqq > 0 && fabs(mse - last_mse) < 1e-3" << std::endl;
			s->outcome = 5;          // 两次误差很小退出
			break;
		}

		//printf("%d %lf\n", i, mse);

		last_mse = mse;
	}
	delete[] r;
	delete[] r_temp;
	delete[] par_temp;
	delete[] Jf;

	s->fnorm = sqrt(mse * m);
	s->nfev = i;


}

//int main() {
//	//a*exp^(b*x) + c*exp^(d*x)
//
//	int totoal_data = 8;
//
//	//double x[8] = { 0,0.740742,1.2037,1.85185,30.2778,45.3704,59.6296,80.9259 };
//	//double x[8] = { 0,1.18182,1.54546,2.45455,33.5455,48.0909,57.3636,84.3636 };
//	//double x[8] = { 0,0.545454,0.81818,1.36364,40.0909,47.8182,63.5455, 87.1818 };
//	//double x[8] = { 0,0.545456,0.818182,1.18182,49.2727,60.6364,77.1818,99.4545 };
//	//double x[8] = { 0,0.545454,0.454544,1.72727,49.7273,62,76.8182,110.455 };
//	double x[8] = { 0,0.545454,1,1.54545,34.8182,45.277,57.8182,86.9091 };
//
//	double y[8] = { 20,90,160,230,230,160,90,20 };
//	double par[4] = { 1873, -0.0448, -1847, -0.1429 };
//
//	LM(Func, x, y, par);
//
//	std::cout << "拟合值: ";
//	for (int i = 0; i < 8; ++i) {
//		std::cout << Func(x, par, i) << "   ";
//	}
//	std::cout << std::endl;
//}
//


#endif // !LMFIT_H

3.2 QR分解

#pragma once
#ifndef MAQR_H
#define MAQR_H


#include "math.h"
#include "stdio.h"

__device__ int maqr(double* a, unsigned m, unsigned n, double* q, lm_status_struct& s)
{
	int i, j, k, l, nn, p, jj;
	double u, alpha, w, t;
	if (m < n)
	{
		printf("fail\n"); return(0);
	}
	for (i = 0; i <= m - 1; i++)
		for (j = 0; j <= m - 1; j++)
		{
			l = i * m + j; q[l] = 0.0;
			if (i == j) q[l] = 1.0;
		}
	nn = n;
	if (m == n) nn = m - 1;
	for (k = 0; k <= nn - 1; k++)
	{
		u = 0.0; l = k * n + k;
		for (i = k; i <= m - 1; i++)
		{
			w = fabs(a[i * n + k]);
			if (w > u) u = w;
		}
		alpha = 0.0;
		for (i = k; i <= m - 1; i++)
		{
			t = a[i * n + k] / u; alpha = alpha + t * t;
		}
		if (a[l] > 0.0) u = -u;
		alpha = u * sqrt(alpha);
		if (fabs(alpha) + 1.0 == 1.0)
		{
			//printf("fail\n"); 
			s.outcome = 2;         // qr分解失败退出
			return(0);
		}
		u = sqrt(2.0 * alpha * (alpha - a[l]));
		if ((u + 1.0) != 1.0)
		{
			a[l] = (a[l] - alpha) / u;
			for (i = k + 1; i <= m - 1; i++)
			{
				p = i * n + k; a[p] = a[p] / u;
			}
			for (j = 0; j <= m - 1; j++)
			{
				t = 0.0;
				for (jj = k; jj <= m - 1; jj++)
					t = t + a[jj * n + k] * q[jj * m + j];
				for (i = k; i <= m - 1; i++)
				{
					p = i * m + j; q[p] = q[p] - 2.0 * t * a[i * n + k];
				}
			}
			for (j = k + 1; j <= n - 1; j++)
			{
				t = 0.0;
				for (jj = k; jj <= m - 1; jj++)
					t = t + a[jj * n + k] * a[jj * n + j];
				for (i = k; i <= m - 1; i++)
				{
					p = i * n + j; a[p] = a[p] - 2.0 * t * a[i * n + k];
				}
			}
			a[l] = alpha;
			for (i = k + 1; i <= m - 1; i++)
				a[i * n + k] = 0.0;
		}
	}
	for (i = 0; i <= m - 2; i++)
		for (j = i + 1; j <= m - 1; j++)
		{
			p = i * m + j; l = j * m + i;
			t = q[p]; q[p] = q[l]; q[l] = t;
		}
	return(1);
}


#endif // !MAQR_H

4. Cuda编程技巧

// 检错函数
void checkCudaError(cudaError_t error， const string& errmsg)
{
    if (error != cudaSuccess)
    {
    	std::cout << "error position: " << errmsg << std::endl;
        printf("Error in CUDA function.\nError: %s\n", cudaGetErrorString(error));
        getchar();
        exit(EXIT_FAILURE);
    }
}


// 分配GFU内存
cudaError_t error;
error = cudaMalloc((void **)&_d_keys, arrayLength * sizeof(*_d_keys));
checkCudaError(error, "_d_keys malloc");

// 进行数据拷贝
error = cudaMemcpy((void *)_d_keys, h_keys, arrayLength * sizeof(*h_keys), cudaMemcpyHostToDevice);
checkCudaError(error, "_d_keys memcpy");

// 数据初始化
error = cudaMemset(_d_keys, 0, arrayLength * sizeof(_d_keys));
checkCudaError(error, "_d_keys memset")

// GPU并行操作
cudaDeviceSynchronize();
error = cudaGetLastError();
checkCudaError(error, "GaussionConv3D_kernel")

// 进行数据拷贝
error = cudaMemcpy(h_keys, (void *)_d_keys, _arrayLength * sizeof(*_h_keys), cudaMemcpyDeviceToHost);
checkCudaError(error, "h_keys memcpy");

// 内存清除
error = cudaFree(_d_keys);
checkCudaError(error, "_d_keys free");

5. 类型定义和指针释放

#ifndef __TYPEDEF_H
#define __TYPEDEF_H

typedef unsigned char             uint8;
typedef unsigned short            uint16;
typedef unsigned int              uint32;
typedef unsigned long long int    uint64;
typedef signed char               int8;
typedef short                     int16;
typedef int                       int32;
typedef signed long long int      int64;


#define RELEASE_POINTER(p) if(p) {delete p; p = nullptr;}
#define RELEASE_ARRAY_POINTER(p) if(p) {delete [] p; p = nullptr;}
#define OFFSET_ADDR(STRUCTURE,FIELD) ((size_t)(&((STRUCTURE*)nullptr)->FIELD))

#endif

6. 多线程调用

void AcquirePET::start() {
	mStopFlag = false;
	std::vector<std::thread> lAcqTaskVec;
	for (uint32 i = 0; i < mBdmNum; ++i) { /// Acquire Receive Tasks
		lAcqTaskVec.emplace_back(std::thread(acqRecvTask, this, i));
	}
	/// Acquire Stop Task
	lAcqTaskVec.emplace_back(std::thread(acqStopTask, this));       // 注意这里只是放入了一个线程
	/// Start Record Task
	lAcqTaskVec.emplace_back(std::thread(acqStartRecordTask, this));       // 注意这里只是放入了一个线程
	for (auto& task : lAcqTaskVec) {
		task.join();
	}
	for (uint32 i = 0; i < mBdmNum; ++i) {
		printf("[Acquire Receive Task %d] had received %d packages in total.\n", i, mRecvPkgNumVec[i]);
	}
}


void AcquirePET::acqRecvTask(AcquirePET* self, uint32 tid) {
	// std::unique_lock<std::mutex> lLock{self->mMtx};
	// Ref from parameters. 
	std::vector<UdpSocket>& lSockets{ self->mUdpSocketVec };               // 12个线程的socket
	DataBuffer<DataFrameV2>* lBuffer{ &self->mDataBufferVec[tid] };        // 该线程对应的缓冲区
	DataBuffer<DataFrameV2>* lBufferBackup{ &self->mDataBufferBackupVec[tid] };         // 该线程对应的辅助缓冲区
	DataBuffer<DataFrameV2>* lBufferTemp;            // 临时缓冲区

	uint32& lTaskIndex{ tid };                         // 线程id
	uint32& lFrameNumOnePkg{ self->mFrameNumOnePkg };           // 一个 UDP 包内的事件帧个数
	uint32& lPkgNum{ self->mPkgNum };                  // 缓冲区可以容纳的udp包的个数
	uint32& lPkgSize{ self->mPkgSize };                // udp包的大小
	std::atomic<uint32>& lSwitchNum{ self->mSwitchNum };  // 使用初始化列表而不使用赋值运算符
	std::atomic<uint32>& lFinishNum{ self->mFinishNum };

	bool& lStopFlag{ self->mStopFlag };
	std::vector<bool>& lSwitches{ self->mIsSwitchableVec };

	/// Local Variables。
	uint32 lSegRecvPkgNum{ 0 };
	uint32 lTtlRecvPkgNum{ 0 };

	char lPkgTemp[lPkgSize];    /* Process one package once */

	printf("[Acquire Receive Task %d] starts.\n", lTaskIndex);

	// 这里是用12个线程去采集12个探测单元的数据
	while (!lStopFlag) {
		if (lSockets[lTaskIndex].recvFrom((char*)lPkgTemp, lPkgSize) == lPkgSize) {
			++lSegRecvPkgNum;
			++lTtlRecvPkgNum;
			++self->mTtlRecvPkgNum;                  // 总接收的UDP包的个数
			lBuffer->writeData((char*)lPkgTemp, lFrameNumOnePkg);
			/// 如果缓冲区处于可交换状态并且已经有其他缓冲区发生了交换（一个缓冲区发生交换，所有缓冲区都需要随之发生交换）
			/// 或者当前缓冲区已满（第一个发生交换的缓冲区）
			if ((lSwitchNum > 0 && lSwitches[lTaskIndex]) || lSegRecvPkgNum >= lPkgNum) {
				uint32 lErrorSwitch{ 0 };
				// 辅助缓存区如果不为空说明，还没有完全写入文件中
				while (lBufferBackup->getCurFrameNum() != 0) {   /// Confirm the backup one is cleared.
					if (lStopFlag) break;                        // 跳出最近一次的循环
					++lErrorSwitch;
				}
				if (lErrorSwitch > 0) {
					printf("[Acquire Receive Task %d] failed switching buffer for %d times!\n", lTaskIndex, lErrorSwitch);
				}
				printf("[Acquire Receive Task %d] has %d frames, now switching buffer...\n", lTaskIndex, lBuffer->getCurFrameNum());
				++lSwitchNum;                          // 交换次数+1
				lSegRecvPkgNum = 0;

				/// Swap buffers and bufferBackups...
				lBufferTemp = lBuffer;
				lBuffer = lBufferBackup;
				lBufferBackup = lBufferTemp;

				lSwitches[lTaskIndex] = false;
			}
		}
	}
	printf("[Acquire Receive Task %d] ended with %d frames.\n", lTaskIndex, lBuffer->getCurFrameNum());
	self->mRecvPkgNumVec[lTaskIndex] = lTtlRecvPkgNum;
	++lFinishNum;        // 当12个线程的采集任务都结束了，也需要停止
}

7. 类的构造

7.1 AcquirePET.h

class AcquirePET {
public:
    explicit AcquirePET(AcquirePetPara&);
    virtual ~AcquirePET();
    
    void start();
private:
    void init();    /// 初始化
    static void acqRecvTask(AcquirePET*, uint32);   /// 采集任务，有 mBdmNum 个
    static void acqStopTask(AcquirePET*);   /// 停止与监控任务
    
    std::mutex mMtx;    /// 互斥量
    std::atomic<uint32> mActiveRecordNum{0};    /// 记录任务中活动记录个数
    std::atomic<uint32> mTtlRecordNum{0};   /// 记录任务中总任务个数


    static void acqStartRecordTask(AcquirePET*); /// 记录任务的启动管理
    static void acqRecordTask(AcquirePET*);   /// 记录任务
    
    AcquirePetPara mAcquirePetPara;
    
    uint32 mBdmNum; /// BDM（基本探测单元）个数
    uint32 mFrameNumOnePkg; /// 一个 UDP 包内的事件帧个数
    uint32 mPkgSize;    /// UDP 包的大小
    uint32 mPkgNum; /// 缓冲区组中一个缓冲区中可容纳 UDP 包的个数
    
    std::string mUdpIp; /// UDP IP 地址
    std::vector<unsigned int> mUdpPortVec;  /// UDP 端口号
    
    std::string mClockIp;   /// 时钟板 IP 地址
    uint32 mClockPort;  /// 时钟板端口号
    
    std::string mPetSavePath;   /// PET 采集数据存储地址（文件夹）
    bool mOnlineCoin;   /// 是否进行在线符合（已弃用）
    uint32 mPetScanTime;    /// PET 采集时间（单位：s）
    
    std::vector<UdpSocket> mUdpSocketVec; /// UDP Socket 数组
    DataBuffer<DataFrameV2>* mDataBufferVec;    /// 主缓冲区组
    
    double mProgress{0.0}; /// 采集进度，0.0 ~ 1.0 代表正在进行采集，2.0 代表结束
    double mKcps{0.0};  /// 千事件 / 秒，衡量采集速率的一种标准
    
    bool mIsBackupDataBufferActive{false};  /// 主从缓冲区交换的标志位
    DataBuffer<DataFrameV2>* mDataBufferBackupVec; /// 从/备用 缓冲区组
    
    bool mStopFlag{false};  /// 采集监听任务结束的标志位（但此时其他任务仍然存在！）
    bool mRecordFlag{false};    /// 发起缓冲区记录的标志位（结束监听任务 -> 开始记录任务）
    bool mAcquireFinishFlag{false}; /// 采集真正结束的标志位
    
    std::atomic<uint32> mSwitchNum{0};  /// 采集任务发生缓冲区切换者个数
    std::atomic<uint32> mFinishNum{0};  /// 采集任务结束个数
    std::vector<bool> mIsSwitchableVec; /// 是否能进行缓冲区切换
    
    std::atomic<uint32> mTtlRecvPkgNum{0};  /// 总接收 UDP 包个数
    std::vector<uint32> mRecvPkgNumVec; /// 每个采集任务接收 UDP 包个数
};

7.2 AcquirePET.cpp

AcquirePET::AcquirePET(AcquirePetPara& app) {
	mAcquirePetPara = app;

	mBdmNum = app.m_nBDM_Num;
	mFrameNumOnePkg = app.m_nFrameNumOnePkg;
	mPkgSize = app.m_nPkgSize;
	mPkgNum = app.m_nPkgNum;

	mUdpIp = app.m_strUdpIp;
	mUdpPortVec = app.m_pUdpPort;
	mClockIp = app.m_strClockIp;
	mClockPort = app.m_nClockPort;

	mPetSavePath = app.m_strPETSavePath;
	mOnlineCoin = app.m_bOnlineCoin;

	mPetScanTime = app.m_nPETScanTime;

	mDataBufferVec = new DataBuffer<DataFrameV2>[mBdmNum];       // 为每一个BdmId创建一个缓冲区
	mDataBufferBackupVec = new DataBuffer<DataFrameV2>[mBdmNum];

	mUdpSocketVec.resize(mBdmNum);
	mIsSwitchableVec.reserve(mBdmNum);                           // 判断是否要进行切换，缓冲区满则进行切换
	mRecvPkgNumVec.reserve(mBdmNum);
	for (uint32 i = 0; i < mBdmNum; ++i) {
		mIsSwitchableVec.push_back(false);
		mRecvPkgNumVec.push_back(0);
	}
	init();
}

AcquirePET::~AcquirePET() {
	RELEASE_ARRAY_POINTER(mDataBufferVec);
	RELEASE_ARRAY_POINTER(mDataBufferBackupVec);
}

8. 内存初始化

memset(mDataBuffer, 0, mBufferSize);