Fix matrix operator bug on - and /. 3

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;

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);
+        //复数时，+和-需要特别操作，只影响实部
+        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);
-        if (aMatrix.getValueType() == Complex) {
-            aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0, output + 1,
-                  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);
+                  2);
-            aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0,
+            //乘法除法需特别操作，影响虚部
-                  aMatrix.getData() + 1, 1);
+            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,8 +60,14 @@ namespace Aurora{
             CalcFuncD aFuncD,
             const int size, double* xC,double* yN,double *output, int DimsStride) {
         aFuncD(size, xC, DimsStride * 2, yN, 1, output, 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(
             CalcFuncD aFuncD,
@@ -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) {
+            const int size, double* xN,double* yC,double *output, int DimsStride) {
-        aFuncD(size, xC, DimsStride, yN, 2, output, 2);
+        aFuncD(size, xN, DimsStride, yC, 2, output, 2);
-        aFuncD(size, xC , DimsStride, yN+ 1, 2, output + 1, 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;
     }
 

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");
+    auto C = sign(A*B);
-    B.printf();
+    double * result = C.getData();
-    printf("sign(A):\r\n");
+    EXPECT_EQ(-1, result[0]);
-    sign(A*B).printf();
+    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 *dataA = new double[9]{1,2,3,4,5,6,7,8,9};
-    double dataB[9]={10,20,30,40,50,40,30,20,10};
+    double *dataB = new double[9]{10,20,30,40,50,40,30,20,10};
-    EXPECT_DOUBLE_EQ(fourDecimalRound(Aurora::immse(dataA,dataB,9)),698.3333)<<" immse error;";
+    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);