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

Function1D bug fix and unit test.

Browse Source
This commit is contained in:
Krad
2023-04-23 10:39:26 +08:00
parent 0cde4460a8
commit 5e85d0f361
3 changed files with 112 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/Function1D.cpp
View File

@@ -26,7 +26,7 @@ Aurora::Matrix Aurora::complex(const Aurora::Matrix &matrix) {
std::cerr<<"complex not support complex value type"<<std::endl;
return matrix;
}
auto output = (std::complex<double> *) mkl_malloc(matrix.getDataSize() * sizeof(std::complex<double>), 64);
auto output = malloc(matrix.getDataSize() ,true);
memset(output, 0, (matrix.getDataSize() * sizeof(std::complex<double>)));
cblas_dcopy(matrix.getDataSize(), matrix.getData(), REAL_STRIDE, (double *) output, COMPLEX_STRIDE);
return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2),
@@ -38,7 +38,7 @@ Aurora::Matrix Aurora::real(const Aurora::Matrix &matrix) {
std::cerr<<"real only support complex value type"<<std::endl;
return matrix;
}
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = (double *) malloc(matrix.getDataSize());
memset(output, 0, (matrix.getDataSize() * sizeof(double)));
cblas_dcopy(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE , (double *) output, REAL_STRIDE);
return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -49,14 +49,14 @@ Aurora::Matrix Aurora::imag(const Aurora::Matrix &matrix) {
std::cerr<<"imag only support complex value type"<<std::endl;
return matrix;
}
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
memset(output, 0, (matrix.getDataSize() * sizeof(double)));
cblas_dcopy(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE , (double *) output, REAL_STRIDE);
return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
}
Aurora::Matrix Aurora::ceil(const Aurora::Matrix &matrix) {
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
//for real part
vdCeilI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
if (matrix.getValueType() == Complex) {
@@ -78,7 +78,7 @@ Aurora::Matrix Aurora::ceil(const Aurora::Matrix &&matrix) {
}
Aurora::Matrix Aurora::round(const Aurora::Matrix &matrix) {
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
//for real part
vdRoundI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
if (matrix.getValueType() == Complex) {
@@ -102,7 +102,7 @@ Aurora::Matrix Aurora::round(const Aurora::Matrix &&matrix) {
Aurora::Matrix Aurora::sqrt(const Aurora::Matrix& matrix) {
if (matrix.getValueType() != Complex) {
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
vdSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
return Aurora::Matrix::New(output, matrix);
}
@@ -110,7 +110,7 @@ Aurora::Matrix Aurora::sqrt(const Aurora::Matrix& matrix) {
return Aurora::Matrix();
}
Aurora::Matrix Aurora::sqrt(const Aurora::Matrix&& matrix) {
Aurora::Matrix Aurora::sqrt(Aurora::Matrix&& matrix) {
std::cout<<"RR sqrt"<<std::endl;
if (matrix.getValueType() != Complex) {
vdSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
@@ -121,7 +121,7 @@ Aurora::Matrix Aurora::sqrt(const Aurora::Matrix&& matrix) {
}
Aurora::Matrix Aurora::abs(const Aurora::Matrix &matrix) {
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
if (matrix.getValueType()==Normal){
vdAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
}
@@ -131,7 +131,7 @@ Aurora::Matrix Aurora::abs(const Aurora::Matrix &matrix) {
return Aurora::Matrix::New(output, matrix);
}
Aurora::Matrix Aurora::abs(const Aurora::Matrix&& matrix) {
Aurora::Matrix Aurora::abs(Aurora::Matrix&& matrix) {
std::cout<<"RR abs"<<std::endl;
if (matrix.getValueType()==Normal){
vdAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
@@ -139,7 +139,7 @@ Aurora::Matrix Aurora::abs(const Aurora::Matrix&& matrix) {
}
//TODO：考虑尝试是不是使用realloc缩短已分配的内存的方式重用matrix
else{
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(double), 64);
auto output = malloc(matrix.getDataSize());
vzAbsI(matrix.getDataSize(), (std::complex<double> *)matrix.getData(), SAME_STRIDE,output, SAME_STRIDE);
return Aurora::Matrix::New(output, matrix);
}
@@ -155,7 +155,7 @@ Aurora::Matrix Aurora::sign(const Aurora::Matrix &matrix) {
}
else{
//sign(x) = x./abs(x)，前提是 x 为复数。
auto output = (double *) mkl_malloc(matrix.getDataSize() * sizeof(std::complex<double>), 64);
auto output = malloc(matrix.getDataSize(),true);
Matrix absMatrix = abs(matrix);
vdDivI(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE,
absMatrix.getData(), REAL_STRIDE,output,COMPLEX_STRIDE);
@@ -165,7 +165,7 @@ Aurora::Matrix Aurora::sign(const Aurora::Matrix &matrix) {
}
}
Aurora::Matrix Aurora::sign(const Aurora::Matrix&& matrix) {
Aurora::Matrix Aurora::sign(Aurora::Matrix&& matrix) {
std::cout<<"RR sign"<<std::endl;
if (matrix.getValueType()==Normal){
Eigen::Map<Eigen::VectorXd> retV(matrix.getData(),matrix.getDataSize());

6
src/Function1D.h
View File

@@ -31,15 +31,15 @@ namespace Aurora {
*/
Matrix sqrt(const Matrix& matrix);
Matrix sqrt(const Matrix&& matrix);
Matrix sqrt(Matrix&& matrix);
Matrix abs(const Matrix& matrix);
Matrix abs(const Matrix&& matrix);
Matrix abs(Matrix&& matrix);
Matrix sign(const Matrix& matrix);
Matrix sign(const Matrix&& matrix);
Matrix sign(Matrix&& matrix);
Matrix interp1(const Matrix& aX, const Matrix& aV, const Matrix& aX1, InterpnMethod aMethod);

97
test/Function1D_Test.cpp
View File

@@ -1,5 +1,8 @@
#include <gtest/gtest.h>
#include <vector>
#include "TestUtility.h"
#include "Matrix.h"
#include "Function.h"
#include "Function1D.h"
@@ -107,3 +110,97 @@ TEST_F(Function1D_Test, polyval){
EXPECT_DOUBLE_EQ(262., resultP[2])<<" polyval error;";
delete [] resultP;
}
TEST_F(Function1D_Test, complexAndEtc){
//complex
{
double *dataP =new double[8]{3,2,1,6, 7, 8 , 19, 13};
auto matrixAN = Aurora::Matrix::fromRawData(dataP,2,4);
auto matrixBC = Aurora::complex(matrixAN);
EXPECT_EQ(8, matrixBC.getDataSize());
EXPECT_EQ(Aurora::Complex, matrixBC.getValueType());
EXPECT_DOUBLE_EQ(3, matrixBC.getData()[0]);
EXPECT_DOUBLE_EQ(0, matrixBC.getData()[3]);
}
//complex & real & imag
{
double *dataP =new double[8]{3,2,1,6, 7, 8 , 19, 13};
auto matrixAC = Aurora::Matrix::fromRawData(dataP,2,2,0,Aurora::Complex);
auto matrixCN = Aurora::real(matrixAC);
EXPECT_EQ(4, matrixCN.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixCN.getValueType());
EXPECT_DOUBLE_EQ(3, matrixCN.getData()[0]);
EXPECT_DOUBLE_EQ(19, matrixCN.getData()[3]);
auto matrixDN = Aurora::imag(matrixAC);
EXPECT_EQ(4, matrixDN.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixDN.getValueType());
EXPECT_DOUBLE_EQ(2, matrixDN.getData()[0]);
EXPECT_DOUBLE_EQ(13, matrixDN.getData()[3]);
}
}
TEST_F(Function1D_Test, ceilAndRound) {
double *dataP = new double[3]{3.1, 2.5, 1.8};
auto matrixA = Aurora::Matrix::fromRawData(dataP, 3);
auto matrixB = Aurora::ceil(matrixA);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(4, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(3, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[2]);
matrixB = Aurora::ceil(matrixA * 0.5) ;
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(2, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(1, matrixB.getData()[2]);
matrixB = Aurora::round(matrixA);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(3, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(3, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[2]);
matrixB = Aurora::round(matrixA * 0.5);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(2, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(1, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(1, matrixB.getData()[2]);
}
TEST_F(Function1D_Test, absAndSqrt) {
double *dataP = new double[3]{1, 4, -3};
auto matrixA = Aurora::Matrix::fromRawData(dataP, 3);
auto matrixB = Aurora::abs(matrixA);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(1, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(4, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(3, matrixB.getData()[2]);
matrixB = Aurora::sqrt(matrixB);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(1, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(1.7321, fourDecimalRound(matrixB.getData()[2]));
matrixB = Aurora::sqrt(Aurora::abs(matrixA*-1));
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(1, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[1]);
EXPECT_DOUBLE_EQ(1.7321, fourDecimalRound(matrixB.getData()[2]));
matrixB = Aurora::sqrt(matrixA);
EXPECT_EQ(3, matrixB.getDataSize());
EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
EXPECT_DOUBLE_EQ(1, matrixB.getData()[0]);
EXPECT_DOUBLE_EQ(2, matrixB.getData()[1]);
EXPECT_TRUE(isnanf(matrixB.getData()[2]));
}