#include <complex>
#include <cstdio>
#include <gtest/gtest.h>
#include "AuroraDefs.h"
#include "Matrix.h"
#include "CudaMatrix.h"
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
#include "Function3D.h"


class CudaMatrix_Test : public ::testing::Test {
protected:
    static void SetUpCudaMatrixTester() {

    }

    static void TearDownTestCase() {
    }

    void SetUp() {
    }

    void TearDown() {
    }
};

TEST_F(CudaMatrix_Test, MatrixAddScalar) {
    {
        auto A = Aurora::zeros(257,257,1);
        auto B = Aurora::zeros(257,257,1);
        size_t size = 257*257;
        for (size_t i = 0; i < size; i++)
        {
            A[i] = -1;
            B[i] =  i;
        }
        printf("Test CudaMatrix operator+(const CudaMatrix &aMatrix) const \r\n");
        //CudaMatrix operator+(const CudaMatrix &aMatrix) const
        auto C = A+B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        auto dC = (dA+dB);
        auto dhC = dC.toHostMatrix();
        EXPECT_EQ(257, dC.getDimSize(0));
        EXPECT_EQ(257, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));

        for (size_t i = 0; i < size; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dhC[i]);
        }
        printf("Test CudaMatrix operator+(float aScalar) const \r\n");
        //CudaMatrix operator+(float aScalar) const
        auto D = C+0.5;
        auto dD = dC+0.5;
        auto dhD = dD.toHostMatrix();
        EXPECT_EQ(257, dD.getDimSize(0));
        EXPECT_EQ(257, dD.getDimSize(1));
        EXPECT_EQ(1, dD.getDimSize(2));
        for (size_t i = 0; i < size; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix operator+(float aScalar, const CudaMatrix &aMatrix) \r\n");
        // CudaMatrix operator+(float aScalar, const CudaMatrix &aMatrix)
        dD = 0.5 + dC;
        dhD = dD.toHostMatrix();
        EXPECT_EQ(257, dD.getDimSize(0));
        EXPECT_EQ(257, dD.getDimSize(1));
        EXPECT_EQ(1, dD.getDimSize(2));
        for (size_t i = 0; i < size; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix &operator+(float aScalar, CudaMatrix &&aMatrix) \r\n");
        // CudaMatrix &operator+(float aScalar, CudaMatrix &&aMatrix)
        {
            auto dD2 = 0.5 + (dA+dB);
            dhD = dD2.toHostMatrix();
            EXPECT_EQ(257, dD2.getDimSize(0));
            EXPECT_EQ(257, dD2.getDimSize(1));
            EXPECT_EQ(1, dD2.getDimSize(2));
            for (size_t i = 0; i < size; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        printf("Test CudaMatrix &operator+(CudaMatrix &&aMatrix, float aScalar) \r\n");
        // CudaMatrix &operator+(CudaMatrix &&aMatrix, float aScalar)
        {
            
            auto dD2 =  (dA+dB)+0.5;
            dhD = dD2.toHostMatrix();
            EXPECT_EQ(257, dD2.getDimSize(0));
            EXPECT_EQ(257, dD2.getDimSize(1));
            EXPECT_EQ(1, dD2.getDimSize(2));
            for (size_t i = 0; i < size; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix) const
        printf("Test CudaMatrix operator+(CudaMatrix &&aMatrix) const \r\n");
        {
            auto D = A+C;
            auto dD2 = dA+(dA+dB);
            dhD = dD2.toHostMatrix();
            EXPECT_EQ(257, dD2.getDimSize(0));
            EXPECT_EQ(257, dD2.getDimSize(1));
            EXPECT_EQ(1, dD2.getDimSize(2));
            for (size_t i = 0; i < size; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix,CudaMatrix &aOther)
        printf("Test CudaMatrix operator+(CudaMatrix &&aMatrix,CudaMatrix &aOther) \r\n");
        {
            auto D = A+C;
            auto dD2 = (dA+dB)+dA;
            dhD = dD2.toHostMatrix();
            EXPECT_EQ(257, dD2.getDimSize(0));
            EXPECT_EQ(257, dD2.getDimSize(1));
            EXPECT_EQ(1, dD2.getDimSize(2));
            for (size_t i = 0; i < size; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
    }
    {
        float *dataA = new float[10000];
        float *dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,100,1,Aurora::Complex);
        auto C = A+B;
        std::complex<float> scalar(-1,-1);
        auto dA = A.toDeviceMatrix();
        //complex matrix + complex scalar
        auto dC = (dA+scalar);
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        
        //complex Matrix&& + complex scalar
        dC = A.toDeviceMatrix()+scalar;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //complex scalar + complex Matrix&&
        dC = scalar+A.toDeviceMatrix();
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //complex scalar + complex matrix
        dC = scalar+dA;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //complex matrix + real scalar
        C = A+1;
        dC = (dA+1);
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index"<<i;
        }
        //complex Matrix&& + real scalar 
        dC = A.toDeviceMatrix()+1;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //real scalar + complex Matrix &&
        dC = 1+A.toDeviceMatrix();
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    {
        float *dataA = new float[5000];
        float *dataB = new float[10000];
        for (size_t i = 0; i < 5000; i++)
        {
            dataA[i] = i;
            dataB[2*i] = -1;
            dataB[2*i+1] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,100,1,Aurora::Complex);
        auto C = A+B;
        std::complex<float> scalar(-1,-1);
        auto dA = A.toDeviceMatrix();
        //real matrix + complex scalar
        auto dC = (dA+scalar);
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        
        //real Matrix&& + complex scalar
        dC = A.toDeviceMatrix()+scalar;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //complex scalar + real Matrix&&
        dC = scalar+A.toDeviceMatrix();
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //complex scalar + real Matrix
        dC = scalar+dA;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(50, dC.getDimSize(0));
        EXPECT_EQ(100, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixAddmm) {
    //real
    {
        auto A = -Aurora::ones(4096,23519,1);
        auto B = Aurora::ones(4096,23519,1)*8;
        auto C = B-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&
        auto dC = (dA+dB);
        EXPECT_TRUE(!dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        #pragma omp parallel for
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&
        dC = A.toDeviceMatrix()+dB;
        EXPECT_TRUE(!dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        // #pragma omp parallel for
        // for (size_t i = 0; i < C.getDataSize(); i++)
        // {
        //     EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        // }
        //Matrix& + Matrix&&
        dC = dB+A.toDeviceMatrix();
        EXPECT_TRUE(!dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        // #pragma omp parallel for
        // for (size_t i = 0; i < C.getDataSize(); i++)
        // {
        //     EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        // }
    }
    //complex
    {
        float *dataA = Aurora::random(4096*23519*2);
        float *dataB = Aurora::random(4096*23519*2);

        auto A = Aurora::Matrix::fromRawData(dataA, 4096,23519,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 4096,23519,1,Aurora::Complex);

        auto C = B+A;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, 
        auto dC = (dA+dB);
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        #pragma omp parallel for
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&,
        dC = A.toDeviceMatrix()+dB;
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        // #pragma omp parallel for
        // for (size_t i = 0; i < C.getDataSize(); i++)
        // {
        //     EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        // }
        //Matrix& + Matrix&&, 
        dC = dB+A.toDeviceMatrix();
        EXPECT_TRUE(dC.isComplex());
        EXPECT_EQ(4096, dC.getDimSize(0));
        EXPECT_EQ(23519, dC.getDimSize(1));
        EXPECT_EQ(1, dC.getDimSize(2));
        // #pragma omp parallel for
        // for (size_t i = 0; i < C.getDataSize(); i++)
        // {
        //     EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        // }
    }
}

TEST_F(CudaMatrix_Test, MatrixAddmv) {
    //real
    {
        auto A = Aurora::zeros(100,100,1);
        auto B = Aurora::ones(100,1,1);
        auto C = A+1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);

        //Matrix& + Matrix&, row mode
        dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + Matrix&&, row mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + vec&&, row mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataC[i] = -1;
            if (i<100)dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto C = A+D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        B.forceReshape(1, 50, 1);
        A.forceReshape(100, 50, 1);
        dA.forceReshape(100, 50, 1);

        //Matrix& + Matrix&, row mode
        dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + Matrix&&, row mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + vec&&, row mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complexM - realv
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            if (i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,1,1);

        auto C = A-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + vec&, col mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec& + Matrix&&, col mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&, col mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        dA.forceReshape(50, 100, 1);
        A.forceReshape(50, 100, 1);
        B.forceReshape(1, 100, 1);

        //Matrix& + vec&, row mode
        dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + vec&, row mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec& + Matrix&&, row mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, row mode
        dC = dA+B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&, row mode
        dC = B.toDeviceMatrix()+dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //realM - complexV
    {
        float *dataA = new float[5000];
        float *dataC = new float[10000];
        float *dataB = new float[100];
        for (size_t i = 0; i < 10000; i++)
        {
            if (i<5000)dataA[i] = 8;
            dataC[i] = -1;
            if (i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);


        auto C = A+D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        dA.forceReshape(100, 50, 1);
        A.forceReshape(100, 50, 1);
        B.forceReshape(1, 50, 1);
        //Matrix& + Matrix&, row mode
        dC = (dA+dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()+dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, row mode
        dC = dB+A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixMul) {
    {
        auto A = Aurora::zeros(1000,1,1);
        auto B = Aurora::zeros(1000,1,1);
        for (size_t i = 0; i < 1000; i++)
        {
            A[i] = -1;
            B[i] =  i;
        }
        printf("Test CudaMatrix operator*(const CudaMatrix &aMatrix) const \r\n");
        //CudaMatrix operator+(const CudaMatrix &aMatrix) const
        auto C = A*B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        auto dC = (dA*dB);
        auto dhC = dC.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dhC[i]);
        }
        printf("Test CudaMatrix operator*(float aScalar) const \r\n");
        //CudaMatrix operator+(float aScalar) const
        auto D = C*0.5;
        auto dD = dC*0.5;
        auto dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix operator*(float aScalar, const CudaMatrix &aMatrix) \r\n");
        // CudaMatrix operator+(float aScalar, const CudaMatrix &aMatrix)
        dD = 0.5 * dC;
        dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix &operator*(float aScalar, CudaMatrix &&aMatrix) \r\n");
        // CudaMatrix &operator+(float aScalar, CudaMatrix &&aMatrix)
        {
            auto dD2 = 0.5 * (dA*dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        printf("Test CudaMatrix &operator*(CudaMatrix &&aMatrix, float aScalar) \r\n");
        // CudaMatrix &operator+(CudaMatrix &&aMatrix, float aScalar)
        {
            
            auto dD2 =  (dA*dB)*0.5;
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix) const
        printf("Test CudaMatrix operator*(CudaMatrix &&aMatrix) const \r\n");
        {
            auto D = A*C;
            auto dD2 = dA*(dA*dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator*(CudaMatrix &&aMatrix,CudaMatrix &aOther)
        printf("Test CudaMatrix operator*(CudaMatrix &&aMatrix,CudaMatrix &aOther) \r\n");
        {
            auto D = A*C;
            auto dD2 = (dA*dB)*dA;
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
    }
    {
        float *dataA = new float[100];
        float *dataB = new float[100];
        for (size_t i = 0; i < 100; i++)
        {
            dataA[i] = 8;
            dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 10,5,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 10,5,1,Aurora::Complex);
        auto C = A*B;
        std::complex<float> scalar(-1,-1);
        auto dA = A.toDeviceMatrix();

        auto dC = (dA*scalar);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = A.toDeviceMatrix()*scalar;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = scalar*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = A*2;
        dC = (dA*2);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index"<<i;
        }

        dC = A.toDeviceMatrix()*2;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = 2*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixMulmm) {
    //real
    {
        auto A = Aurora::zeros(100,100,1);
        auto B = Aurora::zeros(100,100,1);
        for (size_t i = 0; i < 10000; i++)
        {
            A[i] = -1;
            B[i] =  8;
        }
        auto C = B*-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& * Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& * Matrix&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& * Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,50,1,Aurora::Complex);

        auto C = B*A;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixMulmv) {
    //real
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            if(i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,1,1);
        auto C = A*(-1);
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);

        //Matrix& + Matrix&, row mode
        dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + Matrix&&, row mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + vec&&, row mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataC[i] = -1;
            if (i<100)dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto C = A*D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        B.forceReshape(1, 50, 1);
        A.forceReshape(100, 50, 1);
        dA.forceReshape(100, 50, 1);

        //Matrix& + Matrix&, row mode
        dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + Matrix&&, row mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& + vec&&, row mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complexM - realv
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            if (i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,1,1);

        auto C = A*(-1);
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + vec&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec& * Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&, col mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        dA.forceReshape(50, 100, 1);
        A.forceReshape(50, 100, 1);
        B.forceReshape(1, 100, 1);

        //Matrix& + vec&, row mode
        dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + vec&, row mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec& + Matrix&&, row mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, row mode
        dC = dA*B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&, row mode
        dC = B.toDeviceMatrix()*dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //realM - complexV
    {
        float *dataA = new float[5000];
        float *dataC = new float[10000];
        float *dataB = new float[100];
        for (size_t i = 0; i < 10000; i++)
        {
            if (i<5000)dataA[i] = 8;
            dataC[i] = -1;
            if (i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);


        auto C = A*D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, col mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        dA.forceReshape(100, 50, 1);
        A.forceReshape(100, 50, 1);
        B.forceReshape(1, 50, 1);
        //Matrix& + Matrix&, row mode
        dC = (dA*dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()*dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + Matrix&&, row mode
        dC = dB*A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixSub) {
    {
        auto A = Aurora::zeros(1000,1,1);
        auto B = Aurora::zeros(1000,1,1);
        for (size_t i = 0; i < 1000; i++)
        {
            A[i] = -1;
            B[i] =  i;
        }
        printf("Test CudaMatrix operator-(const CudaMatrix &aMatrix) const \r\n");
        //CudaMatrix operator+(const CudaMatrix &aMatrix) const
        auto C = A-B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        auto dC = (dA-dB);
        auto dhC = dC.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dhC[i]);
        }
        printf("Test CudaMatrix operator-(float aScalar) const \r\n");
        //CudaMatrix operator+(float aScalar) const
        auto D = C-0.5;
        auto dD = dC-0.5;
        auto dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix operator-(float aScalar, const CudaMatrix &aMatrix) \r\n");
        // CudaMatrix operator+(float aScalar, const CudaMatrix &aMatrix)
        D = 0.5 - C;
        dD = 0.5 - C.toDeviceMatrix();
        dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix &operator-(float aScalar, CudaMatrix &&aMatrix) \r\n");
        // CudaMatrix &operator+(float aScalar, CudaMatrix &&aMatrix)
        {
            auto dD2 = 0.5 - (dA-dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        D =  C - 0.5;
        printf("Test CudaMatrix &operator-(CudaMatrix &&aMatrix, float aScalar) \r\n");
        // CudaMatrix &operator+(CudaMatrix &&aMatrix, float aScalar)
        {
            
            auto dD2 =  (dA-dB)-0.5;
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix) const
        printf("Test CudaMatrix operator-(CudaMatrix &&aMatrix) const \r\n");
        {
            auto D = A-C;
            auto dD2 = dA-(dA-dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix,CudaMatrix &aOther)
        printf("Test CudaMatrix operator-(CudaMatrix &&aMatrix,CudaMatrix &aOther) \r\n");
        {
            auto D = A-C;
            auto dD2 = dA- (dA-dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
    }
    {
        float *dataA = new float[10000];
        float *dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,100,1,Aurora::Complex);
        auto C = A-B;
        std::complex<float> scalar(-1,-1);
        auto dA = A.toDeviceMatrix();

        auto dC = (dA-scalar);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = A.toDeviceMatrix()-scalar;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B-A;
        dC = scalar-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = A-1;
        dC = (dA-1);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index"<<i;
        }

        dC = A.toDeviceMatrix()-1;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 1-A;
        dC = 1-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixSubmm) {
    //real
    {
        auto A = Aurora::zeros(100,100,1);
        auto B = Aurora::zeros(100,100,1);
        for (size_t i = 0; i < 10000; i++)
        {
            A[i] = 1;
            B[i] =  8;
        }
        auto C = B-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();

        auto dC = (dB-dA);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = dB - A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        auto dataA = new float[10000];
        auto dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 1;
            dataB[i] = 8;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,50,1,Aurora::Complex);

        auto C = A-B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B-A;
        //Matrix& + Matrix&&,
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //real & complex
    {
        auto dataA = new float[10000];
        auto dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 1;
            dataB[i] = 8;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,50,1,Aurora::Complex);

        auto C = A-B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix&(real) - Matrix&(complex)
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i<<", A:"<<A[i]<<", B:"<<B[i/2];
        }
        //Matrix&&(real) - Matrix&(complex)
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&(real) - Matrix&&(complex)
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B-A;
        //Matrix&(complex) - Matrix&&(real)
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&&(complex) - Matrix&(real)
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&(complex) - Matrix&(real)
        dC = dB-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixSubmv) {
    //real
    {
        auto A = Aurora::ones(100,100,1)*8;
        auto B = Aurora::ones(100,1,1);
        auto C = A-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 1- A;
        //Matrix& + Matrix&&, col mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        C = A-1;
        //Matrix& - Matrix&, row mode
        dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& - Matrix&, row mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& - vec&&, row mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 1- A;
        //Matrix& + Matrix&&, row mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataC[i] = -1;
            if (i<100)dataB[i] = -1;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto C = A-D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = D-A;
        //Matrix& + Matrix&&, col mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        B.forceReshape(1, 50, 1);
        A.forceReshape(100, 50, 1);
        D.forceReshape(100, 50, 1);
        dA.forceReshape(100, 50, 1);

        C= A-D;
        //Matrix& - Matrix&, row mode
        dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& - Matrix&, row mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& - vec&&, row mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        C= D - A;
        //Matrix& + Matrix&&, row mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complexM - realv
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            if (i<100)dataB[i] = 1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);

        auto B = Aurora::Matrix::fromRawData(dataB, 100,1,1);

        auto C = A-1;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + vec&, col mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C=1-A;
        //vec& + Matrix&&, col mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& + Matrix&, col mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        dA.forceReshape(50, 100, 1);
        A.forceReshape(50, 100, 1);
        B.forceReshape(1, 100, 1);

        C = A-1;
        //Matrix& + vec&, row mode
        dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + vec&, row mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, row mode
        dC = dA-B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        C = 1-A;
        //vec& + Matrix&&, row mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& + Matrix&, row mode
        dC = B.toDeviceMatrix()-dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //realM - complexV
    {
        float *dataA = new float[5000];
        float *dataC = new float[10000];
        float *dataD = new float[10000];
        float *dataB = new float[100];
        for (size_t i = 0; i < 10000; i++)
        {
            if (i<5000)dataA[i] = 8;
            dataC[i] = -1;
            dataD[i] = i%2==0?8:0;
            if (i<100)dataB[i] = -1;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto E = Aurora::Matrix::fromRawData(dataD, 50,100,1,Aurora::Complex);



        auto C = E-D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C= D-E;
        //Matrix& + Matrix&&, col mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        dA.forceReshape(100, 50, 1);
        A.forceReshape(100, 50, 1);
        B.forceReshape(1, 50, 1);
        C= E-D;
        //Matrix& + Matrix&, row mode
        dC = (dA-dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()-dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C= D-E;
        //Matrix& + Matrix&&, row mode
        dC = dB-A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixDiv) {
    //real
    {
        auto A = Aurora::zeros(1000,1,1);
        auto B = Aurora::zeros(1000,1,1);
        for (size_t i = 0; i < 1000; i++)
        {
            A[i] = -2.0;
            B[i] =  i;
        }
        printf("Test CudaMatrix operator/(const CudaMatrix &aMatrix) const \r\n");
        //CudaMatrix operator+(const CudaMatrix &aMatrix) const
        auto C = A/B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        auto dC = (dA/dB);
        auto dhC = dC.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dhC[i]);
        }
        printf("Test CudaMatrix operator/(float aScalar) const \r\n");
        //CudaMatrix operator+(float aScalar) const
        auto D = C-0.5;
        auto dD = dC-0.5;
        auto dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix operator/(float aScalar, const CudaMatrix &aMatrix) \r\n");
        // CudaMatrix operator+(float aScalar, const CudaMatrix &aMatrix)
        D = 0.5 / C;
        dD = 0.5 / C.toDeviceMatrix();
        dhD = dD.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(D[i],dhD[i]);
        }
        printf("Test CudaMatrix &operator/(float aScalar, CudaMatrix &&aMatrix) \r\n");
        // CudaMatrix &operator+(float aScalar, CudaMatrix &&aMatrix)
        {
            auto dD2 = 0.5 / (dA/dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        D =  C / 0.5;
        printf("Test CudaMatrix &operator/(CudaMatrix &&aMatrix, float aScalar) \r\n");
        // CudaMatrix &operator+(CudaMatrix &&aMatrix, float aScalar)
        {
            
            auto dD2 =  (dA/dB)/0.5;
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix) const
        printf("Test CudaMatrix operator-(CudaMatrix &&aMatrix) const \r\n");
        {
            auto D = A/C;
            auto dD2 = dA/(dA/dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
        //CudaMatrix operator+(CudaMatrix &&aMatrix,CudaMatrix &aOther)
        printf("Test CudaMatrix operator-(CudaMatrix &&aMatrix,CudaMatrix &aOther) \r\n");
        {
            auto D = A/C;
            auto dD2 = dA/ (dA/dB);
            dhD = dD2.toHostMatrix();
            for (size_t i = 0; i < 1000; i++)
            {
                EXPECT_FLOAT_EQ(D[i],dhD[i]);
            }
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataB[i] = -2;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,100,1,Aurora::Complex);
        auto C = A/B;
        std::complex<float> scalar(-2,-2);
        auto dA = A.toDeviceMatrix();

        auto dC = (dA/scalar);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = A.toDeviceMatrix()/scalar;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B/A;
        dC = scalar/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = A/2.0;
        dC = (dA/2.0);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index"<<i;
        }

        dC = A.toDeviceMatrix()/2.0;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 2.0/A;
        dC = 2.0/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixDivmm) {
    //real
    {
        auto A = Aurora::zeros(100,100,1);
        auto B = Aurora::zeros(100,100,1);
        for (size_t i = 0; i < 10000; i++)
        {
            A[i] = 2;
            B[i] =  8;
        }
        auto C = B/2.0;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();

        auto dC = (dB/dA);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = dB / A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        auto dataA = new float[10000];
        auto dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 2;
            dataB[i] = 8;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,50,1,Aurora::Complex);

        auto C = A/B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& / Matrix&
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& / Matrix&
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B/A;
        //Matrix& / Matrix&&,
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //real & complex
    {
        auto dataA = new float[10000];
        auto dataB = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 2;
            dataB[i] = 8;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 100,50,1,Aurora::Complex);

        auto C = A/B;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix&(real) / Matrix&(complex)
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i<<", A:"<<A[i]<<", B:"<<B[i/2];
        }
        //Matrix&&(real) / Matrix&(complex)
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&(real) / Matrix&&(complex)
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = B/A;
        //Matrix&(complex) / Matrix&&(real)
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&&(complex) / Matrix&(real)
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&(complex) / Matrix&(real)
        dC = dB/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixDivmv) {
    //real
    {
        auto A = Aurora::ones(100,100,1)*8;
        auto B = Aurora::ones(100,1,1)*2;
        auto C = A/2;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 2/A;
        //Matrix& + Matrix&&, col mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        C = A/2;
        //Matrix& - Matrix&, row mode
        dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& - Matrix&, row mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& - vec&&, row mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = 2/ A;
        //Matrix& + Matrix&&, row mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //vec&& + Matrix&&, row mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complex
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            dataC[i] = -2;
            if (i<100)dataB[i] = -2;
        }
        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto C = A/D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, col mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C = D/A;
        //Matrix& + Matrix&&, col mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //vec&& + Matrix&&, col mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        B.forceReshape(1, 50, 1);
        A.forceReshape(100, 50, 1);
        D.forceReshape(100, 50, 1);
        dA.forceReshape(100, 50, 1);

        C= A/D;
        //Matrix& / Matrix&, row mode
        dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& / Matrix&, row mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix& / vec&&, row mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        C= D / A;
        //Matrix& / Matrix&&, row mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& / Matrix&&, row mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //complexM - realv
    {
        float *dataA = new float[10000];
        float *dataB = new float[100];
        float *dataC = new float[10000];
        for (size_t i = 0; i < 10000; i++)
        {
            dataA[i] = 8;
            if (i<100)dataB[i] = 2;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 100,50,1,Aurora::Complex);

        auto B = Aurora::Matrix::fromRawData(dataB, 100,1,1);

        auto C = A/2;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& / Matrix&, col mode
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix&& / vec&, col mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& / vec&&, col mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C=2/A;
        //vec& + Matrix&&, col mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& + Matrix&, col mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 100, 1);
        dA.forceReshape(50, 100, 1);
        A.forceReshape(50, 100, 1);
        B.forceReshape(1, 100, 1);

        C = A/2;
        //Matrix& + vec&, row mode
        dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + vec&, row mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        //Matrix& + vec&&, row mode
        dC = dA/B.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        C = 2/A;
        //vec& + Matrix&&, row mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }

        //vec&& + Matrix&, row mode
        dC = B.toDeviceMatrix()/dA;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
    }
    //realM - complexV
    {
        float *dataA = new float[5000];
        float *dataC = new float[10000];
        float *dataD = new float[10000];
        float *dataB = new float[100];
        for (size_t i = 0; i < 10000; i++)
        {
            if (i<5000)dataA[i] = 8;
            dataC[i] = -2;
            dataD[i] = i%2==0?8:0;
            if (i<100)dataB[i] = -2;
        }

        auto A = Aurora::Matrix::fromRawData(dataA, 50,100,1);
        auto B = Aurora::Matrix::fromRawData(dataB, 50,1,1,Aurora::Complex);
        auto D = Aurora::Matrix::fromRawData(dataC, 50,100,1,Aurora::Complex);
        auto E = Aurora::Matrix::fromRawData(dataD, 50,100,1,Aurora::Complex);



        auto C = E/D;
        auto dA = A.toDeviceMatrix();
        auto dB = B.toDeviceMatrix();
        //Matrix& + Matrix&, col mode
        auto dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
        //Matrix&& + Matrix&, col mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C= D/E;
        //Matrix& + Matrix&&, col mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        dB.forceReshape(1, 50, 1);
        dA.forceReshape(100, 50, 1);
        A.forceReshape(100, 50, 1);
        B.forceReshape(1, 50, 1);
        C= E/D;
        //Matrix& + Matrix&, row mode
        dC = (dA/dB);
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }

        //Matrix&& + Matrix&, row mode
        dC = A.toDeviceMatrix()/dB;
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i));
        }
        C= D/E;
        //Matrix& + Matrix&&, row mode
        dC = dB/A.toDeviceMatrix();
        for (size_t i = 0; i < C.getDataSize()*2; i++)
        {
            EXPECT_FLOAT_EQ(C[i],dC.getValue(i))<<"index:"<<i;
        }
    }
}


TEST_F(CudaMatrix_Test, MatrixComplexMulAndDiv){
    float * dataA = new float[100];
    float * dataB = new float[100];
    auto A = Aurora::Matrix::fromRawData(dataA,5,10,1,Aurora::Complex);
    auto B = Aurora::Matrix::fromRawData(dataB,5,10,1,Aurora::Complex);
    for (size_t i = 0; i < 100; i++)
    {
        dataA[i] = (float)1;
        dataB[i] = (float)(-2);
    }
    auto dA = A.toDeviceMatrix();
    auto dB = B.toDeviceMatrix();
    auto ret1 = A*B;
    auto ret2 = dA*dB;
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }
    ret2 = A.toDeviceMatrix()*dB;
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }
    ret2 = dA*B.toDeviceMatrix();
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }
    ret1 = A/B;
    ret2 = dA/dB;
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }
    ret2 = A.toDeviceMatrix()/dB;
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }
    ret2 = dA/B.toDeviceMatrix();
    for (size_t i = 0; i < ret1.getDataSize()*2; i++)
    {
        EXPECT_FLOAT_EQ(ret1[i], ret2.getValue(i))<<"index:"<<i;
    }

}

TEST_F(CudaMatrix_Test, MatrixPow){
    auto A = Aurora::zeros(1000,1,1);
    auto B = Aurora::zeros(1000,1,1);
    for (size_t i = 0; i < 1000; i++)
    {
        A[i] = -1+0.2*i;
    }
    auto dA= A.toDeviceMatrix();
    {
        auto R = A^0;
        auto dR = dA^0;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R = A^1;
        auto dR = dA^1;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
        {
        auto R = A^2;
        auto dR = dA^2;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R = A^3;
        auto dR = dA^3;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R = A^5;
        auto dR = dA^5;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }

}

TEST_F(CudaMatrix_Test, MatrixNeg){
    auto A = Aurora::zeros(1000,1,1);
    for (size_t i = 0; i < 1000; i++)
    {
        A[i] = -1+0.2*i;
    }
    auto dA= A.toDeviceMatrix();
    {
        auto R = -A;
        auto dR = -dA;
        auto dhR = dR.toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
}

TEST_F(CudaMatrix_Test, MatrixCompare){
    auto A = Aurora::zeros(1000,1,1);
    auto B = Aurora::zeros(1000,1,1);
    for (size_t i = 0; i < 1000; i++)
    {
        A[i] = -1+0.2*i;
    }
    auto dA= A.toDeviceMatrix();
    auto dB= B.toDeviceMatrix();
    {
        auto R= (A<B);
        auto dhR = (dA<dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (A>B);
        auto dhR = (dA>dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }    
    {
        auto R= (A<=B);
        auto dhR = (dA<=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (A>=B);
        auto dhR = (dA>=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (A==B);
        auto dhR = (dA==dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (A!=B);
        auto dhR = (dA!=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9<B);
        auto dhR = (9<dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9>B);
        auto dhR = (9>dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9<=B);
        auto dhR = (9<=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9>=B);
        auto dhR = (9>=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9==B);
        auto dhR = (9==dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
    {
        auto R= (9!=B);
        auto dhR = (9!=dB).toHostMatrix();
        for (size_t i = 0; i < 1000; i++)
        {
            EXPECT_FLOAT_EQ(R[i],dhR[i]);
        }
    }
}


TEST_F(CudaMatrix_Test, matrixfunction)
{
    printf("Test CudaMatrix block(int aDim,int aBeginIndx, int aEndIndex) const\r\n");
    {
        float *dataA = new float[9]{1,6,9,4,1,0,5,8,1};
        float *dataB = new float[3]{1,2,3};
        Aurora::CudaMatrix A = Aurora::Matrix::fromRawData(dataA, 3, 3).toDeviceMatrix();
        Aurora::CudaMatrix B = Aurora::Matrix::fromRawData(dataB, 1, 3).toDeviceMatrix();

        Aurora::Matrix block1 = A.block(0, 1, 2).toHostMatrix();
        Aurora::Matrix block2 = A.toHostMatrix().block(0, 1, 2);
        for (size_t i = 0; i < block1.getDataSize(); i++)
        {   
            EXPECT_FLOAT_EQ(block1[i], block2[i]);
        }

        block2 = A.toHostMatrix();
        A.setBlockValue(0, 1, 2,-1);
        block1 = A.toHostMatrix();
        block2.setBlockValue(0, 1, 2,-1);
        for (size_t i = 0; i < block1.getDataSize(); i++)
        {   
            EXPECT_FLOAT_EQ(block1[i], block2[i]);
        }

        Aurora::CudaMatrix C = Aurora::zeros(2,3).toDeviceMatrix();
        EXPECT_TRUE(A.setBlock(0, 0, 1, C));
        block1 = A.block(0, 0, 1).toHostMatrix();
        block2 = C.toHostMatrix();
        for(size_t i = 0; i < C.getDataSize(); i++)
        {
            EXPECT_FLOAT_EQ(block1[i], block2[i]);
        }       
    }

}