Aurora/src/Matrix.cpp

#include "Matrix.h"

#include <string>
#include <cstring>
#include <iostream>
#include <complex>
//必须在mkl.h和Eigen的头之前，<complex>之后
#define MKL_Complex16 std::complex<double>
#include "mkl.h"

#include "Function.h"
namespace Aurora{
    typedef void(*CalcFuncD)(const MKL_INT n, const double a[], const MKL_INT inca, const double b[],
                             const MKL_INT incb, double r[], const MKL_INT incr);

    typedef void(*CalcFuncZ)(const MKL_INT n, const MKL_Complex16 a[], const MKL_INT inca, const MKL_Complex16 b[],
                             const MKL_INT incb, MKL_Complex16 r[], const MKL_INT incr);

    inline Matrix operatorMxA(CalcFuncD aFunc, double aScalar, const Matrix &aMatrix) {
        double *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
        aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, output, 1);
        if (aMatrix.getValueType() == Complex) {
            aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0, output + 1,
                  1);
        }
        return Matrix::New(output, aMatrix.getDimSize(0), aMatrix.getDimSize(1), aMatrix.getDimSize(2),
                                   aMatrix.getValueType());
    }

    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,
                  aMatrix.getData(),
                  1);
            aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 1, &aScalar, 0,
                  aMatrix.getData() + 1, 1);
        } else {
            aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0,
                  aMatrix.getData(),
                  1);
        }
        std::cout << "after operation" << std::endl;
        aMatrix.printf();
        return aMatrix;
    }


    inline void V_MxM_CN_Calc(
            CalcFuncD aFuncD,
            const int size, double* xC,double* yN,double *output, int DimsStride) {
        aFuncD(size, xC, DimsStride * 2, yN, 1, output, 2);
        aFuncD(size, xC + 1, DimsStride * 2, yN, 1, output + 1, 2);
    }

    inline double* _MxM_CN_Calc(
            CalcFuncD aFuncD,
            const int size, double* xC,double* yN, int dimsStride)
    {
        double *output = malloc(size, true);
        V_MxM_CN_Calc(aFuncD, size, xC, yN, output, dimsStride);
        return output;
    }

    inline void V_MxM_NC_Calc(
            CalcFuncD aFuncD,
            const int size, double* xC,double* yN,double *output, int DimsStride) {
        aFuncD(size, xC, DimsStride, yN, 2, output, 2);
        aFuncD(size, xC , DimsStride, yN+ 1, 2, output + 1, 2);
    }

    inline double* _MxM_NC_Calc(
            CalcFuncD aFuncD,
            const int size, double* xN,double* yC, int dimsStride)
    {
        double *output = malloc(size, true);
        V_MxM_NC_Calc(aFuncD, size, xN, yC, output, dimsStride);
        return output;
    }

    inline void V_MxM_NN_Calc(
            CalcFuncD aFuncD,
            const int size, double* x,double* y,double* output, int DimsStride) {
        aFuncD(size, x, DimsStride, y, 1, output,1);
    }

    inline double* _MxM_NN_Calc(
            CalcFuncD aFuncD,
            const int size, double* x,double* y, int DimsStride) {
        double *output = malloc(size);
        V_MxM_NN_Calc(aFuncD, size, x, y, output, DimsStride);
        return output;
    }

    inline void V_MxM_CC_Calc(
            CalcFuncZ aFuncZ, const int size, double* x,double* y,double* output,
            int DimsStride) {
        aFuncZ(size, (std::complex<double> *) x, DimsStride,
               (std::complex<double> *) y, 1, (std::complex<double> *) output, 1);
    }

    inline double* _MxM_CC_Calc(
            CalcFuncZ aFuncZ, const int size, double* x,double* y,
            int DimsStride) {
        double *output = malloc(size, true);
        V_MxM_CC_Calc(aFuncZ, size, x, y, output, DimsStride);
        return output;
    }

    inline Matrix operatorMxM(CalcFuncD aFuncD, CalcFuncZ aFuncZ, const Matrix &aMatrix,
                              const Matrix &aOther) {
        // 2v2,1v1,3v3
        if (aMatrix.compareShape(aOther)) {
            int DimsStride = 1;
            double *data = nullptr;
            if (aMatrix.getValueType() != aOther.getValueType()) {
                if (aMatrix.getValueType() == Normal) {
                    data = _MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
                                        DimsStride);
                    return Matrix::New(data, aOther);
                } else {
                    data = _MxM_CN_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
                                        DimsStride);
                    return Matrix::New(data, aMatrix);
                }

            } else if (aMatrix.getValueType() == Normal) {
                data = _MxM_NN_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), DimsStride);
                return Matrix::New(data, aMatrix);
            } else {
                data = _MxM_CC_Calc(aFuncZ, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), DimsStride);
                return Matrix::New(data, aMatrix);
            }
        }
        //0v3, 0v2
        else if (aMatrix.getDataSize()==1){
            if (aMatrix.getValueType() ==Normal)return operatorMxA(aFuncD,aMatrix.getData()[0],aOther);
            else{
                std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
                return Matrix();
            }
        }
        //3v0, 2v0
        else if (aOther.getDataSize()==1){
            if (aOther.getValueType() ==Normal)return operatorMxA(aFuncD,aOther.getData()[0],aMatrix);
            else{
                std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
                return Matrix();
            }
        }
        //other
        else {
            std::cerr<<"M * M Shape error!If you want do vector * Matrix, use repmat to replicate the vector"<<std::endl;
            return Matrix();
        }
    }

    inline Matrix operatorMxM_RR(CalcFuncD aFuncD, CalcFuncZ aFuncZ, const Aurora::Matrix &aMatrix,
                                 Aurora::Matrix &&aOther) {
        std::cout << "use right ref operation m" << std::endl;
        if (aMatrix.compareShape(aOther)) {
            int DimsStride = 1;
            if (aMatrix.getValueType() != aOther.getValueType()) {
                //aOther is not a complex matrix
                if (aMatrix.getValueType() == Complex) {
                    double *output = _MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
                                                  DimsStride);
                    return Matrix::New(output, aOther);
                }
                //aOther is a complex matrix, use aOther as output
                V_MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), aOther.getData(),
                              DimsStride);
                return aOther;
            } else if (aMatrix.getValueType() == Normal) {
                V_MxM_NN_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), aOther.getData(),
                              DimsStride);
                return aOther;
            } else {
                V_MxM_CC_Calc(aFuncZ, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), aOther.getData(),
                              DimsStride);
                return aOther;
            }
        }
        //0v3, 0v2
        else if (aMatrix.getDataSize()==1){
            if (aMatrix.getValueType() ==Normal){
                return operatorMxA(aFuncD,aMatrix.getData()[0],std::forward<Aurora::Matrix &&>(aOther));
            }
            else{
                std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
                return Matrix();
            }
        }
            //3v0, 2v0
        else if (aOther.getDataSize()==1){
            if (aOther.getValueType() ==Normal)return operatorMxA(aFuncD,aOther.getData()[0],aMatrix);
            else{
                std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
                return Matrix();
            }
        }
            //other
        else {
            std::cerr<<"M * M Shape error!If you want do vector * Matrix, use repmat to replicate the vector"<<std::endl;
            return Matrix();
        }
    }
};


