diff --git a/src/Function2D.cpp b/src/Function2D.cpp index e4872a5..823e2ac 100644 --- a/src/Function2D.cpp +++ b/src/Function2D.cpp @@ -43,7 +43,7 @@ Aurora::Matrix Aurora::inv(const Aurora::Matrix &aMatrix) { int size = aMatrix.getDataSize(); int *ipiv = new int[aMatrix.getDimSize(0)]; auto result = malloc(size); - cblas_dcopy(size,result, 1,aMatrix.getData(), 1); + cblas_dcopy(size,aMatrix.getData(), 1,result, 1); LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv); LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv); delete[] ipiv; diff --git a/src/Matrix.cpp b/src/Matrix.cpp index ecf0133..a66efce 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -18,10 +18,15 @@ namespace Aurora{ inline Matrix operatorMxA(CalcFuncD aFunc, double aScalar, const Matrix &aMatrix) { double *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex); - aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, output, 1); - if (aMatrix.getValueType() == Complex) { - aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0, output + 1, - 1); + //复数时，+和-需要特别操作，只影响实部 + if (aMatrix.getValueType() == Complex && (aFunc == &vdAddI || aFunc == &vdSubI)) { + aFunc(aMatrix.getDataSize(), aMatrix.getData() , 2, &aScalar, 0, output , + 2); + double zero = 0.0; + aFunc(aMatrix.getDataSize(), aMatrix.getData()+1 , 2, &zero, 0, output+1 , + 2); + } else{ + aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, output, 1); } return Matrix::New(output, aMatrix.getDimSize(0), aMatrix.getDimSize(1), aMatrix.getDimSize(2), aMatrix.getValueType()); @@ -30,13 +35,16 @@ namespace Aurora{ inline Matrix &operatorMxA_RR(CalcFuncD aFunc, double aScalar, Aurora::Matrix &&aMatrix) { std::cout << "use right ref operation" << std::endl; std::cout << "before operation" << std::endl; - aMatrix.printf(); if (aMatrix.getValueType() == Complex) { - aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, + //针对实部操作 + aFunc(aMatrix.getDataSize(), aMatrix.getData(), 2, &aScalar, 0, aMatrix.getData(), - 1); - aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0, - aMatrix.getData() + 1, 1); + 2); + //乘法除法需特别操作，影响虚部 + if (aFunc == &vdDivI || aFunc == &vdMulI){ + aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 2, &aScalar, 0, + aMatrix.getData() + 1, 2); + } } else { aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, aMatrix.getData(), @@ -52,7 +60,13 @@ namespace Aurora{ CalcFuncD aFuncD, const int size, double* xC,double* yN,double *output, int DimsStride) { aFuncD(size, xC, DimsStride * 2, yN, 1, output, 2); - aFuncD(size, xC + 1, DimsStride * 2, yN, 1, output + 1, 2); + if (aFuncD == &vdDivI || aFuncD == &vdMulI){ + aFuncD(size, xC + 1, DimsStride * 2, yN, 1, output + 1, 2); + } + else{ + double zero = 0.0; + aFuncD(size, xC+1 , DimsStride*2, &zero, 0, output + 1, 2); + } } inline double* _MxM_CN_Calc( @@ -66,9 +80,15 @@ namespace Aurora{ inline void V_MxM_NC_Calc( CalcFuncD aFuncD, - const int size, double* xC,double* yN,double *output, int DimsStride) { - aFuncD(size, xC, DimsStride, yN, 2, output, 2); - aFuncD(size, xC , DimsStride, yN+ 1, 2, output + 1, 2); + const int size, double* xN,double* yC,double *output, int DimsStride) { + aFuncD(size, xN, DimsStride, yC, 2, output, 2); + if (aFuncD == &vdDivI || aFuncD == &vdMulI){ + aFuncD(size, xN , DimsStride, yC+ 1, 2, output + 1, 2); + } + else{ + double zero = 0.0; + aFuncD(size, yC+ 1, 2, &zero, 0, output + 1, 2); + } } inline double* _MxM_NC_Calc( @@ -290,7 +310,6 @@ namespace Aurora { if (slices > 0)vector.push_back(slices); Matrix ret({data, free}, vector); if (type != Normal)ret.setValueType(type); - ret.printf(); return ret; } diff --git a/test/FunctionTester.cpp b/test/FunctionTester.cpp index aa65dca..ed86582 100644 --- a/test/FunctionTester.cpp +++ b/test/FunctionTester.cpp @@ -2,6 +2,7 @@ #include "Matrix.h" #include "Function.h" #include "Function1D.h" +#include "Function2D.h" #define DISPLAY_MATRIX(Matrix)\ @@ -10,6 +11,9 @@ printf("%s:\r\n", #Matrix);\ Matrix.printf();\ printf("%s end================================\r\n", #Matrix); +#define DISPLAY_MATRIX + + class FunctionTester:public ::testing::Test{ protected: static void SetUpFunctionTester(){ @@ -29,14 +33,12 @@ double fourDecimalRound(double src){ TEST_F(FunctionTester, MatrixCreate) { - printf("1"); double * dataA =Aurora::malloc(9); double * dataB = new double[9]; double * dataC = new double[9]; double * dataD = new double[9]; //mkl matrix { - printf("2"); Aurora::Matrix A = Aurora::Matrix::New(dataA, 3, 3); EXPECT_EQ(dataA, A.getData()); EXPECT_EQ(9, A.getDataSize()); @@ -44,7 +46,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(3, A.getDimSize(0)); EXPECT_EQ(3, A.getDimSize(1)); EXPECT_EQ(Aurora::Normal, A.getValueType()); - DISPLAY_MATRIX(A) + DISPLAY_MATRIX(A); } // double [] matrix { @@ -55,7 +57,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(3, B.getDimSize(0)); EXPECT_EQ(3, B.getDimSize(1)); EXPECT_EQ(Aurora::Normal, B.getValueType()); - DISPLAY_MATRIX(B) + DISPLAY_MATRIX(B); } // copy from double [] to mkl matrix { @@ -66,7 +68,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(3, C.getDimSize(0)); EXPECT_EQ(3, C.getDimSize(1)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } // 2vector { @@ -77,7 +79,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(9, C.getDimSize(0)); EXPECT_EQ(1, C.getDimSize(1)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } // 2d vector { @@ -88,7 +90,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(9, C.getDimSize(0)); EXPECT_EQ(1, C.getDimSize(1)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } // 2d vector column major { @@ -100,7 +102,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(1, C.getDimSize(0)); EXPECT_EQ(9, C.getDimSize(1)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } // 3d matrix { @@ -113,7 +115,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(3, C.getDimSize(1)); EXPECT_EQ(1, C.getDimSize(2)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } // 3d matrix 2 { @@ -126,7 +128,7 @@ TEST_F(FunctionTester, MatrixCreate) { EXPECT_EQ(1, C.getDimSize(1)); EXPECT_EQ(3, C.getDimSize(2)); EXPECT_EQ(Aurora::Normal, C.getValueType()); - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); } } @@ -140,14 +142,8 @@ TEST_F(FunctionTester, matrixSlice) { dataC[i]=(double)(9-i); } Aurora::Matrix A = Aurora::Matrix::New(dataA,2,2,2); - printf("A:\r\n"); - A.printf(); Aurora::Matrix B = Aurora::Matrix::New(dataB,2,2,2); - printf("B:\r\n"); - B.printf(); Aurora::Matrix C = Aurora::Matrix::New(dataC,2,2,2); - printf("C:\r\n"); - C.printf(); //2D slice EXPECT_EQ(4.0,dataA[4]); A(Aurora::$,Aurora::$,1) = B(0,Aurora::$,Aurora::$); @@ -175,23 +171,122 @@ TEST_F(FunctionTester, matrixOpertaor) { { double * dataA =new double[8]; double * dataB =new double[8]; + double * dataD =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,2,2,2); - DISPLAY_MATRIX(A) + DISPLAY_MATRIX(A); Aurora::Matrix B = Aurora::Matrix::fromRawData(dataB,2,2,2); - DISPLAY_MATRIX(B) + DISPLAY_MATRIX(B); + Aurora::Matrix D = Aurora::Matrix::fromRawData(dataD,4,2); + DISPLAY_MATRIX(D); auto C = (A*B)-B; - DISPLAY_MATRIX(C) + DISPLAY_MATRIX(C); EXPECT_EQ(C.getData()[2],2); C = A*B/2.0; EXPECT_EQ(C.getData()[2],2); C = (A*8.0/B) * (B+1.0)-2.0+8; EXPECT_EQ(C.getData()[2],30); + //Error matrix + C = A*D; + EXPECT_EQ(C.getDataSize(),0); } + //complex + { + double * dataA =Aurora::malloc(8,true); + double * dataB =Aurora::malloc(8,true); + double * dataD =new double[8]; + for (int i = 0; i < 16; ++i) { + dataA[i]=(double)(i+1); + dataB[i]=(double)(2); + if(i<8) dataD[i]=(double)(2); + } + Aurora::Matrix A = Aurora::Matrix::New(dataA,2,2,2,Aurora::Complex); + DISPLAY_MATRIX(A); + Aurora::Matrix B = Aurora::Matrix::New(dataB,2,2,2,Aurora::Complex); + DISPLAY_MATRIX(B); + Aurora::Matrix D = Aurora::Matrix::fromRawData(dataD,2,2,2); + DISPLAY_MATRIX(D); + auto C = A + 5.0 -2.0; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],6); + EXPECT_EQ(C.getData()[3],4); + EXPECT_TRUE(C.isComplex()); + C = A*B-B; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],-4); + EXPECT_EQ(C.getData()[3],12); + EXPECT_TRUE(C.isComplex()); + C = A*B-2.0+5.0; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],1); + EXPECT_EQ(C.getData()[3],14); + C = A*B/2.0*3.; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],-3); + EXPECT_EQ(C.getData()[3],21); + C = A*B/B; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],3); + EXPECT_EQ(C.getData()[3],4); + C = A*B/(B*A/A); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],3); + EXPECT_EQ(C.getData()[3],4); + C = B*D; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],4); + C = B/D; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],1); + EXPECT_EQ(C.getData()[3],1); + C = B+D; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],2); + C = B-D; + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],0); + EXPECT_EQ(C.getData()[3],2); + + C = B*(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],4); + C = B/(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],1); + EXPECT_EQ(C.getData()[3],1); + C = B+(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],2); + C = B-(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],0); + EXPECT_EQ(C.getData()[3],2); + + C = (B*1.0)*(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],4); + C = (B*1.0)/(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],1); + EXPECT_EQ(C.getData()[3],1); + C = (B*1.0)+(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],4); + EXPECT_EQ(C.getData()[3],2); + C = (B*1.0)-(D*1.0); + DISPLAY_MATRIX(C); + EXPECT_EQ(C.getData()[2],0); + EXPECT_EQ(C.getData()[3],2); + } } TEST_F(FunctionTester, sign) { @@ -202,20 +297,52 @@ TEST_F(FunctionTester, sign) { dataB[i]=(double)(i+2); } Aurora::Matrix A = Aurora::Matrix::New(dataA,3,3); - printf("A:\r\n"); - A.printf(); Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3); - printf("B:\r\n"); - B.printf(); - printf("sign(A):\r\n"); - sign(A*B).printf(); + auto C = sign(A*B); + double * result = C.getData(); + EXPECT_EQ(-1, result[0]); + EXPECT_EQ(0, result[3]); + EXPECT_EQ(1, result[4]); } TEST_F(FunctionTester, immse){ - double dataA[9]={1,2,3,4,5,6,7,8,9}; - double dataB[9]={10,20,30,40,50,40,30,20,10}; - EXPECT_DOUBLE_EQ(fourDecimalRound(Aurora::immse(dataA,dataB,9)),698.3333)<<" immse error;"; + double *dataA = new double[9]{1,2,3,4,5,6,7,8,9}; + double *dataB = new double[9]{10,20,30,40,50,40,30,20,10}; + Aurora::Matrix A = Aurora::Matrix::fromRawData(dataA,3,3); + Aurora::Matrix B = Aurora::Matrix::fromRawData(dataB,3,3); + double v = immse(A,B); + EXPECT_DOUBLE_EQ(fourDecimalRound(v),698.3333)<<" immse error;"; +} + +TEST_F(FunctionTester, inv){ + //默认是column major排布数据 + double *dataA=new double[9]{2, 0, 2, 2, 3, 0, 3, 3, 3}; + double *dataB=new double[9]{1, 1, 1, 1, 1, 1, 1, 1, 1}; + Aurora::Matrix A = Aurora::Matrix::fromRawData(dataA,3,3); + Aurora::Matrix B = Aurora::Matrix::fromRawData(dataB,3,3); + auto invA = Aurora::inv(A); + double* result = invA.getData(); + EXPECT_DOUBLE_EQ(0.75, fourDecimalRound(result[0])); + EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[1])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[2])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[3])); + EXPECT_DOUBLE_EQ(.0, fourDecimalRound(result[4])); + EXPECT_DOUBLE_EQ(0.3333, fourDecimalRound(result[5])); + EXPECT_DOUBLE_EQ(-0.25, fourDecimalRound(result[6])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[7])); + EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[8])); + invA = Aurora::inv(A*B); + result = invA.getData(); + EXPECT_DOUBLE_EQ(0.75, fourDecimalRound(result[0])); + EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[1])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[2])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[3])); + EXPECT_DOUBLE_EQ(.0, fourDecimalRound(result[4])); + EXPECT_DOUBLE_EQ(0.3333, fourDecimalRound(result[5])); + EXPECT_DOUBLE_EQ(-0.25, fourDecimalRound(result[6])); + EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[7])); + EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[8])); } TEST_F(FunctionTester, polyval){ @@ -258,22 +385,6 @@ TEST_F(FunctionTester, fftAndComplexAndIfft){ delete [] ifftResult; } -TEST_F(FunctionTester, inv){ - //默认是column major排布数据 - double dataMA[9]={2, 0, 2, 2, 3, 0, 3, 3, 3}; - double* result = Aurora::inv(3,dataMA); - EXPECT_DOUBLE_EQ(0.75, fourDecimalRound(result[0]))<<" inv value error"; - EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[1]))<<" inv value error"; - EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[2]))<<" inv value error"; - EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[3]))<<" inv value error"; - EXPECT_DOUBLE_EQ(.0, fourDecimalRound(result[4]))<<" inv value error"; - EXPECT_DOUBLE_EQ(0.3333, fourDecimalRound(result[5]))<<" inv value error"; - EXPECT_DOUBLE_EQ(-0.25, fourDecimalRound(result[6]))<<" inv value error"; - EXPECT_DOUBLE_EQ(-0.5, fourDecimalRound(result[7]))<<" inv value error"; - EXPECT_DOUBLE_EQ(0.5, fourDecimalRound(result[8]))<<" inv value error"; - delete [] result; -} - TEST_F(FunctionTester, hilbert) { double input[10]{1,1,0,2,2,0,1,1,0,2}; auto result = Aurora::hilbert(10,input);