sunwen/Aurora
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add prod, dot and their unit test.

Browse Source
This commit is contained in:
Krad
2023-04-25 17:28:33 +08:00
parent de1c1b0dda
commit cf96572074
3 changed files with 179 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

139
src/Function2D.cpp
View File

@@ -288,37 +288,67 @@ Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& row
}
Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
std::cerr
<< (aMatrix.getDimSize(2) > 1 ? "sum() not support 3D data!" : "sum() not support complex value type!")
if (aMatrix.getDimSize(2)>1 ) {
std::cerr<< "sum() not support 3D data!"
<< std::endl;
return Matrix();
}
switch (direction)
{
case All:
if (aMatrix.isComplex()){
switch (direction)
{
Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
double * ret = malloc(1);
ret[0] = srcV.array().sum();
return Matrix::New(ret,1);
case All:
{
Eigen::Map<Eigen::VectorXcd> srcV((std::complex<double>*)aMatrix.getData(),aMatrix.getDataSize());
std::complex<double>* ret = (std::complex<double>*)malloc(1,true);
ret[0] = srcV.array().sum();
return Matrix::New((double*)ret,1,1,1,Complex);
}
case Row:
{
Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
std::complex<double> * ret = (std::complex<double>*)malloc(aMatrix.getDimSize(0),true);
Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(0));
retV = srcM.rowwise().sum();
return Matrix::New((double*)ret,aMatrix.getDimSize(0),1,1,Complex);
}
case Column:
default:
{
Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
std::complex<double>* ret = (std::complex<double>*)malloc(aMatrix.getDimSize(1),true);
Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(1));
retV = srcM.colwise().sum();
return Matrix::New((double*)ret,1,aMatrix.getDimSize(1),1,Complex);
}
}
case Row:
}
else{
switch (direction)
{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
double * ret = malloc(aMatrix.getDimSize(0));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
retV = srcM.rowwise().sum();
return Matrix::New(ret,aMatrix.getDimSize(0),1);
}
case Column:
default:
{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
double * ret = malloc(aMatrix.getDimSize(1));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
retV = srcM.colwise().sum();
return Matrix::New(ret,1,aMatrix.getDimSize(1));
case All:
{
Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
double * ret = malloc(1);
ret[0] = srcV.array().sum();
return Matrix::New(ret,1);
}
case Row:
{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
double * ret = malloc(aMatrix.getDimSize(0));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
retV = srcM.rowwise().sum();
return Matrix::New(ret,aMatrix.getDimSize(0),1);
}
case Column:
default:
{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
double * ret = malloc(aMatrix.getDimSize(1));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
retV = srcM.colwise().sum();
return Matrix::New(ret,1,aMatrix.getDimSize(1));
}
}
}
}
@@ -624,7 +654,7 @@ Matrix Aurora::ifft(const Matrix &aMatrix) {
Matrix Aurora::hilbert(const Matrix &aMatrix) {
auto x = fft(aMatrix);
auto h = new double[aMatrix.getDimSize(0)];
auto h = malloc(aMatrix.getDimSize(0));
auto two = 2.0;
auto zero = 0.0;
cblas_dcopy(aMatrix.getDimSize(0), &zero, 0, h, 1);
@@ -641,3 +671,60 @@ Matrix Aurora::hilbert(const Matrix &aMatrix) {
return result;
}
Matrix Aurora::prod(const Matrix &aMatrix) {
if (aMatrix.getDimSize(2) > 1 ) {
std::cerr<< "prod() not support 3D data!"
<< std::endl;
return Matrix();
}
if (aMatrix.isComplex()){
Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
auto ret = malloc(aMatrix.getDimSize(1),true);
Eigen::Map<Eigen::VectorXcd> retV((std::complex<double>*)ret,aMatrix.getDimSize(1));
retV = srcM.colwise().prod();
return Matrix::New(ret,1,aMatrix.getDimSize(1),1,Complex);
}
else{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
auto ret = malloc(aMatrix.getDimSize(1));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
retV = srcM.colwise().prod();
return Matrix::New(ret,1,aMatrix.getDimSize(1));
}
}
Matrix Aurora::dot(const Matrix &aMatrix,const Matrix& aOther,FunctionDirection direction ) {
if ( direction == All){
std::cerr<< "dot() not support 3D data!"
<< std::endl;
return Matrix();
}
if (!aMatrix.compareShape(aOther)){
std::cerr<< "dot() matrix must be same shape!"
<< std::endl;
return Matrix();
}
if (aMatrix.isComplex()){
return sum(conj(aMatrix)*aOther,direction);
}
else{
if (direction == Column)
{
auto ret = malloc(aMatrix.getDimSize(1));
for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
ret[i]=cblas_ddot(aMatrix.getDimSize(0),aMatrix.getData()+i*aMatrix.getDimSize(0),1,
aOther.getData()+i*aMatrix.getDimSize(0),1);
}
return Matrix::New(ret,1,aMatrix.getDimSize(1),1);
}
else{
auto ret = malloc(aMatrix.getDimSize(0));
for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
ret[i] = cblas_ddot(aMatrix.getDimSize(1),aMatrix.getData()+i,aMatrix.getDimSize(0),
aOther.getData()+i,aMatrix.getDimSize(0));
}
return Matrix::New(ret,aMatrix.getDimSize(1),1,1);
}
}
}