namespace Aurora {
    Matrix::Matrix(std::shared_ptr<double> aData, std::vector<int> aInfo, ValueType aValueType)
            : mValueType(aValueType)
            , mData(aData)
            , mInfo(aInfo)
    {
        size_t infoSize = mInfo.size();
        for(;infoSize<3;++infoSize)
        {
            mInfo.push_back(1);
        }
    }

    Matrix::Matrix(const Matrix::MatrixSlice& slice) {
        auto temp = slice.toMatrix();
        this->mData = temp.mData;
        this->mInfo = temp.mInfo;
        this->mValueType = temp.mValueType;
    }

    bool Matrix::isNull() const {
        return !mData || mInfo.empty();
    }

    int Matrix::getDims() const {
        if(mInfo[2] > 1)
        {
            return 3;
        }
        return 2;
    }

    double *Matrix::getData() const {
        return mData.get();
    }

    int Matrix::getDimSize(int aIndex) const {
        if (aIndex >= 0 && aIndex < 3) {
            return mInfo.at(aIndex);
        }
        return 0;
    }

    size_t Matrix::getDataSize() const {
        if (!mData.get())return 0;
        size_t ret = 1;
        for (auto v: mInfo) {
            ret *= v;
        }
        return ret;
    }


    bool Matrix::compareShape(const Matrix &other) const {
        if (mInfo.size() != other.mInfo.size()) {
            // all vector compare length
            if (mInfo.size() + other.mInfo.size() == 3){
                return getDataSize()== other.getDataSize();
            }
            // 2 and 3
            else if (mInfo.size() + other.mInfo.size() == 5){
                return (mInfo.size()==3&& mInfo[2]==1)|| (other.mInfo.size()==3&& other.mInfo[2]==1);
            }
            return false;
        }
        for (int i = 0; i < mInfo.size(); ++i) {
            if (mInfo[i] != other.mInfo[i]) return false;
        }
        return true;
    }

