From cf96572074141e6c6390bde25fdf6887658fd290 Mon Sep 17 00:00:00 2001 From: Krad Date: Tue, 25 Apr 2023 17:28:33 +0800 Subject: [PATCH] Add prod, dot and their unit test. --- src/Function2D.cpp | 139 +++++++++++++++++++++++++++++++-------- src/Function2D.h | 11 +++- test/Function2D_Test.cpp | 56 ++++++++++++++++ 3 files changed, 179 insertions(+), 27 deletions(-) diff --git a/src/Function2D.cpp b/src/Function2D.cpp index f12c08c..aa78bf9 100644 --- a/src/Function2D.cpp +++ b/src/Function2D.cpp @@ -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 srcV(aMatrix.getData(),aMatrix.getDataSize()); - double * ret = malloc(1); - ret[0] = srcV.array().sum(); - return Matrix::New(ret,1); + case All: + { + Eigen::Map srcV((std::complex*)aMatrix.getData(),aMatrix.getDataSize()); + std::complex* ret = (std::complex*)malloc(1,true); + ret[0] = srcV.array().sum(); + return Matrix::New((double*)ret,1,1,1,Complex); + } + case Row: + { + Eigen::Map srcM((std::complex*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + std::complex * ret = (std::complex*)malloc(aMatrix.getDimSize(0),true); + Eigen::Map 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 srcM((std::complex*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + std::complex* ret = (std::complex*)malloc(aMatrix.getDimSize(1),true); + Eigen::Map 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 srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); - double * ret = malloc(aMatrix.getDimSize(0)); - Eigen::Map retV(ret,aMatrix.getDimSize(0)); - retV = srcM.rowwise().sum(); - return Matrix::New(ret,aMatrix.getDimSize(0),1); - } - case Column: - default: - { - Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); - double * ret = malloc(aMatrix.getDimSize(1)); - Eigen::Map retV(ret,aMatrix.getDimSize(0)); - retV = srcM.colwise().sum(); - return Matrix::New(ret,1,aMatrix.getDimSize(1)); + case All: + { + Eigen::Map srcV(aMatrix.getData(),aMatrix.getDataSize()); + double * ret = malloc(1); + ret[0] = srcV.array().sum(); + return Matrix::New(ret,1); + } + case Row: + { + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * ret = malloc(aMatrix.getDimSize(0)); + Eigen::Map retV(ret,aMatrix.getDimSize(0)); + retV = srcM.rowwise().sum(); + return Matrix::New(ret,aMatrix.getDimSize(0),1); + } + case Column: + default: + { + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * ret = malloc(aMatrix.getDimSize(1)); + Eigen::Map 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 srcM((std::complex*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + auto ret = malloc(aMatrix.getDimSize(1),true); + Eigen::Map retV((std::complex*)ret,aMatrix.getDimSize(1)); + retV = srcM.colwise().prod(); + return Matrix::New(ret,1,aMatrix.getDimSize(1),1,Complex); + } + else{ + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + auto ret = malloc(aMatrix.getDimSize(1)); + Eigen::Map 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); + } + } +} + diff --git a/src/Function2D.h b/src/Function2D.h index cc9d74e..f87f60a 100644 --- a/src/Function2D.h +++ b/src/Function2D.h @@ -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); }; diff --git a/test/Function2D_Test.cpp b/test/Function2D_Test.cpp index 5ae3785..0c2ff6e 100644 --- a/test/Function2D_Test.cpp +++ b/test/Function2D_Test.cpp @@ -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(xD,std::default_delete()),std::vector{4});