11
src/Function2D.h
View File

@@ -48,7 +48,7 @@ namespace Aurora {
Matrix min(const Matrix& aMatrix,const Matrix& aOther);
/**
* 求矩阵和，可按行、列、单元, 目前不支持三维，不支持复数
* 求矩阵和，可按行、列、单元, 目前不支持三维
* @param aMatrix 目标矩阵
* @param direction 方向，Column, Row, All
* @return 求和结果矩阵
@@ -120,6 +120,15 @@ namespace Aurora {
* @return
*/
Matrix hilbert(const Matrix& aMatrix);
/**
* prod，支持到2维
* @param aMatrix
* @return
*/
Matrix prod(const Matrix& aMatrix);
Matrix dot(const Matrix& aMatrix,const Matrix& aOther,FunctionDirection direction = Column);
};

56
test/Function2D_Test.cpp
View File

@@ -342,6 +342,62 @@ TEST_F(Function2D_Test, hilbert) {
EXPECT_DOUBLE_EQ(fourDecimalRound(result[11].imag()),0.3249);
}
TEST_F(Function2D_Test, prod) {
double *dataB = new double[20]{1.1, 2.6, 3.8, 6.2,
4.3, 5.7, 6.9, 10.6,
7.1, 8.3, 9.7, 11.2,
17.8, 13.3,26.5 , -7.7,
9.9, 8.2, 6.3, 5.1};
auto B = Aurora::Matrix::fromRawData(dataB, 4, 5);
auto ret = Aurora::prod(B);
EXPECT_EQ(1, ret.getDimSize(0));
EXPECT_EQ(5, ret.getDimSize(1));
EXPECT_DOUBLE_EQ(0.0067, fourDecimalRound(ret.getData()[0]/10000));
EXPECT_DOUBLE_EQ(0.1793, fourDecimalRound(ret.getData()[1]/10000));
EXPECT_DOUBLE_EQ(0.6402, fourDecimalRound(ret.getData()[2]/10000));
EXPECT_DOUBLE_EQ(-4.8307, fourDecimalRound(ret.getData()[3]/10000));
EXPECT_DOUBLE_EQ(0.2608, fourDecimalRound(ret.getData()[4]/10000));
}
TEST_F(Function2D_Test, dot) {
{
double *dataA = new double[8];
double *dataB = new double[8];
for (int i = 0; i < 8; ++i) {
dataA[i] = (double) (i);
dataB[i] = (double) (2);
}
Aurora::Matrix A = Aurora::Matrix::fromRawData(dataA, 4, 2);
Aurora::Matrix B = Aurora::Matrix::fromRawData(dataB, 4, 2);
auto C = Aurora::dot(A, B);
EXPECT_EQ(1, C.getDimSize(0));
EXPECT_EQ(2, C.getDimSize(1));
EXPECT_DOUBLE_EQ(12, fourDecimalRound(C.getData()[0]));
EXPECT_DOUBLE_EQ(44, fourDecimalRound(C.getData()[1]));
}
{
double *dataB = new double[20]{1.1, 2.6, 3.8, 6.2,
4.3, 5.7, 6.9, 10.6,
7.1, 8.3, 9.7, 11.2,
17.8, 13.3,26.5 , -7.7,
9.9, 8.2, 6.3, 5.1};
auto B = Aurora::Matrix::fromRawData(dataB, 4, 5);
auto C = Aurora::dot(Aurora::fft(B),Aurora::fft(B));
EXPECT_EQ(1, C.getDimSize(0));
EXPECT_EQ(5, C.getDimSize(1));
EXPECT_DOUBLE_EQ(0.2434, fourDecimalRound(C.getData()[0]/1000));
EXPECT_DOUBLE_EQ(0., fourDecimalRound(C.getData()[1]/1000));
EXPECT_DOUBLE_EQ(0.8438, fourDecimalRound(C.getData()[2]/1000));
EXPECT_DOUBLE_EQ(0., fourDecimalRound(C.getData()[3]/1000));
EXPECT_DOUBLE_EQ(1.3553, fourDecimalRound(C.getData()[4]/1000));
EXPECT_DOUBLE_EQ(0., fourDecimalRound(C.getData()[5]/1000));
EXPECT_DOUBLE_EQ(5.0211, fourDecimalRound(C.getData()[6]/1000));
EXPECT_DOUBLE_EQ(0., fourDecimalRound(C.getData()[7]/1000));
EXPECT_DOUBLE_EQ(0.9238, fourDecimalRound(C.getData()[8]/1000));
EXPECT_DOUBLE_EQ(0., fourDecimalRound(C.getData()[9]/1000));
}
}
TEST_F(Function2D_Test, interp2) {
double* xD= new double[4]{1,2,3,4};
Aurora::Matrix x(std::shared_ptr<double>(xD,std::default_delete<double[]>()),std::vector<int>{4});