    Matrix Matrix::New(double *data, int rows, int cols, int slices, ValueType type) {
        if (!data) return Matrix();
        std::vector<int> vector;
        vector.push_back(rows);
        if (cols > 0)vector.push_back(cols);
        if (slices > 0)vector.push_back(slices);
        Matrix ret({data, free}, vector);
        if (type != Normal)ret.setValueType(type);
        return ret;
    }

    Matrix Matrix::New(double *data, const Matrix &shapeMatrix) {
        return New(data,
                   shapeMatrix.getDimSize(0),
                   shapeMatrix.getDimSize(1),
                   shapeMatrix.getDimSize(2),
                   shapeMatrix.getValueType());
    }

    Matrix Matrix::New(const Matrix &shapeMatrix) {
        double *newBuffer = malloc(shapeMatrix.getDataSize(), shapeMatrix.getValueType());
        return New(newBuffer, shapeMatrix);
    }

    Matrix Matrix::fromRawData(double *data, int rows, int cols, int slices, ValueType type) {
        if (!data) return Matrix();
        std::vector<int> vector;
        vector.push_back(rows);
        if (cols > 0)vector.push_back(cols);
        if (slices > 0)vector.push_back(slices);
        Matrix ret({data, std::default_delete<double[]>()}, vector);
        if (type != Normal)ret.setValueType(type);
        return ret;
    }

    Matrix Matrix::copyFromRawData(double *data, int rows, int cols, int slices, ValueType type) {
        int colsize = cols>0?cols:1;
        int slicesize = slices>0?slices:1;
        int size = rows*colsize*slicesize;
        double *newBuffer = malloc(size, type);
        cblas_dcopy(size*type,data,1,newBuffer,1);
        return New(newBuffer,rows,cols,slices,type);
    }

    Matrix Matrix::deepCopy() const {
        double *newBuffer = malloc(getDataSize(), getValueType()==Complex);
        cblas_dcopy(getDataSize() * getValueType(),getData(),1,newBuffer,1);
        return New(newBuffer,
                   getDimSize(0),
                   getDimSize(1),
                   getDimSize(2),
                   getValueType());
    }

