原代碼
數學思路
point類
//檔案1，類的定義，point.h
#ifndef _POINT_H
#define _POINT_H
class Point{
public:
Point(float x=0, float y=0) :x(x), y(y) {}
float GetX() const {return x;}
float GetY() const {return y;}
private:
float x,y;
};
#endif
主函數
//主函數，main.cpp
#include "point.h"
#include <iostream>
#include <cmath>
using namespace std;
//直線線性拟合，points 為各點，nPoint 為點數
float lineFit(const Point points[], int nPoint) //友元函數體
{
float avgX=0,avgY=0; //定義變量
float lxx=0,lyy=0,lxy=0;
for(int i=0;i<nPoint;i++) //計算X、Y的平均值
{
avgX+=points[i].GetX()/nPoint;
avgY+=points[i].GetY()/nPoint;
}
for(int i=0;i<nPoint;i++) //計算Lxx、Lyy和Lxy
{
lxx+=(points[i].GetX()-avgX)*(points[i].GetX()-avgX);
lyy+=(points[i].GetY()-avgY)*(points[i].GetY()-avgX);
lxy+=(points[i].GetX()-avgX)*(points[i].GetY()-avgY);
}
cout<<"This line can be fitted by y=ax+b."<<endl;
cout<<"a="<<lxy/lxx; //輸出回歸系數a
cout<<" b="<<avgY-lxy*avgX/lxx<<endl; //輸出回歸系數b
return float(lxy/sqrt(lxx*lyy)); //傳回相關系數r
}
int main()
{
Point p[10]={Point(6,10),Point(14,20),Point(26,30),
Point(33,40),Point(46,50),Point(54,60),Point(67,70),
Point(75,80),Point(84,90),Point(100,100)}; //初始化資料點
float r=lineFit(p,10); //進行線性回歸計算
cout<<"Line coefficient r="<<r<<endl; //輸出相關系數
return 0;
}
輸出結果
This line can be fitted by y=ax+b.
a=0.97223 b=5.90237
Line coefficient r=0.998193
Program ended with exit code: 0
`
如上的部分轉載自https://www.jianshu.com/p/128924b2c846?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation
下面的部分為自己修改的代碼
修改為avx+openmp多線程的代碼
```cpp
#include <iostream>
#include <vector>
#include <omp.h>
#include "math.h"
#include <cstdlib>
#include <immintrin.h>
using namespace std;
class Point{
public:
Point(double x=0, double y=0) : x(x), y(y) {}
float GetX() const {return x;}
float GetY() const {return y;}
private:
float x,y;
};
//type為類型，size為p的大小
double add_help(__m256d num, int type, int size, vector<Point> &p, double avgx, double avgy) {
double d[4];
_mm256_storeu_pd(d, num);
double res = 0;
for (int i = 0; i < 4; i++) {
res += d[i];
}
if (type == 0) {
for (int i = size / 4 * 4; i < size; i++) {
res += p[i].GetX();
}
} else if (type == 1) {
for (int i = size / 4 * 4; i < size; i++) {
res += p[i].GetY();
}
} else if (type == 2) {
for (int i = size / 4 * 4; i < size; i++) {
res += (p[i].GetX() - avgx) * (p[i].GetX() - avgy);
}
return res;
}else if( type == 3){
for (int i = size / 4 * 4; i < size; i++) {
res += (p[i].GetY()-avgy)*(p[i].GetY()-avgx);
}
return res;
}else {
for (int i = size / 4 * 4; i < size; i++) {
res += (p[i].GetX()-avgx)*(p[i].GetY()-avgy);
}
return res;
}
return res / size;
}
int main() {
// vector<Point> p(
// {Point(6, 10), Point(14, 20), Point(26, 30), Point(33, 40), Point(46, 50), Point(54, 60), Point(67, 70),
// Point(75, 80), Point(84, 90), Point(100, 100)});
vector<Point> p;
//y = ax + b;
int aa = 15, bb = 45;
//設定p的size大小，産生size對x，y的值
int size = 300;
for(int i = 0; i < size; i++){
//0-99的随機數
double a = rand()%100 * 1.0;
//0-1的浮點數,為噪聲
double random_01 = rand() / double(RAND_MAX);
p.push_back({a,a*aa + bb + random_01});
}
double avgx = 0, avgy = 0;
//double為64位，是以256裝4個double
//初始化sumx為__m256d類型的4個全0
__m256d sumx = _mm256_set1_pd(0.0);
__m256d sumy = _mm256_set1_pd(0.0);
//8線程并行求和
#pragma omp parallel for num_threads(8)
for (int i = 0; i < size / 4 * 4; i += 4) //計算X、Y的和，先算出4的倍數個數的和
{
//裝入4個double數
__m256d a = _mm256_set_pd(p[i + 3].GetX(),
p[i + 2].GetX(),
p[i + 1].GetX(),
p[i + 0].GetX());
//兩個__m256d相加
sumx = _mm256_add_pd(a, sumx);
__m256d b = _mm256_set_pd(p[i + 3].GetY(),
p[i + 2].GetY(),
p[i + 1].GetY(),
p[i + 0].GetY());
sumy = _mm256_add_pd(b, sumy);
}
//補上餘下的4個以内的數
avgx = add_help(sumx, 0, size, p, avgx, avgy);
avgy = add_help(sumy, 1, size, p, avgx, avgy);
cout << "avgx = " << avgx << endl;
cout << "avgy = " << avgy << endl;
//avgx_256設定為4個avgx
__m256d avgx_256 = _mm256_set1_pd(avgx);
__m256d avgy_256 = _mm256_set1_pd(avgy);
double lxx = 0, lyy = 0, lxy = 0;
__m256d lxx_256 = _mm256_set1_pd(0.0);
__m256d lyy_256 = _mm256_set1_pd(0.0);
__m256d lxy_256 = _mm256_set1_pd(0.0);
#pragma omp parallel for num_threads(8)
for (int i = 0; i < size / 4 * 4; i += 4) //計算Lxx、Lyy和Lxy
{
__m256d a = _mm256_set_pd(p[i + 3].GetX(),
p[i + 2].GetX(),
p[i + 1].GetX(),
p[i + 0].GetX());
__m256d b = _mm256_set_pd(p[i + 3].GetY(),
p[i + 2].GetY(),
p[i + 1].GetY(),
p[i + 0].GetY());
__m256d t = _mm256_sub_pd(a, avgx_256);
__m256d f = _mm256_sub_pd(a, avgy_256);
//先相乘再相加
lxx_256 = _mm256_add_pd(_mm256_mul_pd(t, f), lxx_256);
t = _mm256_sub_pd(b, avgy_256);
f = _mm256_sub_pd(b, avgx_256);
lyy_256 = _mm256_add_pd(_mm256_mul_pd(t, f), lyy_256);
t = _mm256_sub_pd(a, avgx_256);
f = _mm256_sub_pd(b, avgy_256);
lxy_256 = _mm256_add_pd(_mm256_mul_pd(t, f), lxy_256);
}
lxx = add_help(lxx_256, 2, size, p, avgx, avgy);
lxy = add_help(lxy_256, 4, size, p, avgx, avgy);
lyy = add_help(lyy_256, 3, size, p, avgx, avgy);
cout<<"a="<<lxy/lxx; //輸出回歸系數a
cout<<" b="<<avgy-lxy*avgx/lxx<<endl; //輸出回歸系數b
cout<<"Line coefficient r="<<float(lxy/sqrt(lxx*lyy))<<endl; //輸出相關系數
return 0;
}
![swmm與lisflood-fp源碼如何一起編譯 CMake指令[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)