    //operation +
    Matrix Matrix::operator+(double aScalar) const { return operatorMxA(&vdAddI, aScalar, *this);}
    Matrix operator+(double aScalar, const Matrix &matrix) {return matrix + aScalar;}
    Matrix Matrix::operator+(const Matrix &matrix) const {return operatorMxM(vdAddI, vzAddI, *this, matrix);}
    Matrix &operator+(double aScalar, Matrix &&matrix) {
        return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix &operator+(Matrix &&matrix,double aScalar) {
        return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix Matrix::operator+(Matrix &&aMatrix) const {
        return operatorMxM_RR(&vdAddI,&vzAddI,*this,std::forward<Matrix&&>(aMatrix));
    }
    Matrix operator+(Matrix &&aMatrix, const Matrix &aOther) {
        return operatorMxM_RR(&vdAddI,&vzAddI,aOther,std::forward<Matrix&&>(aMatrix));
    }

    //operation -
    Matrix Matrix::operator-(double aScalar) const { return operatorMxA(&vdSubI, aScalar, *this);}
    Matrix operator-(double aScalar, const Matrix &matrix) {return matrix - aScalar;}
    Matrix Matrix::operator-(const Matrix &matrix) const {return operatorMxM(vdSubI, vzSubI, *this, matrix);}
    Matrix &operator-(double aScalar, Matrix &&matrix) {
        return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix &operator-(Matrix &&matrix,double aScalar) {
        return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix Matrix::operator-(Matrix &&aMatrix) const {
        return operatorMxM_RR(&vdSubI,&vzSubI,*this,std::forward<Matrix&&>(aMatrix));
    }
    Matrix operator-(Matrix &&aMatrix, const Matrix &aOther) {
        return operatorMxM_RR(&vdSubI,&vzSubI,aOther,std::forward<Matrix&&>(aMatrix));
    }

    //operation *
    Matrix Matrix::operator*(double aScalar) const { return operatorMxA(&vdMulI, aScalar, *this);}
    Matrix operator*(double aScalar, const Matrix &matrix) {return matrix * aScalar;}
    Matrix Matrix::operator*(const Matrix &matrix) const {return operatorMxM(vdMulI, vzMulI, *this, matrix);}
    Matrix &operator*(double aScalar, Matrix &&matrix) {
        return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix &operator*(Matrix &&matrix,double aScalar) {
        return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix Matrix::operator*(Matrix &&aMatrix) const {
        return operatorMxM_RR(&vdMulI,&vzMulI,*this,std::forward<Matrix&&>(aMatrix));
    }
    Matrix operator*(Matrix &&aMatrix, const Matrix &aOther) {
        return operatorMxM_RR(&vdMulI,&vzMulI,aOther,std::forward<Matrix&&>(aMatrix));
    }

    //operation /
    Matrix Matrix::operator/(double aScalar) const { return operatorMxA(&vdDivI, aScalar, *this);}
    Matrix operator/(double aScalar, const Matrix &matrix) {return matrix / aScalar;}
    Matrix Matrix::operator/(const Matrix &matrix) const {return operatorMxM(vdDivI, vzDivI, *this, matrix);}
    Matrix &operator/(double aScalar, Matrix &&matrix) {
        return operatorMxA_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix &operator/(Matrix &&matrix,double aScalar) {
        return operatorMxA_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
    }
    Matrix Matrix::operator/(Matrix &&aMatrix) const {
        return operatorMxM_RR(&vdDivI,&vzDivI,*this,std::forward<Matrix&&>(aMatrix));
    }
    Matrix operator/(Matrix &&aMatrix, const Matrix &aOther) {
        return operatorMxM_RR(&vdDivI,&vzDivI,aOther,std::forward<Matrix&&>(aMatrix));
    }

    //operator ^ (pow)
    Matrix Matrix::operator^(int times) const { return operatorMxA(&vdPowI, times, *this);}
    Matrix operator^( Matrix &&matrix,int times) {
        return operatorMxA(&vdPowI, times, std::forward<Matrix&&>(matrix));
    }
    



    void Matrix::printf() {
        int k_count = getDimSize(2);
        int j_count = getDimSize(1);
        int complexstep = 1;
        const char* mark = "+";
        if(mValueType == Complex)
        {
            complexstep = 2;
        }
        for (int k = 0; k <k_count; ++k) {
        ::printf("slice %d:\r\n[",k);
        for (int i = 0; i < getDimSize(0); ++i) {
            ::printf("[");
            for (int j = 0; j < j_count; ++j) {
                ::printf("%f2",getData()[(k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i)*complexstep]);
                if(mValueType == Complex)
                {
                    double value = getData()[(k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i)*complexstep+1];
                    if(value<0)
                    {
                        mark = "-";
                    }
                    ::printf("%s%f2i",mark,value);
                }
                ::printf(", ");
            }
            ::printf("]\r\n");
        }
        ::printf("]\r\n");
        }
    }

    Matrix::MatrixSlice Matrix::operator()(int aRowIdx, int aColIdx, int aSliceIdx) const {
        std::vector<int> dims = {aRowIdx, aColIdx, aSliceIdx};
        std::vector<int> allDimIndex;
        int mode = 0;
        for (int j = 0; j < 3; ++j) {
            if (dims[j]==$ && this->getDims()>j){
                ++mode;
                allDimIndex.push_back(j);
            }
        }
        int rowStride = 1;
        int colStride =  getDimSize(0);
        int sliceStride =  getDimSize(0)*getDimSize(1);
        int strides[3] = {rowStride, colStride, sliceStride};
        int rowOffset = aRowIdx == $ ? 0 : aRowIdx;
        int colOffset = aColIdx == $ ? 0 : aColIdx;
        int sliceOffset = aSliceIdx == $ ? 0 : aSliceIdx;
        double *startPointer = getData() + (rowStride * rowOffset
                                            + colStride * colOffset
                                            + sliceStride * sliceOffset) * getValueType();
        int size1 = allDimIndex.empty()?1:getDimSize(allDimIndex[0]);
        int stride1 = allDimIndex.empty()?1:strides[allDimIndex[0]];
        switch (mode) {
            //matrix mode
            case 2:{
                int size2 = getDimSize(allDimIndex[1]);
                int stride2 = strides[allDimIndex[1]];
                return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode, size2, stride2);
            }
            //vector mode & default
            case 1:
            {
                return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode);
            }
            //scalar mode or ALL $
            case 0:
            default: {
                return Matrix::MatrixSlice(1 , 1, startPointer,getValueType(), mode);
            }
        }
    }

    Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2):
    mSliceMode(aSliceMode),mData(aData),
    mSize(aSize),mSize2(aSize2),
    mStride(aStride),mStride2(aStride2),
    mType(aType)
    {

    }

    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix::MatrixSlice &slice) {
        if (this==&slice) return *this;
        if(!mData){
            std::cerr <<"Assign value fail!Des data pointer is null!";
            return *this;
        }
        if(!slice.mData){
            std::cerr <<"Assign value fail!Src data pointer is null!"<<std::endl;
            return *this;
        }
        if (slice.mSliceMode!=mSliceMode) {
            std::cerr << "Assign value fail!Src slice(dims count:" << slice.mSliceMode
                      << "), not match of des(dims count:" << mSliceMode << ")!" << std::endl;
            return *this;
        }
        if (slice.mSize!=mSize) {
            std::cerr <<"Assign value fail!Src slice(dim 1 size:"<< slice.mSize <<"), not match of des(dim 1 size:"<<mSize<<")!"<<std::endl;
            return *this;
        }
        if (slice.mSize2!=mSize2) {
            std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< slice.mSize2 <<"), not match of des(dim 2 size:"<<mSize2<<")!"<<std::endl;
            return *this;
        }
        if (slice.mType!=mType) {
            std::cerr <<"Assign value fail!Src slice(value type:"<< slice.mType <<"), not match of des(value type:"<<mType<<")!"<<std::endl;
            return *this;
        }
        switch (mSliceMode) {
            case 2:{
                if (mType== Normal) {
                    cblas_dcopy_batch_strided(mSize, slice.mData, slice.mStride, slice.mStride2, mData, mStride,
                                              mStride2, mSize2);
                }
                else {
                    cblas_zcopy_batch_strided(mSize,(std::complex<double>*)slice.mData,slice.mStride,slice.mStride2,
                            (std::complex<double>*)mData,mStride,mStride2,mSize2);
                }
                break;
            }
            case 1:{
                if (mType== Normal){
                    cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride);
                }
                else {
                    cblas_zcopy(mSize, (std::complex<double> *) slice.mData, slice.mStride,
                                (std::complex<double> *) mData, mStride);
                }
                break;
            }
            case 0:
            default:{
                mData[0] = slice.mData[0];
                if (mType != Normal)mData[1] = slice.mData[1];
            }
        }
        return *this;
    }

    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix &matrix) {
        if(!mData){
            std::cerr <<"Assign value fail!Des data pointer is null!";
            return *this;
        }
        if(!matrix.getData()){
            std::cerr <<"Assign value fail!Src data pointer is null!";
            return *this;
        }
        if (matrix.getDims()!=mSliceMode) {
            std::cerr <<"Assign value fail!Src matrix(dims count:"<< matrix.getDims() <<"), not match of des(dims count:"<<mSliceMode<<")!";
            return *this;
        }
        if (matrix.getDimSize(0)!=mSize) {
            std::cerr <<"Assign value fail!Src matrix(dim 1 size:"<< matrix.getDimSize(0)<<"), not match of des(dim 1 size:"<<mSize<<")!";
            return *this;
        }
        if (matrix.getDimSize(1)!=mSize2) {
            std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< matrix.getDimSize(1) <<"), not match of des(dim 2 size:"<<mSize2<<")!";
            return *this;
        }
        if (matrix.getValueType()!=mType) {
            std::cerr <<"Assign value fail!Src slice(value type:"<< matrix.getValueType() <<"), not match of des(value type:"<<mType<<")!";
            return *this;
        }
        switch (mSliceMode) {
            case 2:{
                if (mType== Normal) {
                    cblas_dcopy_batch_strided(mSize, matrix.getData(), 1, matrix.getDimSize(0), mData, mStride,
                                              mStride2, mSize2);
                }
                else {
                    cblas_zcopy_batch_strided(mSize,(std::complex<double>*)matrix.getData(),1,matrix.getDimSize(0),
                                              (std::complex<double>*)mData,mStride,mStride2,mSize2);
                }
                break;
            }
            case 1:{
                if (mType== Normal){
                    cblas_dcopy(mSize,matrix.getData(),1,mData,mStride);
                }
                else {
                    cblas_zcopy(mSize, (std::complex<double> *) matrix.getData(),1,
                                (std::complex<double> *) mData, mStride);
                }
                break;
            }
            case 0:
            default:{
                mData[0] = matrix.getData()[0];
                if (mType != Normal)mData[1] = matrix.getData()[1];
            }
        }
        return *this;
    }

    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(double value) {
        if(!mData){
            std::cerr <<"Assign value fail!Des data pointer is null!";
            return *this;
        }
        if (mSliceMode!=0) {
            std::cerr <<"Assign value fail!Des slicemode is"<<mSliceMode<<", not scalar mode!";
            return *this;
        }
        if (mSize!=1) {
            std::cerr <<"Assign value fail!Des  size:"<<mSize<<", not scalar mode!";
            return *this;
        }
        if (mType!=Normal) {
            std::cerr <<"Assign value fail!Des type is complex!";
            return *this;
        }
        mData[0]=value;
        return *this;
    }

    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(std::complex<double> value) {
        if(!mData){
            std::cerr <<"Assign value fail!Des data pointer is null!";
            return *this;
        }
        if (mSliceMode!=0) {
            std::cerr <<"Assign value fail!Des slicemode is"<<mSliceMode<<", not scalar mode!";
            return *this;
        }
        if (mSize!=1) {
            std::cerr <<"Assign value fail!Des  size:"<<mSize<<", not scalar mode!";
            return *this;
        }
        if (mType!=Complex) {
            std::cerr <<"Assign value fail!Des type is not complex!";
            return *this;
        }
        mData[0]=value.real();
        mData[1]=value.imag();
        return *this;
    }

    Matrix Matrix::MatrixSlice::toMatrix() const {
        double * data = (double *) mkl_malloc(mSize*(mSize2>0?mSize2:1) * sizeof(double)*mType, 64);

        switch (mSliceMode) {
            case 2:{
                if (mType== Normal) {
                    cblas_dcopy_batch_strided(mSize,  mData, mStride,
                                              mStride2,data, 1, mSize, mSize2);
                }
                else {
                    cblas_zcopy_batch_strided(mSize, (std::complex<double> *) mData, mStride, mStride2,
                                              (std::complex<double> *) data, 1, mSize,
                                              mSize2);
                }
                break;
            }
            case 1:{
                if (mType== Normal){
                    cblas_dcopy(mSize,mData,mStride,data,1);
                }
                else {
                    cblas_zcopy(mSize, (std::complex<double> *) mData, mStride,
                                (std::complex<double> *) data, 1);
                }
                break;
            }
            case 0:
            default:{
                data[0]= mData[0];
                if (mType != Normal) data[1] = mData[1];
            }
        }

        return Matrix::New(data,mSize,mSize2,0,mType);
    }
}