From 8167775d38fba69cd2cd3a7454d0bfb6b871ed4d Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Sun, 8 Oct 2023 15:58:43 +0800
Subject: [PATCH] Make aurora from Double to float

---
 src/AuroraDefs.h     |   2 +-
 src/Function.cpp     | 151 +-------------
 src/Function.h       |  15 +-
 src/Function1D.cpp   | 316 ++++++++++++++---------------
 src/Function1D.h     |  20 +-
 src/Function2D.cpp   | 310 ++++++++++++++--------------
 src/Function2D.h     |   4 +-
 src/Function3D.cpp   |  42 ++--
 src/Function3D.h     |   2 +-
 src/MatlabReader.cpp |  23 ++-
 src/MatlabWriter.cpp |   2 +-
 src/Matrix.cpp       | 466 +++++++++++++++++++++----------------------
 src/Matrix.h         | 100 +++++-----
 13 files changed, 650 insertions(+), 803 deletions(-)

diff --git a/src/AuroraDefs.h b/src/AuroraDefs.h
index a20d163..53fa6a1 100644
--- a/src/AuroraDefs.h
+++ b/src/AuroraDefs.h
@@ -4,7 +4,7 @@
 #define EIGEN_USE_MKL_ALL
 #include <complex>
 //必须在mkl.h和Eigen的头之前，<complex>之后
-#define MKL_Complex16 std::complex<double>
+#define MKL_Complex8 std::complex<float>
 #include "mkl.h"
 
 #define PI 3.141592653589793238462
diff --git a/src/Function.cpp b/src/Function.cpp
index 59ee631..6000ba4 100644
--- a/src/Function.cpp
+++ b/src/Function.cpp
@@ -5,7 +5,7 @@
 #include "Function.h"
 
 //必须在mkl.h和Eigen的头之前，<complex>之后
-#define MKL_Complex16 std::complex<double>
+#define MKL_Complex8 std::complex<float>
 #include "mkl.h"
 #include <string.h>
 #include <algorithm>
@@ -15,155 +15,14 @@
 #include <Eigen/Dense>
 
 namespace Aurora {
-    double immse(double *dataA, double *dataB, int size) {
-        auto temp = new double[size];
-        vdSub(size, dataA, dataB, temp);
-        vdSqr(size, temp, temp);
-        double result = cblas_dasum(size, temp, 1) / (double) size;
-        delete[] temp;
-        return result;
-    }
 
-    double *std(int rows, int cols, double *input) {
-        auto std = new double[cols];
-//    #pragma omp parallel for num_threads(4)
-        for (int i = 0; i < cols; ++i) {
-            double *p = input + i * rows;
-            double mean = cblas_dasum(rows, p, 1) / rows;
-            vdSubI(rows, p, 1, &mean, 0, p, 1);
-            vdSqr(rows, p, p);
-            std[i] = cblas_dasum(rows, p, 1) / (rows - 1);
-        }
-        vdSqrt(cols, std, std);
-        return std;
-    }
-
-    double *inv(int cols, double *pMatrix) {
-        int size = cols * cols;
-        int *ipiv = new int[cols];
-        auto result = new double[size];
-        memcpy(result, pMatrix, sizeof(double) * size);
-        LAPACKE_dgetrf(LAPACK_ROW_MAJOR, cols, cols, result, cols, ipiv);
-        LAPACKE_dgetri(LAPACK_ROW_MAJOR, cols, result, cols, ipiv);
-        delete[] ipiv;
-        return result;
-    }
-
-    std::complex<double> *fft(long int size, std::complex<double> *input) {
-        DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
-
-        auto output = new std::complex<double>[size];
-        //使用ce数据fft
-        MKL_LONG status;
-
-        //创建 Descriptor, 精度 double , 输入类型实数, 维度1
-        status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, size);
-        //通过 setValue 配置Descriptor
-        //使用单独的输出数据缓存
-        status = DftiSetValue(my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
-
-        //提交 修改配置后的Descriptor(实际上会进行FFT的计算初始化)
-        status = DftiCommitDescriptor(my_desc_handle);
-        //执行计算
-        DftiComputeForward(my_desc_handle, input, output);
-        //释放资源
-        status = DftiFreeDescriptor(&my_desc_handle);
-        return output;
-    }
-
-    std::complex<double> *ifft(long int size, std::complex<double> *input) {
-        DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
-        auto output = new std::complex<double>[size];
-        MKL_LONG status;
-
-        //创建 Descriptor, 精度 double , 输入类型实数, 维度1
-        status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, size);
-        //通过 setValue 配置Descriptor
-        //使用单独的输出数据缓存
-        status = DftiSetValue(my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
-        //设置DFTI_BACKWARD_SCALE ！！！很关键，不然值不对
-        status = DftiSetValue(my_desc_handle, DFTI_BACKWARD_SCALE, 1.0f / size);
-        status = DftiCommitDescriptor(my_desc_handle);
-        //提交 修改配置后的Descriptor(实际上会进行FFT的计算初始化)
-        status = DftiCommitDescriptor(my_desc_handle);
-        //执行计算
-        DftiComputeBackward(my_desc_handle, input, output);
-        //释放资源
-        status = DftiFreeDescriptor(&my_desc_handle);
-        return output;
-    }
-
-    double *real(int size, std::complex<double> *input) {
-        auto input_d = (double *) input;
-        auto output = new double[size];
-        const int complex_stride = 2;
-        const int real_stride = 1;
-        cblas_dcopy(size, input_d, complex_stride, output, real_stride);
-        return output;
-    }
-
-    std::complex<double> *complex(int size, double *input) {
-        auto output = (std::complex<double> *) mkl_malloc(size * sizeof(std::complex<double>), 64);
-        memset(output, 0, size * sizeof(std::complex<double>));
-        const int complex_stride = 2;
-        const int real_stride = 1;
-        cblas_dcopy(size, input, real_stride, (double *) output, complex_stride);
-        return output;
-    }
-
-    std::complex<double> *hilbert(int size, double *input) {
-        auto complexInput = complex(size, input);
-        auto x = fft(size, complexInput);
-        auto h = new double[size];
-        auto two = 2.0;
-        auto zero = 0.0;
-        cblas_dcopy(size, &zero, 0, h, 1);
-        cblas_dcopy(size / 2, &two, 0, h, 1);
-        h[size / 2] = ((size << 31) >> 31) ? 2.0 : 1.0;
-        h[0] = 1.0;
-        auto p = (double *) x;
-        vdMulI(size, p, 2, h, 1, p, 2);
-        vdMulI(size, p + 1, 2, h, 1, p + 1, 2);
-        auto result = ifft(size, x);
-        delete[] h;
-        delete[] x;
-        return result;
-    }
-
-    double *malloc(size_t size, bool complex) {
-        if (!complex) return (double *) mkl_malloc(size * sizeof(double), 64);
-        size_t complex_size = size * sizeof(std::complex<double>);
-        return (double *) mkl_malloc(complex_size, 64);
+    float *malloc(size_t size, bool complex) {
+        if (!complex) return (float *) mkl_malloc(size * sizeof(float), 64);
+        size_t complex_size = size * sizeof(std::complex<float>);
+        return (float *) mkl_malloc(complex_size, 64);
     }
 
     void free(void* ptr){
         mkl_free(ptr);
     }
-
-    double * mul(double scalar, double *input, int size) {
-        double* output = malloc(size);
-        vdMulI(size,input,1,&scalar,0,output,1);
-        return output;
-    }
-
-    double *mul(double *inputA, double *inputB, int size) {
-        double* output = malloc(size);
-        vdMulI(size,inputA,1,inputB,1,output,1);
-        return output;
-    }
-
-    double *mulz(std::complex<double> *inputA, std::complex<double> *inputB, int size) {
-        double* output = malloc(size,true);
-        vzMulI(size,inputA,1,inputB,1,(std::complex<double>*)output,1);
-        return output;
-    }
-
-    double *random( int size) {
-        double * data = malloc(size);
-        Eigen::Map<Eigen::VectorXd> srcV(data,size);
-        srcV.setRandom();
-        return data;
-    }
-
-
 }
\ No newline at end of file
diff --git a/src/Function.h b/src/Function.h
index 225741e..4be1195 100644
--- a/src/Function.h
+++ b/src/Function.h
@@ -8,21 +8,8 @@
 #include <complex>
 
 namespace Aurora{
-    double* malloc(size_t size,bool complex = false);
+    float* malloc(size_t size,bool complex = false);
     void free(void* ptr);
-    double* mul( double scalar, double * input, int size);
-    double* mul( double* inputA, double * inputB, int size);
-    double* mulz( std::complex<double> *inputA, std::complex<double> * inputB, int size);
-    double immse(double * dataA, double * dataB, int size);
-    double* std(int rows, int cols, double * input);
-    double* random(int size);
-    double* inv(int cols,double *pMatrix);
-    double* real(int size, std::complex<double> * input);
-    std::complex<double> * complex(int size, double * input);
-    std::complex<double>* fft(long int size, std::complex<double> * input);
-//  ic std::complex<double>* fft(long int size, double * input);
-    std::complex<double>* ifft(long int size, std::complex<double> * input);
-    std::complex<double>* hilbert(int size, double * input);
 };
 
 
diff --git a/src/Function1D.cpp b/src/Function1D.cpp
index 11542a0..4a853a3 100644
--- a/src/Function1D.cpp
+++ b/src/Function1D.cpp
@@ -28,26 +28,26 @@ namespace {
     const int COMPLEX_STRIDE = 2;
     const int REAL_STRIDE = 1;
     const int SAME_STRIDE = 1;
-    const double VALUE_ONE = 1.0;
+    const float VALUE_ONE = 1.0;
     const ushort CONVERT_AND_VALUE = 15;
     const ushort CONVERT_AND_VALUE_2 = 2047;
     const ushort CONVERT_MUL_VALUE = 2048;
 
     uint CONVERT_ADD_VALUE = UINT32_MAX - 4095;
 
-    inline void convertValue(double aValue ,double* des){
-      double value = aValue;
+    inline void convertValue(float aValue ,float* des){
+      float value = aValue;
       ushort *exponentPtr = (ushort *)&value;
       exponentPtr[0] = (exponentPtr[0] >> 11) & CONVERT_AND_VALUE;
       exponentPtr[1] = (exponentPtr[1] >> 11) & CONVERT_AND_VALUE;
       exponentPtr[2] = (exponentPtr[2] >> 11) & CONVERT_AND_VALUE;
       exponentPtr[3] = (exponentPtr[3] >> 11) & CONVERT_AND_VALUE;
-      double signValue = aValue;
+      float signValue = aValue;
       short *signPtr = (short *)&signValue;
       uint sign_bit[4] = {
           (uint)(signPtr[0] < 0 ? 1 : 0), (uint)(signPtr[1] < 0 ? 1 : 0),
           (uint)(signPtr[2] < 0 ? 1 : 0), (uint)(signPtr[3] < 0 ? 1 : 0)};
-      double fraction3Value = aValue;
+      float fraction3Value = aValue;
       ushort *fraction3Ptr = (ushort *)&fraction3Value;
       fraction3Ptr[0] &= CONVERT_AND_VALUE_2;
       fraction3Ptr[1] &= CONVERT_AND_VALUE_2;
@@ -83,7 +83,7 @@ namespace {
 
     }
 
-    inline void convertValue2(short* aValue ,double* des){
+    inline void convertValue2(short* aValue ,float* des){
       ushort exponentPtr[4] = {(ushort)aValue[0],(ushort)aValue[1],(ushort)aValue[2],(ushort)aValue[3]};
       exponentPtr[0] = (exponentPtr[0] >> 11) & CONVERT_AND_VALUE;
       exponentPtr[1] = (exponentPtr[1] >> 11) & CONVERT_AND_VALUE;
@@ -171,9 +171,9 @@ Aurora::Matrix Aurora::complex(const Aurora::Matrix &matrix) {
         return matrix;
     }
     auto output = malloc(matrix.getDataSize() ,true);
-    memset(output, 0, (matrix.getDataSize() * sizeof(std::complex<double>)));
-    cblas_dcopy(matrix.getDataSize(), matrix.getData(), REAL_STRIDE, (double *) output, COMPLEX_STRIDE);
-    return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2),
+    memset(output, 0, (matrix.getDataSize() * sizeof(std::complex<float>)));
+    cblas_scopy(matrix.getDataSize(), matrix.getData(), REAL_STRIDE, (float *) output, COMPLEX_STRIDE);
+    return Aurora::Matrix::New((float *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2),
                                Complex);
 }
 
@@ -182,10 +182,10 @@ Aurora::Matrix Aurora::real(const Aurora::Matrix &matrix) {
         std::cerr<<"real only support complex value type"<<std::endl;
         return matrix;
     }
-    auto output = (double *) malloc(matrix.getDataSize());
-    memset(output, 0, (matrix.getDataSize() * sizeof(double)));
-    cblas_dcopy(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE , (double *) output, REAL_STRIDE);
-    return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
+    auto output = (float *) malloc(matrix.getDataSize());
+    memset(output, 0, (matrix.getDataSize() * sizeof(float)));
+    cblas_scopy(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE , (float *) output, REAL_STRIDE);
+    return Aurora::Matrix::New((float *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
 }
 
 Aurora::Matrix Aurora::imag(const Aurora::Matrix &matrix) {
@@ -194,28 +194,28 @@ Aurora::Matrix Aurora::imag(const Aurora::Matrix &matrix) {
         return matrix;
     }
     auto output = malloc(matrix.getDataSize());
-    memset(output, 0, (matrix.getDataSize() * sizeof(double)));
-    cblas_dcopy(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE , (double *) output, REAL_STRIDE);
-    return Aurora::Matrix::New((double *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
+    memset(output, 0, (matrix.getDataSize() * sizeof(float)));
+    cblas_scopy(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE , (float *) output, REAL_STRIDE);
+    return Aurora::Matrix::New((float *) output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
 }
 
 Aurora::Matrix Aurora::ceil(const Aurora::Matrix &matrix) {
     auto output = malloc(matrix.getDataSize());
     //for real part
-    vdCeilI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
+    vsCeilI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdCeilI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
+        vsCeilI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
     }
     return Aurora::Matrix::New(output, matrix);
 }
 
 Aurora::Matrix Aurora::ceil(const Aurora::Matrix &&matrix) {
     //for real part
-    vdCeilI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
+    vsCeilI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdCeilI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
+        vsCeilI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
     }
     return matrix;
 }
@@ -223,20 +223,20 @@ Aurora::Matrix Aurora::ceil(const Aurora::Matrix &&matrix) {
 Aurora::Matrix Aurora::round(const Aurora::Matrix &matrix) {
     auto output = malloc(matrix.getDataSize());
     //for real part
-    vdRoundI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
+    vsRoundI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdRoundI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
+        vsRoundI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
     }
     return Aurora::Matrix::New(output, matrix);
 }
 
 Aurora::Matrix Aurora::round(const Aurora::Matrix &&matrix) {
     //for real part
-    vdRoundI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
+    vsRoundI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdRoundI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
+        vsRoundI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
     }
     return matrix;
 }
@@ -244,20 +244,20 @@ Aurora::Matrix Aurora::round(const Aurora::Matrix &&matrix) {
 Aurora::Matrix Aurora::floor(const Aurora::Matrix &matrix) {
     auto output = malloc(matrix.getDataSize());
     //for real part
-    vdFloorI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
+    vsFloorI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdFloorI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
+        vsFloorI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, output + 1, SAME_STRIDE);
     }
     return Aurora::Matrix::New(output, matrix);
 }
 
 Aurora::Matrix Aurora::floor(const Aurora::Matrix &&matrix) {
     //for real part
-    vdFloorI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
+    vsFloorI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
     if (matrix.getValueType() == Complex) {
         //for imag part
-        vdFloorI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
+        vsFloorI(matrix.getDataSize(), matrix.getData() + 1, SAME_STRIDE, matrix.getData() + 1, SAME_STRIDE);
     }
     return matrix;
 }
@@ -267,7 +267,7 @@ Matrix Aurora::auroraNot(const Matrix& aMatrix){
 }
 
 Matrix Aurora::auroraNot(Matrix&& aMatrix){
-    Eigen::Map<Eigen::VectorXd> v2(aMatrix.getData(), aMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> v2(aMatrix.getData(), aMatrix.getDataSize());
     v2 = (v2.array()>0).select(1,v2);
     v2 = (v2.array()<0).select(0,v2);
     v2 = v2.array()+1.0;
@@ -279,7 +279,7 @@ Aurora::Matrix Aurora::sqrt(const Aurora::Matrix& matrix) {
 
     if (matrix.getValueType() != Complex) {
         auto output = malloc(matrix.getDataSize());
-        vdSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
+        vsSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
         return Aurora::Matrix::New(output, matrix);
     }
     std::cerr<<"sqrt not support complex"<<std::endl;
@@ -288,7 +288,7 @@ Aurora::Matrix Aurora::sqrt(const Aurora::Matrix& matrix) {
 
 Aurora::Matrix Aurora::sqrt(Aurora::Matrix&& matrix) {
     if (matrix.getValueType() != Complex) {
-        vdSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
+        vsSqrtI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
         return matrix;
     }
     std::cerr<<"sqrt not support complex"<<std::endl;
@@ -298,23 +298,23 @@ Aurora::Matrix Aurora::sqrt(Aurora::Matrix&& matrix) {
 Aurora::Matrix Aurora::abs(const Aurora::Matrix &matrix) {
     auto output = malloc(matrix.getDataSize());
     if (matrix.getValueType()==Normal){
-        vdAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
+        vsAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, output, SAME_STRIDE);
     }
     else{
-        vzAbsI(matrix.getDataSize(), (std::complex<double> *)matrix.getData(), SAME_STRIDE,output, SAME_STRIDE);
+        vcAbsI(matrix.getDataSize(), (std::complex<float> *)matrix.getData(), SAME_STRIDE,output, SAME_STRIDE);
     }
     return Aurora::Matrix::New(output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
 }
 
 Aurora::Matrix Aurora::abs(Aurora::Matrix&& matrix) {
     if (matrix.getValueType()==Normal){
-        vdAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
+        vsAbsI(matrix.getDataSize(), matrix.getData(), SAME_STRIDE, matrix.getData(), SAME_STRIDE);
         return matrix;
     }
     //TODO：考虑尝试是不是使用realloc缩短已分配的内存的方式重用matrix
     else{
         auto output = malloc(matrix.getDataSize());
-        vzAbsI(matrix.getDataSize(), (std::complex<double> *)matrix.getData(), SAME_STRIDE,output, SAME_STRIDE);
+        vcAbsI(matrix.getDataSize(), (std::complex<float> *)matrix.getData(), SAME_STRIDE,output, SAME_STRIDE);
         return Aurora::Matrix::New(output, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
@@ -323,7 +323,7 @@ Aurora::Matrix Aurora::abs(Aurora::Matrix&& matrix) {
 Aurora::Matrix Aurora::sign(const Aurora::Matrix &matrix) {
     if (matrix.getValueType()==Normal){
         auto ret = matrix.deepCopy();
-        Eigen::Map<Eigen::VectorXd> retV(ret.getData(),ret.getDataSize());
+        Eigen::Map<Eigen::VectorXf> retV(ret.getData(),ret.getDataSize());
         retV = retV.array().sign();
         return ret;
     }
@@ -331,9 +331,9 @@ Aurora::Matrix Aurora::sign(const Aurora::Matrix &matrix) {
         //sign(x) = x./abs(x)，前提是 x 为复数。
         auto output = malloc(matrix.getDataSize(),true);
         Matrix absMatrix = abs(matrix);
-        vdDivI(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE,
+        vsDivI(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE,
                 absMatrix.getData(), REAL_STRIDE,output,COMPLEX_STRIDE);
-        vdDivI(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE,
+        vsDivI(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE,
                absMatrix.getData(), REAL_STRIDE,output+1,COMPLEX_STRIDE);
         return Aurora::Matrix::New(output, matrix);
     }
@@ -341,16 +341,16 @@ Aurora::Matrix Aurora::sign(const Aurora::Matrix &matrix) {
 
 Aurora::Matrix Aurora::sign(Aurora::Matrix&& matrix) {
     if (matrix.getValueType()==Normal){
-        Eigen::Map<Eigen::VectorXd> retV(matrix.getData(),matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> retV(matrix.getData(),matrix.getDataSize());
         retV = retV.array().sign();
         return matrix;
     }
     else{
         //sign(x) = x./abs(x)，前提是 x 为复数。
         Matrix absMatrix = abs(matrix);
-        vdDivI(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE,
+        vsDivI(matrix.getDataSize(), matrix.getData(),COMPLEX_STRIDE,
                absMatrix.getData(), REAL_STRIDE,matrix.getData(),COMPLEX_STRIDE);
-        vdDivI(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE,
+        vsDivI(matrix.getDataSize(), matrix.getData()+1,COMPLEX_STRIDE,
                absMatrix.getData(), REAL_STRIDE,matrix.getData()+1,COMPLEX_STRIDE);
         return matrix;
     }
@@ -361,11 +361,11 @@ Matrix Aurora::interp1(const Matrix& aX, const Matrix& aV, const Matrix& aX1, In
     const int nx = aX.getDimSize(0);
     const int ny = 1;
     int nx1 = aX1.getDimSize(0);
-    std::shared_ptr<double> resultData = std::shared_ptr<double>(Aurora::malloc(nx1), Aurora::free);
+    std::shared_ptr<float> resultData = std::shared_ptr<float>(Aurora::malloc(nx1), Aurora::free);
     std::vector<int> resultInfo = {nx1};
     Aurora::Matrix result(resultData, resultInfo);
     DFTaskPtr task ;
-    int status = dfdNewTask1D(&task, nx, aX.getData(), DF_NO_HINT, ny, aV.getData(), DF_NO_HINT);
+    int status = dfsNewTask1D(&task, nx, aX.getData(), DF_NO_HINT, ny, aV.getData(), DF_NO_HINT);
     if (status != DF_STATUS_OK)
     {
         return Matrix();
@@ -382,20 +382,20 @@ Matrix Aurora::interp1(const Matrix& aX, const Matrix& aV, const Matrix& aX1, In
         sorder = DF_PP_LINEAR;
         stype = DF_PP_BESSEL;
     }
-    double* scoeffs = Aurora::malloc(ny * (nx-1) * sorder);
-     status = dfdEditPPSpline1D(task, sorder,DF_PP_NATURAL , DF_BC_NOT_A_KNOT, 0, DF_NO_IC, 0, scoeffs, DF_NO_HINT);
+    float* scoeffs = Aurora::malloc(ny * (nx-1) * sorder);
+     status = dfsEditPPSpline1D(task, sorder,DF_PP_NATURAL , DF_BC_NOT_A_KNOT, 0, DF_NO_IC, 0, scoeffs, DF_NO_HINT);
     if (status != DF_STATUS_OK)
     {
         return Matrix();
     }
-    status = dfdConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD );
+    status = dfsConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD );
     if (status != DF_STATUS_OK)
     {
         return Matrix();
     }
 
     int dorder = 1;
-    status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, nx1, aX1.getData(), DF_NO_HINT, 1, &dorder,
+    status = dfsInterpolate1D(task, DF_INTERP, DF_METHOD_PP, nx1, aX1.getData(), DF_NO_HINT, 1, &dorder,
                               DF_NO_APRIORI_INFO, resultData.get(), DF_MATRIX_STORAGE_ROWS, nullptr);
 
     status = dfDeleteTask(&task);
@@ -412,20 +412,20 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes)
     }
     int complexStep = aMatrix.getValueType();
     int originalDataSize = aMatrix.getDataSize() * complexStep;
-    double* resultData = Aurora::malloc(originalDataSize * aRowTimes * aColumnTimes);
+    float* resultData = Aurora::malloc(originalDataSize * aRowTimes * aColumnTimes);
     int row = aMatrix.getDimSize(0);
     int column = aMatrix.getDimSize(1);
-    double* originalData = aMatrix.getData();
-    double* resultDataTemp = resultData;
+    float* originalData = aMatrix.getData();
+    float* resultDataTemp = resultData;
     size_t step = row*complexStep;
     #pragma omp parallel for
     for(int i=0; i<column; ++i)
     {
-        double* origninalStart = originalData + i * step;
+        float* origninalStart = originalData + i * step;
         for(int j=0; j<aRowTimes; ++j)
         {
-            double* copyStart = resultData + step * (i*aRowTimes + j);
-            cblas_dcopy(step, origninalStart, 1, copyStart, 1);
+            float* copyStart = resultData + step * (i*aRowTimes + j);
+            cblas_scopy(step, origninalStart, 1, copyStart, 1);
             //resultDataTemp += row*complexStep;
         }
         //originalData += step;
@@ -436,7 +436,7 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes)
     #pragma omp parallel for
     for(int i=1; i<aColumnTimes; ++i)
     {
-        cblas_dcopy(step, resultData, 1, resultData + i*step, 1);
+        cblas_scopy(step, resultData, 1, resultData + i*step, 1);
     }
 
     std::vector<int> resultInfo;
@@ -447,7 +447,7 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes)
         resultInfo.push_back(column);
     }
 
-    return Matrix(std::shared_ptr<double>(resultData, Aurora::free),resultInfo, aMatrix.getValueType());
+    return Matrix(std::shared_ptr<float>(resultData, Aurora::free),resultInfo, aMatrix.getValueType());
 }
 
 Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int aSliceTimes)
@@ -460,11 +460,11 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int
     int complexStep = aMatrix.getValueType();
     Matrix resultTemp = Aurora::repmat(aMatrix, aRowTimes, aColumnTimes);
     int resultTempDataSize = resultTemp.getDataSize() * complexStep;
-    double* resultData = Aurora::malloc(resultTempDataSize * aSliceTimes);
+    float* resultData = Aurora::malloc(resultTempDataSize * aSliceTimes);
     std::copy(resultTemp.getData(), resultTemp.getData() + resultTempDataSize, resultData);
     for(int i=1; i<aSliceTimes; ++i)
     {
-        cblas_dcopy(resultTempDataSize, resultData, 1, resultData + i*resultTempDataSize, 1);
+        cblas_scopy(resultTempDataSize, resultData, 1, resultData + i*resultTempDataSize, 1);
     }
     std::vector<int> resultInfo;
     int row = resultTemp.getDimSize(0);
@@ -479,7 +479,7 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int
         resultInfo.push_back(aSliceTimes);
     }
 
-    return Matrix(std::shared_ptr<double>(resultData, Aurora::free), resultInfo, aMatrix.getValueType());
+    return Matrix(std::shared_ptr<float>(resultData, Aurora::free), resultInfo, aMatrix.getValueType());
 }
 
 Matrix Aurora::repmat3d(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int aSliceTimes)
@@ -488,27 +488,27 @@ Matrix Aurora::repmat3d(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, i
     {
         return Matrix();
     }
-    double* start = aMatrix.getData();
+    float* start = aMatrix.getData();
     int rows = aMatrix.getDimSize(0);
     int columns = aMatrix.getDimSize(1);
     int slices = aMatrix.getDimSize(2);
-    double* extended2DimsData = Aurora::malloc(rows * columns * aRowTimes * aColumnTimes * slices);
+    float* extended2DimsData = Aurora::malloc(rows * columns * aRowTimes * aColumnTimes * slices);
     Matrix extended2DimsMatrix = Matrix::New(extended2DimsData, aRowTimes*rows, aColumnTimes*columns, slices);
 
     for(int i=0; i<aMatrix.getDimSize(2); ++i)
     {
         Matrix dim2Matrix = Matrix::copyFromRawData(start, rows, columns);
          Matrix extendedTemp = repmat(dim2Matrix, aRowTimes, aColumnTimes);
-         cblas_dcopy(extendedTemp.getDataSize(), extendedTemp.getData(), 1, extended2DimsData, 1);
+         cblas_scopy(extendedTemp.getDataSize(), extendedTemp.getData(), 1, extended2DimsData, 1);
          extended2DimsData += extendedTemp.getDataSize();
         start += columns * rows;
     }
 
-    double* extended3DimsData = Aurora::malloc(rows * columns * aRowTimes * aColumnTimes * aSliceTimes * slices);
+    float* extended3DimsData = Aurora::malloc(rows * columns * aRowTimes * aColumnTimes * aSliceTimes * slices);
     Matrix result = Matrix::New(extended3DimsData, aRowTimes*rows, aColumnTimes*columns, slices * aSliceTimes);
     for(int i=0;i<aSliceTimes;++i)
     {
-        cblas_dcopy(extended2DimsMatrix.getDataSize(), extended2DimsMatrix.getData(), 1, extended3DimsData, 1);
+        cblas_scopy(extended2DimsMatrix.getDataSize(), extended2DimsMatrix.getData(), 1, extended3DimsData, 1);
         extended3DimsData+=extended2DimsMatrix.getDataSize();
     }
     return result;
@@ -517,16 +517,16 @@ Matrix Aurora::repmat3d(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, i
 
 Matrix Aurora::polyval(const Matrix &aP, const Matrix &aX) {
     auto result = malloc(aX.getDataSize());
-    auto powArg = new double[aP.getDataSize()];
+    auto powArg = new float[aP.getDataSize()];
     for (int j = aP.getDataSize(), i = 0; j > 0; --j, ++i) {
-        powArg[i] = (double) (j - 1);
+        powArg[i] = (float) (j - 1);
     }
-    auto temp = new double[aP.getDataSize()];
+    auto temp = new float[aP.getDataSize()];
     for (int i = 0; i < aX.getDataSize(); ++i) {
-        vdPowI(aP.getDataSize(), aX.getData() + i, 0, powArg, 1, temp, 1);
-        vdMul(aP.getDataSize(), aP.getData(), temp, temp);
-        Eigen::Map<Eigen::VectorXd> vd(temp,aP.getDataSize());
-        result[i] = vd.array().sum();
+        vsPowI(aP.getDataSize(), aX.getData() + i, 0, powArg, 1, temp, 1);
+        vsMul(aP.getDataSize(), aP.getData(), temp, temp);
+        Eigen::Map<Eigen::VectorXf> vs(temp,aP.getDataSize());
+        result[i] = vs.array().sum();
     }
     delete[] powArg;
     delete[] temp;
@@ -537,13 +537,13 @@ Matrix Aurora::polyval(const Matrix &aP, const Matrix &aX) {
 Matrix Aurora::log(const Matrix& aMatrix, int aBaseNum)
 {
     size_t size = aMatrix.getDataSize();
-    double* data = Aurora::malloc(size);
-    vdLn(size, aMatrix.getData(), data);
+    float* data = Aurora::malloc(size);
+    vsLn(size, aMatrix.getData(), data);
     if(aBaseNum != -1)
     {
-        double baseNum = aBaseNum;
-        double temp;
-        vdLn(1, &baseNum, &temp);
+        float baseNum = aBaseNum;
+        float temp;
+        vsLn(1, &baseNum, &temp);
         for (size_t i = 0; i < size; i++)
         {
             data[i] /= temp;
@@ -555,21 +555,21 @@ Matrix Aurora::log(const Matrix& aMatrix, int aBaseNum)
 Matrix Aurora::exp(const Matrix& aMatrix)
 {
     size_t size = aMatrix.getDataSize();
-    double* data;
+    float* data;
     if (aMatrix.isComplex())
     {
         data = Aurora::malloc(size, true);
-        vzExp(size, (MKL_Complex16*)aMatrix.getData(), (MKL_Complex16*)data);
+        vcExp(size, (MKL_Complex8*)aMatrix.getData(), (MKL_Complex8*)data);
     }
     else
     {
         data = Aurora::malloc(size);
-        vdExp(size, aMatrix.getData(), data);
+        vsExp(size, aMatrix.getData(), data);
     }
     return Matrix::New(data, aMatrix);
 }
 
-Matrix Aurora::mod(const Matrix& aMatrix, double aValue)
+Matrix Aurora::mod(const Matrix& aMatrix, float aValue)
 {
     if(aMatrix.isComplex() || aMatrix.isNull())
     {
@@ -577,8 +577,8 @@ Matrix Aurora::mod(const Matrix& aMatrix, double aValue)
     }
 
     size_t size = aMatrix.getDataSize();
-    double* matrixData = aMatrix.getData();
-    double* resultData = Aurora::malloc(size);
+    float* matrixData = aMatrix.getData();
+    float* resultData = Aurora::malloc(size);
     for(size_t i=0; i<size; ++i)
     {
         resultData[i] = fmod(matrixData[i], aValue);
@@ -595,9 +595,9 @@ Matrix Aurora::acos(const Matrix& aMatrix)
     }
 
     size_t size = aMatrix.getDataSize();
-    double* matrixData = aMatrix.getData();
-    double* resultData = Aurora::malloc(size);
-    vdAcos(size, matrixData, resultData);
+    float* matrixData = aMatrix.getData();
+    float* resultData = Aurora::malloc(size);
+    vsAcos(size, matrixData, resultData);
     return Matrix::New(resultData, aMatrix);
 }
 
@@ -609,9 +609,9 @@ Matrix Aurora::acosd(const Matrix& aMatrix)
     }
 
     size_t size = aMatrix.getDataSize();
-    double* matrixData = aMatrix.getData();
-    double* resultData = Aurora::malloc(size);
-    vdAcos(size, matrixData, resultData);
+    float* matrixData = aMatrix.getData();
+    float* resultData = Aurora::malloc(size);
+    vsAcos(size, matrixData, resultData);
     for(size_t i=0; i<size; ++i)
     {
         resultData[i] = resultData[i] * 180 / PI;
@@ -626,12 +626,12 @@ Matrix Aurora::conj(const Matrix& aMatrix)
         return Matrix::copyFromRawData(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1),aMatrix.getDimSize(2));
     }
     size_t size = aMatrix.getDataSize();
-    double* data = malloc(size,true);
-    vzConj(size,(MKL_Complex16*)aMatrix.getData(), (MKL_Complex16*)data);
+    float* data = malloc(size,true);
+    vcConj(size,(MKL_Complex8*)aMatrix.getData(), (MKL_Complex8*)data);
     return Matrix::New(data, aMatrix);
 }
 
-double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
+float Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
 {
     if(aMatrix.isNull())
     {
@@ -647,10 +647,10 @@ double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
     int row = aMatrix.getDimSize(0);
     if (aNormMethod == NormMethod::Norm1)
     {
-        double value = 0;
+        float value = 0;
         for(int i=0; i<column; ++i)
         {
-            double temp = Aurora::sum(abs(aMatrix($,i,$).toMatrix())).getData()[0];
+            float temp = Aurora::sum(abs(aMatrix($,i,$).toMatrix())).getData()[0];
             if(temp > value)
             {
                 value = temp;
@@ -660,7 +660,7 @@ double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
     }
     else if(aNormMethod == NormMethod::NormF)
     {
-        return cblas_dnrm2(size, aMatrix.getData(), 1);
+        return cblas_snrm2(size, aMatrix.getData(), 1);
     }
     else if(aNormMethod == NormMethod::Norm2)
     {
@@ -669,20 +669,20 @@ double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
         {
             if(aMatrix.isComplex())
             {
-                Eigen::Map<Eigen::MatrixXcd> eMatrix((MKL_Complex16*)aMatrix.getData(), row, column);
-                Eigen::JacobiSVD<Eigen::MatrixXcd> svd(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
-                return svd.singularValues()(0);
+                Eigen::Map<Eigen::MatrixXcf> eMatrix((MKL_Complex8*)aMatrix.getData(), row, column);
+                Eigen::JacobiSVD<Eigen::MatrixXcf> svs(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
+                return svs.singularValues()(0);
             }
             else
             {
-                Eigen::Map<Eigen::MatrixXd> eMatrix(aMatrix.getData(), row, column);
-                Eigen::JacobiSVD<Eigen::MatrixXd> svd(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
-                return svd.singularValues()(0);
+                Eigen::Map<Eigen::MatrixXf> eMatrix(aMatrix.getData(), row, column);
+                Eigen::JacobiSVD<Eigen::MatrixXf> svs(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
+                return svs.singularValues()(0);
             }
         }
         else
         {
-            return cblas_dnrm2(size, aMatrix.getData(), 1);
+            return cblas_snrm2(size, aMatrix.getData(), 1);
         }
     }
 
@@ -699,17 +699,17 @@ Matrix Aurora::transpose(const Matrix& aMatrix)
     size_t size = aMatrix.getDataSize();
     int row = aMatrix.getDimSize(0);
     int col = aMatrix.getDimSize(1);
-    double* resultData;
-    double* data = aMatrix.getData();
+    float* resultData;
+    float* data = aMatrix.getData();
     if(aMatrix.isComplex())
     {
         resultData = Aurora::malloc(size, true);
-        mkl_zomatcopy('C', 'T',row ,col, 1.0, (MKL_Complex16*)data, row, (MKL_Complex16*)resultData, col);
+        mkl_comatcopy('C', 'T',row ,col, 1.0, (MKL_Complex8*)data, row, (MKL_Complex8*)resultData, col);
     }
     else
     {
         resultData = Aurora::malloc(size);
-        mkl_domatcopy('C', 'T',row ,col, 1.0, data, row, resultData, col);
+        mkl_somatcopy('C', 'T',row ,col, 1.0, data, row, resultData, col);
     }
 
     return Matrix::New(resultData,col,row,1,aMatrix.getValueType());
@@ -728,12 +728,12 @@ Matrix Aurora::horzcat(const Matrix& aMatrix1, const Matrix& aMatrix2)
     int row = aMatrix1.getDimSize(0);
     size_t size1= row*column1;
     size_t size2= row*column2;
-    double* resultData = Aurora::malloc(aMatrix1.getDataSize() + aMatrix2.getDataSize(),aMatrix1.getValueType());
+    float* resultData = Aurora::malloc(aMatrix1.getDataSize() + aMatrix2.getDataSize(),aMatrix1.getValueType());
     size_t sliceStride = row*(column1+column2);
     for (size_t i = 0; i < slice; i++)
     {
-        cblas_dcopy(size1, aMatrix1.getData()+i*size1 , 1, resultData + i*sliceStride, 1);
-        cblas_dcopy(size2, aMatrix2.getData()+i*size2, 1, resultData + i*sliceStride + size1, 1);
+        cblas_scopy(size1, aMatrix1.getData()+i*size1 , 1, resultData + i*sliceStride, 1);
+        cblas_scopy(size2, aMatrix2.getData()+i*size2, 1, resultData + i*sliceStride + size1, 1);
     }
     return Matrix::New(resultData, row, column1+column2, slice, aMatrix1.getValueType());
 }
@@ -751,9 +751,9 @@ Matrix Aurora::vertcat(const Matrix& aMatrix1, const Matrix& aMatrix2){
     int column = aMatrix1.getDimSize(1);
     size_t size1= aMatrix1.getDataSize();
     size_t size2= aMatrix2.getDataSize();
-    double* resultData = Aurora::malloc(size1 + size2,aMatrix1.getValueType());
-    cblas_dcopy_batch_strided(row1, aMatrix1.getData(), 1,row1, resultData, 1, row1+row2, column*slice);
-    cblas_dcopy_batch_strided(row2, aMatrix2.getData(), 1,row2, resultData + row1, 1, row1+row2, column*slice);
+    float* resultData = Aurora::malloc(size1 + size2,aMatrix1.getValueType());
+    cblas_scopy_batch_strided(row1, aMatrix1.getData(), 1,row1, resultData, 1, row1+row2, column*slice);
+    cblas_scopy_batch_strided(row2, aMatrix2.getData(), 1,row2, resultData + row1, 1, row1+row2, column*slice);
 
     return Matrix::New(resultData, row1+row2, column, slice, aMatrix1.getValueType());
 }
@@ -766,7 +766,7 @@ Matrix Aurora::vecnorm(const Matrix& aMatrix, NormMethod aNormMethod, int aDim)
         return Matrix();
     }
     int column = aMatrix.getDimSize(1);
-    double* resultData = Aurora::malloc(column);
+    float* resultData = Aurora::malloc(column);
     for(int i=0; i<column; ++i)
     {
         resultData[i] = norm(aMatrix($,i,$).toMatrix(), aNormMethod);
@@ -775,10 +775,10 @@ Matrix Aurora::vecnorm(const Matrix& aMatrix, NormMethod aNormMethod, int aDim)
     return Matrix::New(resultData,column);
 }
 
-Matrix Aurora::linspace(double aStart, double aEnd, int aNum)
+Matrix Aurora::linspace(float aStart, float aEnd, int aNum)
 {
-    double step = (aEnd - aStart) / (aNum - 1);
-    double* resultData = Aurora::malloc(aNum);
+    float step = (aEnd - aStart) / (aNum - 1);
+    float* resultData = Aurora::malloc(aNum);
     for (int i = 0; i < aNum; i++)
     {
         resultData[i] = aStart + step * i;
@@ -796,10 +796,10 @@ Matrix Aurora::auroraUnion(const Matrix& aMatrix1, const Matrix& aMatrix2)
 
     size_t size1= aMatrix1.getDataSize();
     size_t size2= aMatrix2.getDataSize();
-    double* resultData = Aurora::malloc(size1 + size2);
-    cblas_dcopy(size1, aMatrix1.getData(), 1, resultData, 1);
-    cblas_dcopy(size2, aMatrix2.getData(), 1, resultData + size1, 1);
-    std::vector<double> vector(resultData,  resultData + size1 + size2);
+    float* resultData = Aurora::malloc(size1 + size2);
+    cblas_scopy(size1, aMatrix1.getData(), 1, resultData, 1);
+    cblas_scopy(size2, aMatrix2.getData(), 1, resultData + size1, 1);
+    std::vector<float> vector(resultData,  resultData + size1 + size2);
     Aurora::free(resultData);
     std::sort(vector.begin(), vector.end());
     auto last = std::unique(vector.begin(), vector.end());
@@ -817,12 +817,12 @@ Matrix Aurora::intersect(const Matrix& aMatrix1, const Matrix& aMatrix2)
 
     size_t size1= aMatrix1.getDataSize();
     size_t size2= aMatrix2.getDataSize();
-    std::vector<double> vector1(aMatrix1.getData(),  aMatrix1.getData() + size1);
-    std::vector<double> vector2(aMatrix2.getData(),  aMatrix2.getData() + size2);
+    std::vector<float> vector1(aMatrix1.getData(),  aMatrix1.getData() + size1);
+    std::vector<float> vector2(aMatrix2.getData(),  aMatrix2.getData() + size2);
     std::sort(vector1.begin(), vector1.end());
     std::sort(vector2.begin(), vector2.end());
 
-    std::vector<double> intersection;
+    std::vector<float> intersection;
     std::set_intersection(vector1.begin(), vector1.end(),
                           vector2.begin(), vector2.end(),
                           std::back_inserter(intersection));
@@ -840,7 +840,7 @@ Matrix Aurora::intersect(const Matrix& aMatrix1, const Matrix& aMatrix2, Matrix&
     Matrix result = intersect(aMatrix1,aMatrix2);
 
     size_t size = result.getDataSize();
-    double* iaResult = Aurora::malloc(size);
+    float* iaResult = Aurora::malloc(size);
     for(size_t i=0; i<size; ++i)
     {
         for(size_t j=0; j<aMatrix1.getDataSize(); ++j)
@@ -866,10 +866,10 @@ Matrix Aurora::xcorr(const Matrix& aMatrix1, const Matrix& aMatrix2)
     }
     size_t matrixSize = aMatrix1.getDataSize();
     size_t resultSize = 2 * matrixSize -1;
-    double* resultData = Aurora::malloc(resultSize);
+    float* resultData = Aurora::malloc(resultSize);
     for(int i=0;i<matrixSize;++i)
     {
-        double data = 0;
+        float data = 0;
         for(int j=0;j<i+1;++j)
         {
             data+= aMatrix1[j] * aMatrix2[matrixSize-i-1+j];
@@ -879,7 +879,7 @@ Matrix Aurora::xcorr(const Matrix& aMatrix1, const Matrix& aMatrix2)
 
     for(int i=0;i<matrixSize-1;++i)
     {
-        double result = 0;
+        float result = 0;
         for(int j=0;j<i+1;++j)
         {
             result+= aMatrix1[matrixSize-i-1+j]*aMatrix2[j];
@@ -900,29 +900,29 @@ Matrix Aurora::deleteColumn(const Matrix& aMatrix, int aColumnIndex)
         return aMatrix;
     }
 
-    double* resultData = Aurora::malloc(rows* (columns-1));
+    float* resultData = Aurora::malloc(rows* (columns-1));
     if(aColumnIndex == 0)
     {
-        cblas_dcopy(rows* (columns-1), aMatrix.getData() + rows, 1, resultData, 1);
+        cblas_scopy(rows* (columns-1), aMatrix.getData() + rows, 1, resultData, 1);
     }
 
     else if(aColumnIndex == (columns - 1))
     {
-        cblas_dcopy(rows* (columns-1), aMatrix.getData(), 1, resultData, 1);
+        cblas_scopy(rows* (columns-1), aMatrix.getData(), 1, resultData, 1);
     }
     else
     {
-        cblas_dcopy(rows * aColumnIndex, aMatrix.getData(), 1, resultData, 1);
-        cblas_dcopy(rows * (columns - aColumnIndex - 1), aMatrix.getData() + rows * (aColumnIndex + 1), 1, resultData + rows * aColumnIndex, 1);
+        cblas_scopy(rows * aColumnIndex, aMatrix.getData(), 1, resultData, 1);
+        cblas_scopy(rows * (columns - aColumnIndex - 1), aMatrix.getData() + rows * (aColumnIndex + 1), 1, resultData + rows * aColumnIndex, 1);
     }
 
     return Matrix::New(resultData, rows, columns-1);
 }
 
-Matrix Aurora::createVectorMatrix(double aStartValue, double aStepValue, double aEndValue)
+Matrix Aurora::createVectorMatrix(float aStartValue, float aStepValue, float aEndValue)
 {
-    std::vector<double> matrixData;
-    double tempValue = aStartValue;
+    std::vector<float> matrixData;
+    float tempValue = aStartValue;
     matrixData.push_back(tempValue);
     long long compare1 = std::round(aEndValue * 10e13);
     long long compare2 = std::round(tempValue * 10e13);
@@ -945,28 +945,28 @@ Matrix Aurora::reshape(const Matrix& aMatrix, int aRows, int aColumns, int aSlic
     }
     return Matrix::copyFromRawData(aMatrix.getData(),aRows,aColumns,aSlices);
 }
-void Aurora::nantoval(Matrix& aMatrix, double val2) {
-    Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
+void Aurora::nantoval(Matrix& aMatrix, float val2) {
+    Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
     srcV = srcV.array().isNaN().select(val2,srcV);
 }
 
 Matrix Aurora::isnan(const Matrix& aMatrix){
-    Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
     auto result = zeros(aMatrix.getDimSize(0),aMatrix.getDimSize(1),aMatrix.getDimSize(2));
-    Eigen::Map<Eigen::VectorXd> resultV(result.getData(),result.getDataSize());
+    Eigen::Map<Eigen::VectorXf> resultV(result.getData(),result.getDataSize());
     resultV = srcV.array().isNaN().select(1.0,resultV);
     return result;
 }
 
 Matrix Aurora::isfinite(const Matrix& aMatrix){
-    Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
     auto result = zeros(aMatrix.getDimSize(0),aMatrix.getDimSize(1),aMatrix.getDimSize(2));
-    Eigen::Map<Eigen::VectorXd> resultV(result.getData(),result.getDataSize());
+    Eigen::Map<Eigen::VectorXf> resultV(result.getData(),result.getDataSize());
     resultV = srcV.array().isFinite().select(1.0,resultV);
     return result;
 }
 
-void Aurora::padding(Matrix &aMatrix, int aIndex, double aValue)
+void Aurora::padding(Matrix &aMatrix, int aIndex, float aValue)
 {
     if(aMatrix.isNull() || !aMatrix.isVector())
     {
@@ -978,16 +978,16 @@ void Aurora::padding(Matrix &aMatrix, int aIndex, double aValue)
         return;
     }
     int size = (aIndex+1);
-    double* newData = malloc(size,aMatrix.isComplex());
-    cblas_dcopy(aMatrix.getDataSize()*aMatrix.getValueType(),
+    float* newData = malloc(size,aMatrix.isComplex());
+    cblas_scopy(aMatrix.getDataSize()*aMatrix.getValueType(),
         aMatrix.getData(),1,newData,1);
-    cblas_dcopy((size-aMatrix.getDataSize())*aMatrix.getValueType(),
+    cblas_scopy((size-aMatrix.getDataSize())*aMatrix.getValueType(),
         &aValue,0,newData+aMatrix.getDataSize()*aMatrix.getValueType(),1);
     aMatrix = Matrix::New(newData,size,1,1,aMatrix.getValueType());
 }
 
-void Aurora::compareSet(Matrix& aMatrix,double compareValue, double newValue,CompareOp op){
-    Eigen::Map<Eigen::VectorXd> v(aMatrix.getData(),aMatrix.getDataSize());
+void Aurora::compareSet(Matrix& aMatrix,float compareValue, float newValue,CompareOp op){
+    Eigen::Map<Eigen::VectorXf> v(aMatrix.getData(),aMatrix.getDataSize());
     switch (op) {
         case EQ:
             v = (v.array() == compareValue).select(newValue, v);
@@ -1010,9 +1010,9 @@ void Aurora::compareSet(Matrix& aMatrix,double compareValue, double newValue,Com
     }
 }
 
-void Aurora::compareSet(Matrix& aMatrix,Matrix& aCompareMatrix,double compareValue, double newValue,CompareOp op){
-    Eigen::Map<Eigen::VectorXd> v(aMatrix.getData(),aMatrix.getDataSize());
-    Eigen::Map<Eigen::VectorXd> c(aCompareMatrix.getData(),aCompareMatrix.getDataSize());
+void Aurora::compareSet(Matrix& aMatrix,Matrix& aCompareMatrix,float compareValue, float newValue,CompareOp op){
+    Eigen::Map<Eigen::VectorXf> v(aMatrix.getData(),aMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> c(aCompareMatrix.getData(),aCompareMatrix.getDataSize());
     switch (op) {
         case EQ:
             v = (c.array() == compareValue).select(newValue, v);
@@ -1035,9 +1035,9 @@ void Aurora::compareSet(Matrix& aMatrix,Matrix& aCompareMatrix,double compareVal
     }
 }
 
-void Aurora::compareSet(Matrix& aValueAndCompareMatrix,Matrix& aOtherCompareMatrix, double newValue,CompareOp op){
-    Eigen::Map<Eigen::VectorXd> v(aValueAndCompareMatrix.getData(),aValueAndCompareMatrix.getDataSize());
-    Eigen::Map<Eigen::VectorXd> c(aOtherCompareMatrix.getData(),aOtherCompareMatrix.getDataSize());
+void Aurora::compareSet(Matrix& aValueAndCompareMatrix,Matrix& aOtherCompareMatrix, float newValue,CompareOp op){
+    Eigen::Map<Eigen::VectorXf> v(aValueAndCompareMatrix.getData(),aValueAndCompareMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> c(aOtherCompareMatrix.getData(),aOtherCompareMatrix.getDataSize());
     switch (op) {
         case EQ:
             v = (v.array() == c.array()).select(newValue, v);
@@ -1059,9 +1059,9 @@ void Aurora::compareSet(Matrix& aValueAndCompareMatrix,Matrix& aOtherCompareMatr
             break;
     }
 }
-void Aurora::compareSet(Matrix& aCompareMatrix,double compareValue, Matrix& aNewValueMatrix,CompareOp op){
-    Eigen::Map<Eigen::VectorXd> v(aCompareMatrix.getData(),aCompareMatrix.getDataSize());
-    Eigen::Map<Eigen::VectorXd> nv(aNewValueMatrix.getData(),aNewValueMatrix.getDataSize());
+void Aurora::compareSet(Matrix& aCompareMatrix,float compareValue, Matrix& aNewValueMatrix,CompareOp op){
+    Eigen::Map<Eigen::VectorXf> v(aCompareMatrix.getData(),aCompareMatrix.getDataSize());
+    Eigen::Map<Eigen::VectorXf> nv(aNewValueMatrix.getData(),aNewValueMatrix.getDataSize());
     switch (op) {
         case EQ:
             v = (v.array() == compareValue).select(nv, v);
@@ -1115,13 +1115,13 @@ Matrix Aurora::uniqueByRows(const Matrix& aMatrix, Matrix& aIndexResult)
         return Matrix();
     }
     Matrix transposeMatrix = transpose(aMatrix);
-    double* transposeMatrixData = transposeMatrix.getData();
+    float* transposeMatrixData = transposeMatrix.getData();
     int rows = transposeMatrix.getDimSize(0);
     int columns = transposeMatrix.getDimSize(1);
     std::vector<Matrix> uniqueResult;
     for(int i=1; i<=columns; ++i)
     {
-        Matrix rowData = Matrix::fromRawData(new double[rows], rows);
+        Matrix rowData = Matrix::fromRawData(new float[rows], rows);
         std::copy(transposeMatrixData, transposeMatrixData + rows, rowData.getData());
         uniqueResult.push_back(rowData);
         transposeMatrixData += rows;
@@ -1132,7 +1132,7 @@ Matrix Aurora::uniqueByRows(const Matrix& aMatrix, Matrix& aIndexResult)
     auto uniqueIndex = std::unique(uniqueResult.begin(), uniqueResult.end(), isEqual);
     uniqueResult.erase(uniqueIndex, uniqueResult.end());
     
-    std::vector<double> indexResultData;
+    std::vector<float> indexResultData;
     for(int i=0; i<matrixsCopy.size();++i)
     {
         auto index = lower_bound(uniqueResult.begin(), uniqueResult.end(), matrixsCopy[i], compare);
@@ -1142,7 +1142,7 @@ Matrix Aurora::uniqueByRows(const Matrix& aMatrix, Matrix& aIndexResult)
 
     aIndexResult = Matrix::copyFromRawData(indexResultData.data(), indexResultData.size());
     size_t resultSize = rows * uniqueResult.size();
-    double* resultData = new double[resultSize];
+    float* resultData = new float[resultSize];
     Matrix result = Matrix::fromRawData(resultData, rows, uniqueResult.size()) ;
     for(size_t i=0; i<uniqueResult.size(); ++i)
     {
diff --git a/src/Function1D.h b/src/Function1D.h
index 959c005..ccb10f6 100644
--- a/src/Function1D.h
+++ b/src/Function1D.h
@@ -63,7 +63,7 @@ namespace Aurora {
 
     Matrix exp(const Matrix& aMatrix);
 
-    Matrix mod(const Matrix& aMatrix, double aValue);
+    Matrix mod(const Matrix& aMatrix, float aValue);
 
     Matrix acos(const Matrix& aMatrix);
 
@@ -71,7 +71,7 @@ namespace Aurora {
 
     Matrix conj(const Matrix& aMatrix);
 
-    double norm(const Matrix& aMatrix, NormMethod aNormMethod);
+    float norm(const Matrix& aMatrix, NormMethod aNormMethod);
 
     Matrix transpose(const Matrix& aMatrix);
 
@@ -81,7 +81,7 @@ namespace Aurora {
 
     Matrix vecnorm(const Matrix& aMatrix, NormMethod aNormMethod, int aDim);
 
-    Matrix linspace(double aStart, double aEnd, int aNum);
+    Matrix linspace(float aStart, float aEnd, int aNum);
 
     Matrix auroraUnion(const Matrix& aMatrix1, const Matrix& aMatrix2);
 
@@ -93,7 +93,7 @@ namespace Aurora {
 
     Matrix deleteColumn(const Matrix& aMatrix, int aColumnIndex);
 
-    Matrix createVectorMatrix(double aStartValue, double aStepValue, double aEndValue);
+    Matrix createVectorMatrix(float aStartValue, float aStepValue, float aEndValue);
 
     Matrix uniqueByRows(const Matrix& aMatrix, Matrix& aIndexResult);
     /**
@@ -119,7 +119,7 @@ namespace Aurora {
      * @param aMatrix 向量
      * @param val 指定值
      */
-    void nantoval(Matrix& aMatrix,double val);
+    void nantoval(Matrix& aMatrix,float val);
 
 
     Matrix isnan(const Matrix& aMatrix);
@@ -133,17 +133,17 @@ namespace Aurora {
      * @param aIndex 长度索引
      * @param aValue 指定值
      */
-    void padding(Matrix& aMatrix, int aIndex, double aValue);
+    void padding(Matrix& aMatrix, int aIndex, float aValue);
 
     Matrix auroraNot(const Matrix& aMatrix);
     Matrix auroraNot(Matrix&& aMatrix);
     enum CompareOp{
         EQ,GT,LT,NG,NL,NE
     };
-    void compareSet(Matrix& aValueMatrix,double compareValue, double newValue,CompareOp op);
-    void compareSet(Matrix& aValueMatrix,Matrix& aCompareMatrix,double compareValue, double newValue,CompareOp op);
-    void compareSet(Matrix& aDesAndCompareMatrix,Matrix& aOtherCompareMatrix, double newValue,CompareOp op);
-    void compareSet(Matrix& aCompareMatrix,double compareValue, Matrix& aNewValueMatrix,CompareOp op);
+    void compareSet(Matrix& aValueMatrix,float compareValue, float newValue,CompareOp op);
+    void compareSet(Matrix& aValueMatrix,Matrix& aCompareMatrix,float compareValue, float newValue,CompareOp op);
+    void compareSet(Matrix& aDesAndCompareMatrix,Matrix& aOtherCompareMatrix, float newValue,CompareOp op);
+    void compareSet(Matrix& aCompareMatrix,float compareValue, Matrix& aNewValueMatrix,CompareOp op);
     
     Matrix convertfp16tofloat(short* aData, int aRows, int aColumns);
 };
diff --git a/src/Function2D.cpp b/src/Function2D.cpp
index 7ec63ee..b4cc1d2 100644
--- a/src/Function2D.cpp
+++ b/src/Function2D.cpp
@@ -18,7 +18,7 @@
 
 using namespace Aurora;
 
-double Aurora::immse(const Aurora::Matrix &aImageA, const Aurora::Matrix &aImageB) {
+float Aurora::immse(const Aurora::Matrix &aImageA, const Aurora::Matrix &aImageB) {
     if (aImageA.getDims()!=2|| aImageB.getDims()!=2){
         std::cerr<<"Fail!immse args must all 2d matrix!";
         return 0.0;
@@ -33,9 +33,9 @@ double Aurora::immse(const Aurora::Matrix &aImageA, const Aurora::Matrix &aImage
     }
     int size = aImageA.getDataSize();
     auto temp = malloc(size);
-    vdSub(size, aImageA.getData(), aImageB.getData(), temp);
-    vdSqr(size, temp, temp);
-    double result = cblas_dasum(size, temp, 1) / (double) size;
+    vsSub(size, aImageA.getData(), aImageB.getData(), temp);
+    vsSqr(size, temp, temp);
+    float result = cblas_sasum(size, temp, 1) / (float) size;
     free(temp);
     return result;
 }
@@ -56,9 +56,9 @@ Aurora::Matrix Aurora::inv(const Aurora::Matrix &aMatrix) {
     int size = aMatrix.getDataSize();
     int *ipiv = new int[aMatrix.getDimSize(0)];
     auto result = malloc(size);
-    cblas_dcopy(size,aMatrix.getData(), 1,result, 1);
-    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
-    LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
+    cblas_scopy(size,aMatrix.getData(), 1,result, 1);
+    LAPACKE_sgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
+    LAPACKE_sgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
     delete[] ipiv;
     return Matrix::New(result,aMatrix);
 }
@@ -77,8 +77,8 @@ Aurora::Matrix Aurora::inv(Aurora::Matrix&& aMatrix) {
         return aMatrix;
     }
     int *ipiv = new int[aMatrix.getDimSize(0)];
-    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
-    LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
+    LAPACKE_sgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
+    LAPACKE_sgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
     delete[] ipiv;
     return aMatrix;
 }
@@ -104,10 +104,10 @@ Matrix Aurora::interp2(const Matrix& aX, const Matrix& aY, const Matrix& aV, con
     }
 
     int nx1 = aX1.getDimSize(0);
-    std::shared_ptr<double> resultData = std::shared_ptr<double>(Aurora::malloc(nx1), Aurora::free);
+    std::shared_ptr<float> resultData = std::shared_ptr<float>(Aurora::malloc(nx1), Aurora::free);
     for (int i = 0; i < nx1; ++i)
     {
-        std::shared_ptr<double> xResultData = std::shared_ptr<double>(Aurora::malloc(columnNum), Aurora::free);
+        std::shared_ptr<float> xResultData = std::shared_ptr<float>(Aurora::malloc(columnNum), Aurora::free);
         for(int j =0; j < columnNum; ++j)
         {
             xResultData.get()[j] = interp1(aY,aV($,j).toMatrix(),aY1(i).toMatrix(),aMethod).getData()[0];
@@ -137,14 +137,14 @@ Matrix Aurora::std(const Matrix &aMatrix) {
     auto std = Aurora::malloc(col);
     auto meanM = Aurora::mean(aMatrix);
     for (int i = 0; i < col; ++i) {
-        double *p = src.getData() + i * calc_size;
+        float *p = src.getData() + i * calc_size;
         
-        double mean = meanM[i];
-        vdSubI(calc_size, p, 1, &mean, 0, p, 1);
-        vdSqr(calc_size, p, p);
-        std[i] = cblas_dasum(calc_size, p, 1) / (calc_size - 1);
+        float mean = meanM[i];
+        vsSubI(calc_size, p, 1, &mean, 0, p, 1);
+        vsSqr(calc_size, p, p);
+        std[i] = cblas_sasum(calc_size, p, 1) / (calc_size - 1);
     }
-    vdSqrt(col, std, std);
+    vsSqrt(col, std, std);
     return Matrix::New(std, 1, col, aMatrix.getDimSize(2));
 }
 
@@ -162,21 +162,21 @@ Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& row
     switch (direction)
     {
         case All: {
-            Eigen::Map<Eigen::MatrixXd> retV(aMatrix.getData(), aMatrix.getDimSize(0), aMatrix.getDimSize(1));
-            double *ret = malloc(1);
+            Eigen::Map<Eigen::MatrixXf> retV(aMatrix.getData(), aMatrix.getDimSize(0), aMatrix.getDimSize(1));
+            float *ret = malloc(1);
             ret[0] = retV.array().minCoeff(&rowIdx, &colIdx);
             return Matrix::New(ret, 1);
         }
         case Row:
         {
-            Eigen::Map<Eigen::MatrixXd> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-            double * ret = malloc(aMatrix.getDimSize(0));
+            Eigen::Map<Eigen::MatrixXf> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+            float * ret = malloc(aMatrix.getDimSize(0));
 
             if (aMatrix.getDimSize(0) == 1){
                 ret[0] = srcMatrix.topRows(0).minCoeff(&rowIdx, &colIdx);
             }
             else{
-                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,aMatrix.getDimSize(0),1);
+                Eigen::Map<Eigen::MatrixXf> retMatrix(ret,aMatrix.getDimSize(0),1);
                 retMatrix = srcMatrix.rowwise().minCoeff();
             }
             return Matrix::New(ret,aMatrix.getDimSize(0),1);
@@ -184,14 +184,14 @@ Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& row
         case Column:
         default:
         {
-            Eigen::Map<Eigen::MatrixXd> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-            double * ret = malloc(aMatrix.getDimSize(1));
+            Eigen::Map<Eigen::MatrixXf> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+            float * ret = malloc(aMatrix.getDimSize(1));
 
             if (aMatrix.getDimSize(1) == 1){
                 ret[0] = srcMatrix.col(0).minCoeff(&rowIdx, &colIdx);
             }
             else {
-                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,1,aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> retMatrix(ret,1,aMatrix.getDimSize(1));
                 retMatrix = srcMatrix.colwise().minCoeff();
             }
             return Matrix::New(ret,1,aMatrix.getDimSize(1));
@@ -219,31 +219,31 @@ Matrix Aurora::min(const Matrix &aMatrix, const Matrix &aOther) {
     }
     //same shape
     if (aMatrix.compareShape(aOther)){
-        double* output = malloc(aMatrix.getDataSize());
-        vdFminI(aMatrix.getDataSize(),aMatrix.getData(),1,aOther.getData(),1,output,1);
+        float* output = malloc(aMatrix.getDataSize());
+        vsFminI(aMatrix.getDataSize(),aMatrix.getData(),1,aOther.getData(),1,output,1);
         return Matrix::New(output,aMatrix);
     }
     // one is scalar
     else if (aMatrix.getDataSize() == 1 || aOther.getDataSize() == 1){
-            double scalar = (aMatrix.getDataSize() == 1)?aMatrix.getData()[0]:aOther.getData()[0];
+            float scalar = (aMatrix.getDataSize() == 1)?aMatrix.getData()[0]:aOther.getData()[0];
             auto matrix = (aMatrix.getDataSize() == 1)?aOther:aMatrix;
-            double* output = malloc(matrix.getDataSize());
-            vdFminI(matrix.getDataSize(),matrix.getData(),1,&scalar,0,output,1);
+            float* output = malloc(matrix.getDataSize());
+            vsFminI(matrix.getDataSize(),matrix.getData(),1,&scalar,0,output,1);
             return Matrix::New(output,matrix);
     }
     else if (aMatrix.getDimSize(1) == 1 || aOther.getDimSize(0) == 1) {
         if (aMatrix.getDimSize(1) == 1){
-            double* output = malloc(aOther.getDataSize());
+            float* output = malloc(aOther.getDataSize());
             for (int i = 0; i < aOther.getDimSize(1); ++i) {
-                vdFminI(aMatrix.getDataSize(), aMatrix.getData(), 1, aOther.getData() + aOther.getDimSize(0) * i, 1,
+                vsFminI(aMatrix.getDataSize(), aMatrix.getData(), 1, aOther.getData() + aOther.getDimSize(0) * i, 1,
                         output + aOther.getDimSize(0) * i, 1);
             }
             return Matrix::New(output,aOther);
         }
         else{
-            double* output = malloc(aMatrix.getDataSize());
+            float* output = malloc(aMatrix.getDataSize());
             for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
-                vdFminI(aOther.getDataSize(), aOther.getData(), 1, aMatrix.getData() +  i, aMatrix.getDimSize(0),
+                vsFminI(aOther.getDataSize(), aOther.getData(), 1, aMatrix.getData() +  i, aMatrix.getDimSize(0),
                         output +  i, aOther.getDimSize(0));
             }
             return Matrix::New(output,aMatrix);
@@ -279,20 +279,20 @@ Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& row
     {
         case All:
         {
-            Eigen::Map<Eigen::MatrixXd> retV(calcMatrix.getData(), calcMatrix.getDimSize(0), calcMatrix.getDimSize(1));
-            double *ret = malloc(1);
+            Eigen::Map<Eigen::MatrixXf> retV(calcMatrix.getData(), calcMatrix.getDimSize(0), calcMatrix.getDimSize(1));
+            float *ret = malloc(1);
             ret[0] = retV.array().maxCoeff(&rowIdx, &colIdx);
             return Matrix::New(ret,1);
         }
         case Row:
         {
-            Eigen::Map<Eigen::MatrixXd> srcMatrix(calcMatrix.getData(),calcMatrix.getDimSize(0),calcMatrix.getDimSize(1));
-            double * ret = malloc(calcMatrix.getDimSize(0));
+            Eigen::Map<Eigen::MatrixXf> srcMatrix(calcMatrix.getData(),calcMatrix.getDimSize(0),calcMatrix.getDimSize(1));
+            float * ret = malloc(calcMatrix.getDimSize(0));
             if (calcMatrix.getDimSize(0) == 1){
                 ret[0] = srcMatrix.topRows(0).maxCoeff(&rowIdx, &colIdx);
             }
             else{
-                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,calcMatrix.getDimSize(0),1);
+                Eigen::Map<Eigen::MatrixXf> retMatrix(ret,calcMatrix.getDimSize(0),1);
                 retMatrix = srcMatrix.rowwise().maxCoeff();
             }
             return Matrix::New(ret,calcMatrix.getDimSize(0),1);
@@ -300,13 +300,13 @@ Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& row
         case Column:
         default:
         {
-            Eigen::Map<Eigen::MatrixXd> srcMatrix(calcMatrix.getData(),calcMatrix.getDimSize(0),calcMatrix.getDimSize(1));
-            double * ret = malloc(calcMatrix.getDimSize(1));
+            Eigen::Map<Eigen::MatrixXf> srcMatrix(calcMatrix.getData(),calcMatrix.getDimSize(0),calcMatrix.getDimSize(1));
+            float * ret = malloc(calcMatrix.getDimSize(1));
             if (calcMatrix.getDimSize(1) == 1){
                 ret[0] = srcMatrix.col(0).maxCoeff(&rowIdx, &colIdx);
             }
             else {
-                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,1,calcMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> retMatrix(ret,1,calcMatrix.getDimSize(1));
                 retMatrix = srcMatrix.colwise().maxCoeff();
             }
             return Matrix::New(ret,1,calcMatrix.getDimSize(1));
@@ -329,31 +329,31 @@ Matrix Aurora::max(const Matrix &aMatrix, const Matrix &aOther) {
     }
     //same shape
     if (aMatrix.compareShape(aOther)){
-        double* output = malloc(aMatrix.getDataSize());
-        vdFmaxI(aMatrix.getDataSize(),aMatrix.getData(),1,aOther.getData(),1,output,1);
+        float* output = malloc(aMatrix.getDataSize());
+        vsFmaxI(aMatrix.getDataSize(),aMatrix.getData(),1,aOther.getData(),1,output,1);
         return Matrix::New(output,aMatrix);
     }
     // one is scalar
     else if (aMatrix.getDataSize() == 1 || aOther.getDataSize() == 1){
-            double scalar = (aMatrix.getDataSize() == 1)?aMatrix.getData()[0]:aOther.getData()[0];
+            float scalar = (aMatrix.getDataSize() == 1)?aMatrix.getData()[0]:aOther.getData()[0];
             auto matrix = (aMatrix.getDataSize() == 1)?aOther:aMatrix;
-            double* output = malloc(matrix.getDataSize());
-            vdFmaxI(matrix.getDataSize(),matrix.getData(),1,&scalar,0,output,1);
+            float* output = malloc(matrix.getDataSize());
+            vsFmaxI(matrix.getDataSize(),matrix.getData(),1,&scalar,0,output,1);
             return Matrix::New(output,matrix);
     }
     else if (aMatrix.getDimSize(1) == 1 || aOther.getDimSize(0) == 1) {
         if (aMatrix.getDimSize(1) == 1){
-            double* output = malloc(aOther.getDataSize());
+            float* output = malloc(aOther.getDataSize());
             for (int i = 0; i < aOther.getDimSize(1); ++i) {
-                vdFmaxI(aMatrix.getDataSize(), aMatrix.getData(), 1, aOther.getData() + aOther.getDimSize(0) * i, 1,
+                vsFmaxI(aMatrix.getDataSize(), aMatrix.getData(), 1, aOther.getData() + aOther.getDimSize(0) * i, 1,
                         output + aOther.getDimSize(0) * i, 1);
             }
             return Matrix::New(output,aOther);
         }
         else{
-            double* output = malloc(aMatrix.getDataSize());
+            float* output = malloc(aMatrix.getDataSize());
             for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
-                vdFmaxI(aOther.getDataSize(), aOther.getData(), 1, aMatrix.getData() +  i, aMatrix.getDimSize(0),
+                vsFmaxI(aOther.getDataSize(), aOther.getData(), 1, aMatrix.getData() +  i, aMatrix.getDimSize(0),
                         output +  i, aOther.getDimSize(0));
             }
             return Matrix::New(output,aMatrix);
@@ -367,8 +367,8 @@ Matrix Aurora::max(const Matrix &aMatrix, const Matrix &aOther) {
     }
 }
 
-Matrix Aurora::max(const Matrix &aMatrix, const double aValue){
-    double *output = malloc(1);
+Matrix Aurora::max(const Matrix &aMatrix, const float aValue){
+    float *output = malloc(1);
     output[0] = aValue;
     return max(aMatrix,Matrix::New(output, 1,1,1));
 }
@@ -388,27 +388,27 @@ Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
         {
             case All:
             {
-                Eigen::Map<Eigen::VectorXcd> srcV((std::complex<double>*)aMatrix.getData(),aMatrix.getDataSize());
-                std::complex<double>* ret = (std::complex<double>*)malloc(1,true);
+                Eigen::Map<Eigen::VectorXcf> srcV((std::complex<float>*)aMatrix.getData(),aMatrix.getDataSize());
+                std::complex<float>* ret = (std::complex<float>*)malloc(1,true);
                 ret[0] = srcV.array().sum();
-                return Matrix::New((double*)ret,1,1,1,Complex);
+                return Matrix::New((float*)ret,1,1,1,Complex);
             }
             case Row:
             {
-                Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                std::complex<double> * ret = (std::complex<double>*)malloc(aMatrix.getDimSize(0),true);
-                Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::MatrixXcf> srcM((std::complex<float>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                std::complex<float> * ret = (std::complex<float>*)malloc(aMatrix.getDimSize(0),true);
+                Eigen::Map<Eigen::VectorXcf> retV(ret,aMatrix.getDimSize(0));
                 retV = srcM.rowwise().sum();
-                return Matrix::New((double*)ret,aMatrix.getDimSize(0),1,1,Complex);
+                return Matrix::New((float*)ret,aMatrix.getDimSize(0),1,1,Complex);
             }
             case Column:
             default:
             {
-                Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                std::complex<double>* ret = (std::complex<double>*)malloc(aMatrix.getDimSize(1),true);
-                Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXcf> srcM((std::complex<float>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                std::complex<float>* ret = (std::complex<float>*)malloc(aMatrix.getDimSize(1),true);
+                Eigen::Map<Eigen::VectorXcf> retV(ret,aMatrix.getDimSize(1));
                 retV = srcM.colwise().sum();
-                return Matrix::New((double*)ret,1,aMatrix.getDimSize(1),1,Complex);
+                return Matrix::New((float*)ret,1,aMatrix.getDimSize(1),1,Complex);
             }
         }
     }
@@ -417,25 +417,25 @@ Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
         {
             case All:
             {
-                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
-                double * ret = malloc(1);
+                Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
+                float * ret = malloc(1);
                 ret[0] = srcV.array().sum();
                 return Matrix::New(ret,1);
             }
             case Row:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * ret = malloc(aMatrix.getDimSize(0));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * ret = malloc(aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(0));
                 retV = srcM.rowwise().sum();
                 return Matrix::New(ret,aMatrix.getDimSize(0),1);
             }
             case Column:
             default:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * ret = malloc(aMatrix.getDimSize(1));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * ret = malloc(aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(1));
                 retV = srcM.colwise().sum();
                 return Matrix::New(ret,1,aMatrix.getDimSize(1));
             }
@@ -459,25 +459,25 @@ Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIn
         {
             case All:
             {
-                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
-                double * ret = malloc(64);
+                Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
+                float * ret = malloc(64);
                 ret[0] = srcV.array().mean();
                 return Matrix::New(ret,1);
             }
             case Row:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * ret = malloc(aMatrix.getDimSize(0));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * ret = malloc(aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(0));
                 retV = srcM.rowwise().mean();
                 return Matrix::New(ret,aMatrix.getDimSize(0),1);
             }
             case Column:
             default:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * ret = malloc(aMatrix.getDimSize(1));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * ret = malloc(aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(1));
                 retV = srcM.colwise().mean();
                 return Matrix::New(ret,1,aMatrix.getDimSize(1));
             }
@@ -488,40 +488,40 @@ Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIn
         {
             case All:
             {
-                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
-                double * retVd = malloc(aMatrix.getDataSize());
-                Eigen::Map<Eigen::VectorXd> retV(retVd,aMatrix.getDataSize());
-                Eigen::VectorXd ones = Eigen::VectorXd(aMatrix.getDataSize());
+                Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData(),aMatrix.getDataSize());
+                float * retvs = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::VectorXf> retV(retvs,aMatrix.getDataSize());
+                Eigen::VectorXf ones = Eigen::VectorXf(aMatrix.getDataSize());
                 ones.setConstant(1.0);
                 retV = srcV.array().isNaN().select(0.0,ones);
                 int count = retV.sum();
                 if (count == 0){
-                    free(retVd);
-                    double *ret = malloc(1);
+                    free(retvs);
+                    float *ret = malloc(1);
                     ret[0]=0;
                     return Matrix::New(ret,1);
                 }
                 else {
-                    double *ret = malloc(1);
+                    float *ret = malloc(1);
                     retV = srcV.array().isNaN().select(0.0,srcV);
                     ret[0] = retV.sum() / count;
-                    free(retVd);
+                    free(retvs);
                     return Matrix::New(ret, 1);
                 }
             }
             case Row:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * retMd = malloc(aMatrix.getDataSize());
-                Eigen::Map<Eigen::MatrixXd> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                Eigen::MatrixXd zeros = Eigen::MatrixXd(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * retMd = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::MatrixXf> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                Eigen::MatrixXf zeros = Eigen::MatrixXf(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
                 zeros.setConstant(0.0);
                 retM = srcM.array().isNaN().select(zeros,1.0);
-                Eigen::VectorXd countM = retM.rowwise().sum();
+                Eigen::VectorXf countM = retM.rowwise().sum();
                 countM = (countM.array()==0.0).select(1.0,countM);
                 retM = srcM.array().isNaN().select(0.0,srcM);
-                double * ret = malloc(aMatrix.getDimSize(0));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                float * ret = malloc(aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(0));
                 retV = retM.rowwise().sum();
                 retV =retV.array()/countM.array();
                 free(retMd);
@@ -530,17 +530,17 @@ Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIn
             case Column:
             default:
             {
-                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                double * retMd = malloc(aMatrix.getDataSize());
-                Eigen::Map<Eigen::MatrixXd> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-                Eigen::MatrixXd zeros = Eigen::MatrixXd(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                float * retMd = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::MatrixXf> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                Eigen::MatrixXf zeros = Eigen::MatrixXf(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
                 zeros.setConstant(0.0);
                 retM = srcM.array().isNaN().select(zeros,1.0);
-                Eigen::VectorXd countM = retM.colwise().sum();
+                Eigen::VectorXf countM = retM.colwise().sum();
                 countM = (countM.array()==0).select(1.0,countM);
                 retM = srcM.array().isNaN().select(0.0,srcM);
-                double * ret = malloc(aMatrix.getDimSize(1));
-                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
+                float * ret = malloc(aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::VectorXf> retV(ret,aMatrix.getDimSize(1));
                 retV = retM.colwise().sum();
                 retV = retV.array()/countM.array();
                 free(retMd);
@@ -576,20 +576,20 @@ Matrix Aurora::sort(Matrix &&aMatrix, FunctionDirection direction) {
         if (aMatrix.getDimSize(0)>=100000){
             #pragma omp parallel for
             for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
-                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
                 std::sort(srcV.array().begin(),srcV.array().end());
             }
         }
         else
         {
             for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
-                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXf> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
                 std::sort(srcV.array().begin(),srcV.array().end());
             }
         }
     }
     else if(direction == Row){
-        Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+        Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
         if (aMatrix.getDimSize(1)>=100000){
             #pragma omp parallel for
             for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
@@ -618,8 +618,8 @@ Matrix Aurora::sortrows(const Matrix &aMatrix, Matrix* indexMatrix) {
     }
     auto result = aMatrix.deepCopy();
     int rows = aMatrix.getDimSize(0);
-    std::vector<std::pair<double,int>> col;
-    std::vector<std::pair<double,int>>::iterator last;
+    std::vector<std::pair<float,int>> col;
+    std::vector<std::pair<float,int>>::iterator last;
     for (size_t j = 0; j < rows; j++)
     {
         col.push_back(std::make_pair(aMatrix[j], j));
@@ -658,11 +658,11 @@ Matrix Aurora::sortrows(const Matrix &aMatrix, Matrix* indexMatrix) {
         //未发现不同值 break 循环
         if(!sameFlag) break;
         //按照新一列刷新数组值
-        std::for_each(col.begin(), col.end() , [&aMatrix,i,rows](std::pair<double,int>& a){return a.first=aMatrix[a.second+i*rows];});
+        std::for_each(col.begin(), col.end() , [&aMatrix,i,rows](std::pair<float,int>& a){return a.first=aMatrix[a.second+i*rows];});
     }
     int i=0;
     (*indexMatrix) = zeros(aMatrix.getDimSize(0),1);
-    std::for_each(col.begin(),col.end(), [&aMatrix,&result,&i,indexMatrix](std::pair<double,int>& a){ 
+    std::for_each(col.begin(),col.end(), [&aMatrix,&result,&i,indexMatrix](std::pair<float,int>& a){ 
         result(i,$) = aMatrix(a.second,$);
         (*indexMatrix)[i] = a.second;
         i++;
@@ -682,10 +682,10 @@ Matrix Aurora::median(const Matrix &aMatrix) {
     Matrix sorted = horVector?sort(aMatrix,Row):sort(aMatrix);
     int rows = horVector?sorted.getDimSize(1):sorted.getDimSize(0);
     int cols = horVector?sorted.getDimSize(0):sorted.getDimSize(1);
-    Eigen::Map<Eigen::MatrixXd> srcM(sorted.getData(),rows,cols);
+    Eigen::Map<Eigen::MatrixXf> srcM(sorted.getData(),rows,cols);
     bool flag = rows % 2 == 1;
-    double* ret = malloc(cols);
-    Eigen::Map<Eigen::VectorXd> retV(ret,cols);
+    float* ret = malloc(cols);
+    Eigen::Map<Eigen::VectorXf> retV(ret,cols);
     if (flag) {
         retV = srcM.row(rows/2);
         return Matrix::New(ret,1,cols);
@@ -696,20 +696,20 @@ Matrix Aurora::median(const Matrix &aMatrix) {
 }
 
 Matrix Aurora::fft(const Matrix &aMatrix, long aFFTSize) {
-    double *output = nullptr;
+    float *output = nullptr;
     mkl_free_buffers();
     MKL_LONG rowSize = (aFFTSize>0)?aFFTSize:aMatrix.getDimSize(0);
     //实际需要copy赋值的非0值
     MKL_LONG needCopySize = (rowSize<aMatrix.getDimSize(0))?rowSize:aMatrix.getDimSize(0);
     MKL_LONG bufferSize = rowSize*aMatrix.getDimSize(1);
     output = malloc(bufferSize, true);
-    double zero = 0.0;
+    float zero = 0.0;
     //先全部置为0
-    cblas_dcopy(bufferSize*2, &zero, 0, output, 1);
+    cblas_scopy(bufferSize*2, &zero, 0, output, 1);
     if (!aMatrix.isComplex()) {
         //按列copy原值
         for (int i = 0 ; i < aMatrix.getDimSize(1); ++i) {
-            cblas_dcopy(needCopySize, aMatrix.getData()+i*aMatrix.getDimSize(0), 1, output+i*rowSize*2, 2);
+            cblas_scopy(needCopySize, aMatrix.getData()+i*aMatrix.getDimSize(0), 1, output+i*rowSize*2, 2);
         }
     } else {
         //按列copy原值
@@ -721,8 +721,8 @@ Matrix Aurora::fft(const Matrix &aMatrix, long aFFTSize) {
     DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
     MKL_LONG status;
 
-    //创建 Descriptor, 精度 double , 输入类型实数, 维度1
-    status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, rowSize);
+    //创建 Descriptor, 精度 float , 输入类型实数, 维度1
+    status = DftiCreateDescriptor(&my_desc_handle, DFTI_SINGLE, DFTI_COMPLEX, 1, rowSize);
     if (status != DFTI_NO_ERROR) goto error;
 
     //通过 setValue 配置Descriptor
@@ -769,20 +769,20 @@ Matrix Aurora::ifft(const Matrix &aMatrix, long aFFTSize ) {
     // mkl_free_buffers();
     // auto output = malloc(aMatrix.getDataSize(),true);
     MKL_LONG status;
-    double *output = nullptr;
+    float *output = nullptr;
     mkl_free_buffers();
     MKL_LONG rowSize = (aFFTSize>0)?aFFTSize:aMatrix.getDimSize(0);
     //实际需要copy赋值的非0值
     MKL_LONG needCopySize = (rowSize<aMatrix.getDimSize(0))?rowSize:aMatrix.getDimSize(0);
     MKL_LONG bufferSize = rowSize*aMatrix.getDimSize(1);
     output = malloc(bufferSize, true);
-    double zero = 0.0;
+    float zero = 0.0;
     //先全部置为0
-    cblas_dcopy(bufferSize*2, &zero, 0, output, 1);
+    cblas_scopy(bufferSize*2, &zero, 0, output, 1);
     if (!aMatrix.isComplex()) {
         //按列copy原值
         for (int i = 0 ; i < aMatrix.getDimSize(1); ++i) {
-            cblas_dcopy(needCopySize, aMatrix.getData()+i*aMatrix.getDimSize(0), 1, output+i*rowSize*2, 2);
+            cblas_scopy(needCopySize, aMatrix.getData()+i*aMatrix.getDimSize(0), 1, output+i*rowSize*2, 2);
         }
     } else {
         //按列copy原值
@@ -790,9 +790,9 @@ Matrix Aurora::ifft(const Matrix &aMatrix, long aFFTSize ) {
             cblas_zcopy(needCopySize, aMatrix.getData()+i*aMatrix.getDimSize(0)*2, 1, output+i*rowSize*2, 1);
         }
     }
-    //创建 Descriptor, 精度 double , 输入类型实数, 维度1
+    //创建 Descriptor, 精度 float , 输入类型实数, 维度1
     int size = aMatrix.getDimSize(0);
-    status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, rowSize);
+    status = DftiCreateDescriptor(&my_desc_handle, DFTI_SINGLE, DFTI_COMPLEX, 1, rowSize);
     if (status != DFTI_NO_ERROR) goto error;
     //通过 setValue 配置Descriptor
     //使用单独的输出数据缓存
@@ -843,17 +843,17 @@ Matrix Aurora::ifft_symmetric(const Matrix &aMatrix,long length)
         return Matrix();
     }
     int matrixLength = aMatrix.getDataSize();
-    int resultHalfLength = (int)std::ceil(((double)length*0.5));
+    int resultHalfLength = (int)std::ceil(((float)length*0.5));
     int copyLength = resultHalfLength<matrixLength?resultHalfLength:matrixLength;
-    double* calcData = malloc(length,true);
-    double zero = 0.0;
+    float* calcData = malloc(length,true);
+    float zero = 0.0;
     //所有数据统一置0
-    cblas_dcopy(length*2,&zero,0,calcData,1);
+    cblas_scopy(length*2,&zero,0,calcData,1);
     //copy前半段数据
     cblas_zcopy(copyLength,aMatrix.getData(),1,calcData,1);
     //copy后半段数据,跳过index 0的值,并设置虚部共轭
-    vdAddI(copyLength-1,&zero,0,(aMatrix.getData()+2),2,(calcData+(length-1)*2),-2);
-    vdSubI(copyLength-1,&zero,0,(aMatrix.getData()+2+1),2,(calcData+(length-1)*2+1),-2);
+    vsAddI(copyLength-1,&zero,0,(aMatrix.getData()+2),2,(calcData+(length-1)*2),-2);
+    vsSubI(copyLength-1,&zero,0,(aMatrix.getData()+2+1),2,(calcData+(length-1)*2+1),-2);
 
     return real(ifft(Matrix::New(calcData,length,1,1,Complex)));
 }
@@ -862,13 +862,13 @@ void Aurora::fftshift(Matrix &aMatrix){
     int backwardLength = aMatrix.getDimSize(0)/2;
     if (aMatrix.getDimSize(0)%2!=0)++backwardLength;
     int forwardLength = aMatrix.getDimSize(0) - backwardLength;
-    double* buffer = malloc(forwardLength,true);
+    float* buffer = malloc(forwardLength,true);
     int valueStep = aMatrix.isComplex()?2:1;
     for (int i = 0; i<aMatrix.getDimSize(1); ++i) {
-        double* dataPtr = aMatrix.getData()+aMatrix.getDimSize(0)*i*valueStep;
-        cblas_dcopy(forwardLength*valueStep, dataPtr+backwardLength*valueStep, 1, buffer, 1);
-        cblas_dcopy(backwardLength*valueStep, dataPtr, 1, dataPtr+forwardLength*valueStep, 1);
-        cblas_dcopy(forwardLength*valueStep, buffer, 1, dataPtr, 1);
+        float* dataPtr = aMatrix.getData()+aMatrix.getDimSize(0)*i*valueStep;
+        cblas_scopy(forwardLength*valueStep, dataPtr+backwardLength*valueStep, 1, buffer, 1);
+        cblas_scopy(backwardLength*valueStep, dataPtr, 1, dataPtr+forwardLength*valueStep, 1);
+        cblas_scopy(forwardLength*valueStep, buffer, 1, dataPtr, 1);
     }
     Aurora::free(buffer);
 }
@@ -877,13 +877,13 @@ void Aurora::ifftshift(Matrix &aMatrix){
     int forwardLength= aMatrix.getDimSize(0)/2;
     if (aMatrix.getDimSize(0)%2!=0)++forwardLength;
     int backwardLength = aMatrix.getDimSize(0) - forwardLength;
-    double* buffer = malloc(forwardLength,true);
+    float* buffer = malloc(forwardLength,true);
     int valueStep = aMatrix.isComplex()?2:1;
     for (int i = 0; i<aMatrix.getDimSize(1); ++i) {
-        double* dataPtr = aMatrix.getData()+aMatrix.getDimSize(0)*i*valueStep;
-        cblas_dcopy(forwardLength*valueStep, dataPtr+backwardLength*valueStep, 1, buffer, 1);
-        cblas_dcopy(backwardLength*valueStep, dataPtr, 1, dataPtr+forwardLength*valueStep, 1);
-        cblas_dcopy(forwardLength*valueStep, buffer, 1, dataPtr, 1);
+        float* dataPtr = aMatrix.getData()+aMatrix.getDimSize(0)*i*valueStep;
+        cblas_scopy(forwardLength*valueStep, dataPtr+backwardLength*valueStep, 1, buffer, 1);
+        cblas_scopy(backwardLength*valueStep, dataPtr, 1, dataPtr+forwardLength*valueStep, 1);
+        cblas_scopy(forwardLength*valueStep, buffer, 1, dataPtr, 1);
     }
     Aurora::free(buffer);
 }
@@ -891,16 +891,16 @@ void Aurora::ifftshift(Matrix &aMatrix){
 Matrix Aurora::hilbert(const Matrix &aMatrix) {
     auto x = fft(aMatrix);
     auto h = malloc(aMatrix.getDimSize(0));
-    auto two = 2.0;
-    auto zero = 0.0;
-    cblas_dcopy(aMatrix.getDimSize(0), &zero, 0, h, 1);
-    cblas_dcopy(aMatrix.getDimSize(0) / 2, &two, 0, h, 1);
+    float two = 2.0;
+    float zero = 0.0;
+    cblas_scopy(aMatrix.getDimSize(0), &zero, 0, h, 1);
+    cblas_scopy(aMatrix.getDimSize(0) / 2, &two, 0, h, 1);
     h[aMatrix.getDimSize(0) / 2] = ((aMatrix.getDimSize(0) << 31) >> 31) ? 2.0 : 1.0;
     h[0] = 1.0;
     for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
-        auto p = (double *)(x.getData() + aMatrix.getDimSize(0)* i*2);
-        vdMulI(aMatrix.getDimSize(0), p, 2, h, 1, p, 2);
-        vdMulI(aMatrix.getDimSize(0), p + 1, 2, h, 1, p + 1, 2);
+        auto p = (float *)(x.getData() + aMatrix.getDimSize(0)* i*2);
+        vsMulI(aMatrix.getDimSize(0), p, 2, h, 1, p, 2);
+        vsMulI(aMatrix.getDimSize(0), p + 1, 2, h, 1, p + 1, 2);
     }
     auto result = ifft( x);
     free(h);
@@ -915,16 +915,16 @@ Matrix Aurora::prod(const Matrix &aMatrix) {
     int row = aMatrix.getDimSize(0)==1?aMatrix.getDimSize(1):aMatrix.getDimSize(0);
     int col = aMatrix.getDimSize(0)==1?1:aMatrix.getDimSize(1);
     if (aMatrix.isComplex()){
-        Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),row,col);
+        Eigen::Map<Eigen::MatrixXcf> srcM((std::complex<float>*)aMatrix.getData(),row,col);
         auto ret = malloc(col,true);
-        Eigen::Map<Eigen::VectorXcd> retV((std::complex<double>*)ret,col);
+        Eigen::Map<Eigen::VectorXcf> retV((std::complex<float>*)ret,col);
         retV = srcM.colwise().prod();
         return Matrix::New(ret,1,col,1,Complex);
     }
     else{
-        Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),row,col);
+        Eigen::Map<Eigen::MatrixXf> srcM(aMatrix.getData(),row,col);
         auto ret = malloc(col);
-        Eigen::Map<Eigen::VectorXd> retV(ret,col);
+        Eigen::Map<Eigen::VectorXf> retV(ret,col);
         retV = srcM.colwise().prod();
         return Matrix::New(ret,1,col);
     }
@@ -951,7 +951,7 @@ Matrix Aurora::dot(const Matrix &aMatrix,const Matrix& aOther,FunctionDirection
         {
             auto ret = malloc(aMatrix.getDimSize(1));
             for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
-                ret[i]=cblas_ddot(aMatrix.getDimSize(0),aMatrix.getData()+i*aMatrix.getDimSize(0),1,
+                ret[i]=cblas_sdot(aMatrix.getDimSize(0),aMatrix.getData()+i*aMatrix.getDimSize(0),1,
                                 aOther.getData()+i*aMatrix.getDimSize(0),1);
             }
             return Matrix::New(ret,1,aMatrix.getDimSize(1),1);
@@ -959,7 +959,7 @@ Matrix Aurora::dot(const Matrix &aMatrix,const Matrix& aOther,FunctionDirection
         else{
             auto ret = malloc(aMatrix.getDimSize(0));
             for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
-                ret[i] = cblas_ddot(aMatrix.getDimSize(1),aMatrix.getData()+i,aMatrix.getDimSize(0),
+                ret[i] = cblas_sdot(aMatrix.getDimSize(1),aMatrix.getData()+i,aMatrix.getDimSize(0),
                                 aOther.getData()+i,aMatrix.getDimSize(0));
             }
             return Matrix::New(ret,aMatrix.getDimSize(1),1,1);
@@ -979,12 +979,12 @@ Matrix Aurora::sub2ind(const Matrix &aVMatrixSize, std::initializer_list<Matrix>
         return Matrix();
     }
     size_t returnVectorSize = aSliceIdxsList.begin()->getDataSize();
-    double *strides =new double[aVMatrixSize.getDataSize()+1];
+    float *strides =new float[aVMatrixSize.getDataSize()+1];
     strides[0] = 1;
     for (int i = 1; i<=aVMatrixSize.getDataSize(); ++i) {
         strides[i] = strides[i-1]*aVMatrixSize.getData()[i-1];
     } 
-    double* output = Aurora::malloc(returnVectorSize);
+    float* output = Aurora::malloc(returnVectorSize);
     for (size_t i = 0; i<returnVectorSize; ++i) {
         size_t j = 0;
         std::for_each(aSliceIdxsList.begin(), aSliceIdxsList.end(), [strides,output,i,&aVMatrixSize,&j](const Matrix& matrix){
diff --git a/src/Function2D.h b/src/Function2D.h
index 4255296..f056a3b 100644
--- a/src/Function2D.h
+++ b/src/Function2D.h
@@ -15,7 +15,7 @@ namespace Aurora
         Row,
         All
     };
-    double immse(const Matrix &aImageA, const Matrix &aImageB);
+    float immse(const Matrix &aImageA, const Matrix &aImageB);
     Matrix inv(const Matrix &aMatrix);
     Matrix inv(Matrix &&aMatrix);
     Matrix interp2(const Matrix &aX, const Matrix &aY, const Matrix &aV, const Matrix &aX1, const Matrix &aY1, InterpnMethod aMethod);
@@ -67,7 +67,7 @@ namespace Aurora
      * @param aOther 目标矩阵2
      * @return 最大值矩阵
      */
-    Matrix max(const Matrix &aMatrix, const double aValue);
+    Matrix max(const Matrix &aMatrix, const float aValue);
 
     /**
      * 求矩阵和，可按行、列、单元, 目前不支持三维
diff --git a/src/Function3D.cpp b/src/Function3D.cpp
index 7619c50..01d08af 100644
--- a/src/Function3D.cpp
+++ b/src/Function3D.cpp
@@ -30,10 +30,10 @@ Matrix Aurora::interp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, con
     {
       return Matrix();
     }
-    std::shared_ptr<double> resultData = std::shared_ptr<double>(new double[nx1], std::default_delete<double[]>());
+    std::shared_ptr<float> resultData = std::shared_ptr<float>(new float[nx1], std::default_delete<float[]>());
     for (int i = 0; i < nx1; ++i)
     {
-        std::shared_ptr<double> zResultData = std::shared_ptr<double>(new double[zAxisNum], std::default_delete<double[]>());
+        std::shared_ptr<float> zResultData = std::shared_ptr<float>(new float[zAxisNum], std::default_delete<float[]>());
         for(int j =0; j < zAxisNum; ++j)
         {
             zResultData.get()[j] = interp2(aX,aY,aV($,$,j).toMatrix(),aX1(i).toMatrix(),aY1(i).toMatrix(),aMethod).getData()[0];
@@ -59,9 +59,9 @@ Matrix Aurora::ones(int aRow, int aColumn, int aSlice) {
     int colSize =  aColumn;
     int sliceSize = aSlice == 0 ? 1 : aSlice;
     size_t arraySize  = rowSize * colSize* sliceSize;
-    double* data = malloc(arraySize);
-    double one = 1.0;
-    cblas_dcopy(arraySize,&one,0,data,1);
+    float* data = malloc(arraySize);
+    float one = 1.0;
+    cblas_scopy(arraySize,&one,0,data,1);
     return Matrix::New(data,rowSize,colSize,aSlice);
 }
 
@@ -79,9 +79,9 @@ Matrix Aurora::zeros(int aRow, int aColumn, int aSlice) {
     int colSize =  aColumn;
     int sliceSize = aSlice == 0 ? 1 : aSlice;
     size_t arraySize  = rowSize * colSize* sliceSize;
-    double* data = malloc(arraySize);
-    double zero = 0.0;
-    cblas_dcopy(arraySize,&zero,0,data,1);
+    float* data = malloc(arraySize);
+    float zero = 0.0;
+    cblas_scopy(arraySize,&zero,0,data,1);
     return Matrix::New(data,rowSize,colSize,sliceSize);
 }
 
@@ -92,19 +92,19 @@ Matrix Aurora::zeros(int aSquareRow) {
 Matrix Aurora::size(const Matrix &aMatrix)
 {
     if (aMatrix.isScalar()){
-        double * output = Aurora::malloc(1);
+        float * output = Aurora::malloc(1);
         output[0]=1;
         return Matrix::New(output,1,1,1);
     }
     else if (aMatrix.isVector()){
-        double * output = Aurora::malloc(2);
+        float * output = Aurora::malloc(2);
         output[0]=aMatrix.getDimSize(0);
         output[1]=aMatrix.getDimSize(1);
         return Matrix::New(output,2,1,1);
     }
     //3D
     else if (aMatrix.getDimSize(2)>1){
-        double * output = Aurora::malloc(3);
+        float * output = Aurora::malloc(3);
         output[0]=aMatrix.getDimSize(0);
         output[1]=aMatrix.getDimSize(1);
         output[2]=aMatrix.getDimSize(2);
@@ -112,7 +112,7 @@ Matrix Aurora::size(const Matrix &aMatrix)
     }
     //2D matrix
     else{
-        double * output = Aurora::malloc(2);
+        float * output = Aurora::malloc(2);
         output[0]=aMatrix.getDimSize(0);
         output[1]=aMatrix.getDimSize(1);
         return Matrix::New(output,2,1,1);
@@ -124,7 +124,7 @@ int Aurora::size(const Matrix &aMatrix,int dims)
     return aMatrix.getDimSize(dims-1);
 }
 
-Matrix Aurora::meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod, double aExtrapval)
+Matrix Aurora::meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod, float aExtrapval)
 {
     std::vector<Matrix> zTemps;
     for(int rIndex=0; rIndex<aV.getDimSize(2); ++rIndex)
@@ -140,7 +140,7 @@ Matrix Aurora::meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix&
         std::vector<Matrix> xTemps;
         for(int yIndex=0; yIndex<yTemps[0].getDataSize(); ++yIndex)
         {
-            double* xTempData = new double[yTemps.size()];
+            float* xTempData = new float[yTemps.size()];
             for(int i=0;i<yTemps.size();++i)
             {
                 xTempData[i] =  yTemps[i][yIndex];
@@ -162,7 +162,7 @@ Matrix Aurora::meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix&
     std::vector<Matrix>  resultVector;
     for(int i=0; i<zTemps[0].getDataSize(); ++i)
     {
-        double* resultTemp = new double[zTemps.size()]; 
+        float* resultTemp = new float[zTemps.size()]; 
         for(int j=0; j<zTemps.size(); ++j)
         {
             resultTemp[j] = zTemps[j][i];
@@ -180,12 +180,12 @@ Matrix Aurora::meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix&
     result = transpose(result);
     result.forceReshape(aX1.getDataSize(), aY1.getDataSize(), aZ1.getDataSize());
 
-    double minXSOS = aX[0];
-    double maxXSOS = aX[aX.getDataSize() - 1];
-    double minYSOS = aY[0];
-    double maxYSOS = aY[aY.getDataSize() - 1];
-    double minZSOS = aZ[0];
-    double maxZSOS = aZ[aZ.getDataSize() - 1];
+    float minXSOS = aX[0];
+    float maxXSOS = aX[aX.getDataSize() - 1];
+    float minYSOS = aY[0];
+    float maxYSOS = aY[aY.getDataSize() - 1];
+    float minZSOS = aZ[0];
+    float maxZSOS = aZ[aZ.getDataSize() - 1];
     for(size_t zIndex=0; zIndex<aZ1.getDataSize(); ++zIndex)
     {
         for(size_t yIndex=0; yIndex<aY1.getDataSize(); ++yIndex)
diff --git a/src/Function3D.h b/src/Function3D.h
index 1a8a3f6..f57f734 100644
--- a/src/Function3D.h
+++ b/src/Function3D.h
@@ -39,7 +39,7 @@ namespace Aurora {
      */
     Matrix zeros(int aSquareRow);
     Matrix interp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod);
-    Matrix meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod, double aExtrapval);
+    Matrix meshgridInterp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod, float aExtrapval);
     Matrix interpn(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod);
 
     Matrix size(const Matrix &aMatrix);
diff --git a/src/MatlabReader.cpp b/src/MatlabReader.cpp
index 761c388..f6fba2c 100644
--- a/src/MatlabReader.cpp
+++ b/src/MatlabReader.cpp
@@ -1,13 +1,14 @@
 #include "MatlabReader.h"
 #include "Function.h"
+#include <matio.h>
 
 
-double* calculateArray(matvar_t* aMatVar, const char* aFieldName)
+float* calculateArray(matvar_t* aMatVar, const char* aFieldName)
 {
     matvar_t* matVar = Mat_VarGetStructFieldByName(aMatVar,aFieldName,0);
-    double* data = reinterpret_cast<double*>(matVar->data);
+    float* data = reinterpret_cast<float*>(matVar->data);
     unsigned long long dataLength = matVar->nbytes / static_cast<unsigned long long>(matVar->data_size);
-    double* result = new double[dataLength];
+    float* result = new float[dataLength];
     std::copy(data,data+dataLength,result);
     Mat_VarFree(matVar);
     return result;
@@ -68,7 +69,7 @@ Aurora::Matrix MatlabReader::read(const std::string& aFieldName)
     if(!isComplex)
     {
         size_t matrixSize =  matvar->nbytes / static_cast<unsigned long long>(matvar->data_size);
-        double* matrixData = new double[matrixSize];
+        float* matrixData = new float[matrixSize];
         switch(matvar->data_type )
         {
             case MAT_T_INT16:
@@ -86,9 +87,9 @@ Aurora::Matrix MatlabReader::read(const std::string& aFieldName)
                 std::copy((int*)matvar->data, (int*)matvar->data + matrixSize, matrixData);
                 break;
             }
-            case MAT_T_DOUBLE:
+            case MAT_T_SINGLE:
             {
-                std::copy((double*)matvar->data, (double*)matvar->data + matrixSize, matrixData);
+                std::copy((float*)matvar->data, (float*)matvar->data + matrixSize, matrixData);
                 break;
             }
         }
@@ -96,11 +97,11 @@ Aurora::Matrix MatlabReader::read(const std::string& aFieldName)
     }
     else
     {
-        mat_complex_split_t* complexDouble = reinterpret_cast<mat_complex_split_t*>(matvar->data);
+        mat_complex_split_t* complexfloat = reinterpret_cast<mat_complex_split_t*>(matvar->data);
         unsigned long long dataLength =  matvar->nbytes / static_cast<unsigned long long>(matvar->data_size);
-        double* data = Aurora::malloc(dataLength, true);
-        double* re = reinterpret_cast<double*>(complexDouble->Re);
-        double* im = reinterpret_cast<double*>(complexDouble->Im);
+        float* data = Aurora::malloc(dataLength, true);
+        float* re = reinterpret_cast<float*>(complexfloat->Re);
+        float* im = reinterpret_cast<float*>(complexfloat->Im);
         for ( unsigned long long i = 0; i < dataLength; i++)
         {
             data[2*i] = re[i];
@@ -138,7 +139,7 @@ std::vector<Aurora::Matrix> MatlabReader::read4d(const std::string& aFieldName)
     std::vector<Aurora::Matrix> result;
     if(!isComplex)
     {
-        double* matrix3d = reinterpret_cast<double*>(matvar->data);
+        float* matrix3d = reinterpret_cast<float*>(matvar->data);
         for(int i=0; i<dims4; ++i)
         {
             result.push_back(Aurora::Matrix::copyFromRawData(matrix3d,rows,columns,slices));
diff --git a/src/MatlabWriter.cpp b/src/MatlabWriter.cpp
index 49bb18d..f6ebc72 100644
--- a/src/MatlabWriter.cpp
+++ b/src/MatlabWriter.cpp
@@ -27,7 +27,7 @@ void MatlabWriter::write()
         std::string name = i->first;
         int dimsSize = matrix.getDims();
         size_t dims[3] ={(size_t)matrix.getDimSize(0),(size_t)matrix.getDimSize(1), (size_t)matrix.getDimSize(2)};
-        matvar_t* var =  Mat_VarCreate(name.c_str(),MAT_C_DOUBLE,MAT_T_DOUBLE,dimsSize,dims,matrix.getData(),0);
+        matvar_t* var =  Mat_VarCreate(name.c_str(),MAT_C_SINGLE,MAT_T_SINGLE,dimsSize,dims,matrix.getData(),0);
         Mat_VarWrite(mMat,var,MAT_COMPRESSION_NONE);
         Mat_VarFree(var);
     }
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 005e208..3932c45 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -17,19 +17,19 @@
 #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(*CalcFuncD)(const MKL_INT n, const float a[], const MKL_INT inca, const float b[],
+                             const MKL_INT incb, float 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);
+    typedef void(*CalcFuncZ)(const MKL_INT n, const MKL_Complex8 a[], const MKL_INT inca, const MKL_Complex8 b[],
+                             const MKL_INT incb, MKL_Complex8 r[], const MKL_INT incr);
 
-    inline Matrix operatorAxM(CalcFuncD aFunc, double aScalar, const Matrix &aMatrix) {
-        double *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
+    inline Matrix operatorAxM(CalcFuncD aFunc, float aScalar, const Matrix &aMatrix) {
+        float *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
         //复数时，+和-需要特别操作，只影响实部
-        if (aMatrix.getValueType() == Complex && (aFunc == &vdAddI || aFunc == &vdSubI)) {
+        if (aMatrix.getValueType() == Complex && (aFunc == &vsAddI || aFunc == &vsSubI)) {
             aFunc(aMatrix.getDataSize(),  &aScalar, 0,aMatrix.getData() , 2, output ,
                   2);
-            double zero = 0.0;
+            float zero = 0.0;
             aFunc(aMatrix.getDataSize(),  &zero, 0,aMatrix.getData()+1 , 2, output+1 ,
                   2);
         } else{
@@ -39,13 +39,13 @@ namespace Aurora{
                                    aMatrix.getValueType());
     }
 
-    inline Matrix operatorMxA(CalcFuncD aFunc, double aScalar, const Matrix &aMatrix) {
-        double *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
+    inline Matrix operatorMxA(CalcFuncD aFunc, float aScalar, const Matrix &aMatrix) {
+        float *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
         //复数时，+和-需要特别操作，只影响实部
-        if (aMatrix.getValueType() == Complex && (aFunc == &vdAddI || aFunc == &vdSubI)) {
+        if (aMatrix.getValueType() == Complex && (aFunc == &vsAddI || aFunc == &vsSubI)) {
             aFunc(aMatrix.getDataSize(), aMatrix.getData() , 2, &aScalar, 0, output ,
                   2);
-            double zero = 0.0;
+            float zero = 0.0;
             aFunc(aMatrix.getDataSize(), aMatrix.getData()+1 , 2, &zero, 0, output+1 ,
                   2);
         } else{
@@ -55,14 +55,14 @@ namespace Aurora{
                                    aMatrix.getValueType());
     }
 
-    inline Matrix &operatorAxM_RR(CalcFuncD aFunc, double aScalar, Aurora::Matrix &&aMatrix) {
+    inline Matrix &operatorAxM_RR(CalcFuncD aFunc, float aScalar, Aurora::Matrix &&aMatrix) {
         if (aMatrix.getValueType() == Complex) {
             //针对实部操作
             aFunc(aMatrix.getDataSize(),  &aScalar, 0,aMatrix.getData(), 2,
                   aMatrix.getData(),
                   2);
             //乘法除法需特别操作，影响虚部
-            if (aFunc == &vdDivI || aFunc == &vdMulI){
+            if (aFunc == &vsDivI || aFunc == &vsMulI){
                 aFunc(aMatrix.getDataSize(),  &aScalar, 0,aMatrix.getData() + 1, 2,
                   aMatrix.getData() + 1, 2);
             }
@@ -74,14 +74,14 @@ namespace Aurora{
         return aMatrix;
     }
 
-    inline Matrix &operatorMxA_RR(CalcFuncD aFunc, double aScalar, Aurora::Matrix &&aMatrix) {
+    inline Matrix &operatorMxA_RR(CalcFuncD aFunc, float aScalar, Aurora::Matrix &&aMatrix) {
         if (aMatrix.getValueType() == Complex) {
             //针对实部操作
             aFunc(aMatrix.getDataSize(), aMatrix.getData(), 2, &aScalar, 0,
                   aMatrix.getData(),
                   2);
             //乘法除法需特别操作，影响虚部
-            if (aFunc == &vdDivI || aFunc == &vdMulI){
+            if (aFunc == &vsDivI || aFunc == &vsMulI){
                 aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 2, &aScalar, 0,
                   aMatrix.getData() + 1, 2);
             }
@@ -96,73 +96,73 @@ namespace Aurora{
 
     inline void V_MxM_CN_Calc(
             CalcFuncD aFuncD,
-            const int size, double* xC,double* yN,double *output, int DimsStride) {
+            const int size, float* xC,float* yN,float *output, int DimsStride) {
         aFuncD(size, xC, DimsStride * 2, yN, 1, output, 2);
-        if (aFuncD == &vdDivI || aFuncD == &vdMulI){
+        if (aFuncD == &vsDivI || aFuncD == &vsMulI){
             aFuncD(size, xC + 1, DimsStride * 2, yN, 1, output + 1, 2);
         }
         else{
-            double zero = 0.0;
+            float zero = 0.0;
             aFuncD(size, xC+1 , DimsStride*2, &zero, 0, output + 1, 2);
         }
     }
 
-    inline double* _MxM_CN_Calc(
+    inline float* _MxM_CN_Calc(
             CalcFuncD aFuncD,
-            const int size, double* xC,double* yN, int dimsStride)
+            const int size, float* xC,float* yN, int dimsStride)
     {
-        double *output = malloc(size, true);
+        float *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* xN,double* yC,double *output, int DimsStride) {
+            const int size, float* xN,float* yC,float *output, int DimsStride) {
         aFuncD(size, xN, DimsStride, yC, 2, output, 2);
-        if (aFuncD == &vdDivI || aFuncD == &vdMulI){
+        if (aFuncD == &vsDivI || aFuncD == &vsMulI){
             aFuncD(size, xN , DimsStride, yC+ 1, 2, output + 1, 2);
         }
         else{
-            double zero = 0.0;
+            float zero = 0.0;
             aFuncD(size, yC+ 1, 2, &zero, 0, output + 1, 2);
         }
     }
 
-    inline double* _MxM_NC_Calc(
+    inline float* _MxM_NC_Calc(
             CalcFuncD aFuncD,
-            const int size, double* xN,double* yC, int dimsStride)
+            const int size, float* xN,float* yC, int dimsStride)
     {
-        double *output = malloc(size, true);
+        float *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) {
+            const int size, float* x,float* y,float* output, int DimsStride) {
         aFuncD(size, x, DimsStride, y, 1, output,1);
     }
 
-    inline double* _MxM_NN_Calc(
+    inline float* _MxM_NN_Calc(
             CalcFuncD aFuncD,
-            const int size, double* x,double* y, int DimsStride) {
-        double *output = malloc(size);
+            const int size, float* x,float* y, int DimsStride) {
+        float *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,
+            CalcFuncZ aFuncZ, const int size, float* x,float* y,float* output,
             int DimsStride) {
-        aFuncZ(size, (std::complex<double> *) x, DimsStride,
-               (std::complex<double> *) y, 1, (std::complex<double> *) output, 1);
+        aFuncZ(size, (std::complex<float> *) x, DimsStride,
+               (std::complex<float> *) y, 1, (std::complex<float> *) output, 1);
     }
 
-    inline double* _MxM_CC_Calc(
-            CalcFuncZ aFuncZ, const int size, double* x,double* y,
+    inline float* _MxM_CC_Calc(
+            CalcFuncZ aFuncZ, const int size, float* x,float* y,
             int DimsStride) {
-        double *output = malloc(size, true);
+        float *output = malloc(size, true);
         V_MxM_CC_Calc(aFuncZ, size, x, y, output, DimsStride);
         return output;
     }
@@ -172,7 +172,7 @@ namespace Aurora{
         // 2v2,1v1,3v3
         if (aMatrix.compareShape(aOther)) {
             int DimsStride = 1;
-            double *data = nullptr;
+            float *data = nullptr;
             if (aMatrix.getValueType() != aOther.getValueType()) {
                 if (aMatrix.getValueType() == Normal) {
                     data = _MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
@@ -219,12 +219,12 @@ namespace Aurora{
                                  Aurora::Matrix &&aOther,bool oppside = false) {
         if (aMatrix.compareShape(aOther)) {
             int DimsStride = 1;
-            double* X = oppside?aOther.getData():aMatrix.getData();
-            double* Y = oppside?aMatrix.getData():aOther.getData();
+            float* X = oppside?aOther.getData():aMatrix.getData();
+            float* Y = oppside?aMatrix.getData():aOther.getData();
             if (aMatrix.getValueType() != aOther.getValueType()) {
                 //aOther is not a complex matrix
                 if (aMatrix.getValueType() == Complex) {
-                    double *output = oppside?_MxM_NC_Calc(aFuncD,aMatrix.getDataSize(),  aOther.getData(),aMatrix.getData(),DimsStride)
+                    float *output = oppside?_MxM_NC_Calc(aFuncD,aMatrix.getDataSize(),  aOther.getData(),aMatrix.getData(),DimsStride)
                         :_MxM_CN_Calc(aFuncD,aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),DimsStride);
                     return Matrix::New(output, aMatrix);
                 }
@@ -275,7 +275,7 @@ namespace Aurora{
 
 
 namespace Aurora {
-    Matrix::Matrix(std::shared_ptr<double> aData, std::vector<int> aInfo, ValueType aValueType)
+    Matrix::Matrix(std::shared_ptr<float> aData, std::vector<int> aInfo, ValueType aValueType)
             : mValueType(aValueType)
             , mData(aData)
             , mInfo(aInfo)
@@ -316,7 +316,7 @@ namespace Aurora {
                 getDimSize(2) < 2);
     }
 
-    double Matrix::getScalar() const {
+    float Matrix::getScalar() const {
         if (isNull()) return 0.0;
         if (isNull()) return 0.0;
         return getData()[0];
@@ -342,7 +342,7 @@ namespace Aurora {
         return 2;
     }
 
-    double *Matrix::getData() const {
+    float *Matrix::getData() const {
         return mData.get();
     }
 
@@ -379,7 +379,7 @@ namespace Aurora {
         return true;
     }
 
-    Matrix Matrix::New(double *data, int rows, int cols, int slices, ValueType type) {
+    Matrix Matrix::New(float *data, int rows, int cols, int slices, ValueType type) {
         if (!data) return Matrix();
         std::vector<int> vector(3);
         vector[0]=rows;
@@ -391,7 +391,7 @@ namespace Aurora {
         return ret;
     }
 
-    Matrix Matrix::New(double *data, const Matrix &shapeMatrix) {
+    Matrix Matrix::New(float *data, const Matrix &shapeMatrix) {
         return New(data,
                    shapeMatrix.getDimSize(0),
                    shapeMatrix.getDimSize(1),
@@ -400,33 +400,33 @@ namespace Aurora {
     }
 
     Matrix Matrix::New(const Matrix &shapeMatrix) {
-        double *newBuffer = malloc(shapeMatrix.getDataSize(), shapeMatrix.getValueType());
+        float *newBuffer = malloc(shapeMatrix.getDataSize(), shapeMatrix.getValueType());
         return New(newBuffer, shapeMatrix);
     }
 
-    Matrix Matrix::fromRawData(double *data, int rows, int cols, int slices, ValueType type) {
+    Matrix Matrix::fromRawData(float *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);
+        Matrix ret({data, std::default_delete<float[]>()}, vector);
         if (type != Normal)ret.setValueType(type);
         return ret;
     }
 
-    Matrix Matrix::copyFromRawData(double *data, int rows, int cols, int slices, ValueType type) {
+    Matrix Matrix::copyFromRawData(float *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 = Aurora::malloc(size, type);
-        cblas_dcopy(size*type,data,1,newBuffer,1);
+        float *newBuffer = Aurora::malloc(size, type);
+        cblas_scopy(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);
+        float *newBuffer = malloc(getDataSize(), getValueType()==Complex);
+        cblas_scopy(getDataSize() * getValueType(),getData(),1,newBuffer,1);
         return New(newBuffer,
                    getDimSize(0),
                    getDimSize(1),
@@ -434,8 +434,8 @@ namespace Aurora {
                    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+(float aScalar) const { return operatorMxA(&vsAddI, aScalar, *this);}
+    Matrix operator+(float aScalar, const Matrix &matrix) {return matrix + aScalar;}
     Matrix Matrix::operator+(const Matrix &matrix) const {
         if (isScalar()){
             return getScalar()+matrix;
@@ -443,13 +443,13 @@ namespace Aurora {
         if (matrix.isScalar()){
             return (*this)+matrix.getScalar();
         }
-        return operatorMxM(vdAddI, vzAddI, *this, matrix);
+        return operatorMxM(vsAddI, vcAddI, *this, matrix);
     }
-    Matrix &operator+(double aScalar, Matrix &&matrix) {
-        return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator+(float aScalar, Matrix &&matrix) {
+        return operatorMxA_RR(&vsAddI,aScalar, std::forward<Matrix&&>(matrix));
     }
-    Matrix &operator+(Matrix &&matrix,double aScalar) {
-        return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator+(Matrix &&matrix,float aScalar) {
+        return operatorMxA_RR(&vsAddI,aScalar, std::forward<Matrix&&>(matrix));
     }
     Matrix Matrix::operator+(Matrix &&aMatrix) const {
         if (isScalar()){
@@ -458,7 +458,7 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)+aMatrix.getScalar();
         }
-        return operatorMxM_RR(&vdAddI,&vzAddI,*this,std::forward<Matrix&&>(aMatrix));
+        return operatorMxM_RR(&vsAddI,&vcAddI,*this,std::forward<Matrix&&>(aMatrix));
     }
 
     Matrix operator+(Matrix &&aMatrix, Matrix &aOther) {
@@ -468,11 +468,11 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return aOther+aMatrix.getScalar();
         }
-        return operatorMxM_RR(&vdAddI,&vzAddI,aOther,std::forward<Matrix&&>(aMatrix),true);
+        return operatorMxM_RR(&vsAddI,&vcAddI,aOther,std::forward<Matrix&&>(aMatrix),true);
     }
     //operation -
-    Matrix Matrix::operator-(double aScalar) const { return operatorMxA(&vdSubI, aScalar, *this);}
-    Matrix operator-(double aScalar, const Matrix &matrix) {return operatorAxM(&vdSubI, aScalar, matrix);}
+    Matrix Matrix::operator-(float aScalar) const { return operatorMxA(&vsSubI, aScalar, *this);}
+    Matrix operator-(float aScalar, const Matrix &matrix) {return operatorAxM(&vsSubI, aScalar, matrix);}
     Matrix Matrix::operator-(const Matrix &aMatrix) const {
         if (isScalar()){
             return getScalar()-aMatrix;
@@ -480,13 +480,13 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)-aMatrix.getScalar();
         }
-        return operatorMxM(vdSubI, vzSubI, *this, aMatrix);
+        return operatorMxM(vsSubI, vcSubI, *this, aMatrix);
     }
-    Matrix &operator-(double aScalar, Matrix &&matrix) {
-        return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator-(float aScalar, Matrix &&matrix) {
+        return operatorMxA_RR(&vsSubI,aScalar, std::forward<Matrix&&>(matrix));
     }
-    Matrix &operator-(Matrix &&matrix,double aScalar) {
-        return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator-(Matrix &&matrix,float aScalar) {
+        return operatorMxA_RR(&vsSubI,aScalar, std::forward<Matrix&&>(matrix));
     }
     Matrix Matrix::operator-(Matrix &&aMatrix) const {
         if (isScalar()){
@@ -495,7 +495,7 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)-aMatrix.getScalar();
         }
-        return operatorMxM_RR(&vdSubI,&vzSubI,*this,std::forward<Matrix&&>(aMatrix));
+        return operatorMxM_RR(&vsSubI,&vcSubI,*this,std::forward<Matrix&&>(aMatrix));
     }
 
     Matrix operator-(Matrix &&aMatrix, Matrix &aOther) {
@@ -505,11 +505,11 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return aMatrix.getScalar()-aOther;
         }
-        return operatorMxM_RR(&vdSubI,&vzSubI,aOther,std::forward<Matrix&&>(aMatrix),true);
+        return operatorMxM_RR(&vsSubI,&vcSubI,aOther,std::forward<Matrix&&>(aMatrix),true);
     }
     //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*(float aScalar) const { return operatorMxA(&vsMulI, aScalar, *this);}
+    Matrix operator*(float aScalar, const Matrix &matrix) {return matrix * aScalar;}
     Matrix Matrix::operator*(const Matrix &aMatrix) const {
         if (isScalar()){
             return getScalar()*aMatrix;
@@ -517,13 +517,13 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)*aMatrix.getScalar();
         }
-        return operatorMxM(vdMulI, vzMulI, *this, aMatrix);
+        return operatorMxM(vsMulI, vcMulI, *this, aMatrix);
     }
-    Matrix &operator*(double aScalar, Matrix &&matrix) {
-        return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator*(float aScalar, Matrix &&matrix) {
+        return operatorMxA_RR(&vsMulI,aScalar, std::forward<Matrix&&>(matrix));
     }
-    Matrix &operator*(Matrix &&aMatrix,double aScalar) {
-        return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(aMatrix));
+    Matrix &operator*(Matrix &&aMatrix,float aScalar) {
+        return operatorMxA_RR(&vsMulI,aScalar, std::forward<Matrix&&>(aMatrix));
     }
     Matrix Matrix::operator*(Matrix &&aMatrix) const {
         if (isScalar()){
@@ -532,7 +532,7 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)*aMatrix.getScalar();
         }
-        return operatorMxM_RR(&vdMulI,&vzMulI,*this,std::forward<Matrix&&>(aMatrix));
+        return operatorMxM_RR(&vsMulI,&vcMulI,*this,std::forward<Matrix&&>(aMatrix));
     }
 
     Matrix operator*(Matrix &&aMatrix, Matrix &aOther) {
@@ -542,11 +542,11 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return aMatrix.getScalar()*aOther;
         }
-        return operatorMxM_RR(&vdMulI,&vzMulI,aOther,std::forward<Matrix&&>(aMatrix),true);
+        return operatorMxM_RR(&vsMulI,&vcMulI,aOther,std::forward<Matrix&&>(aMatrix),true);
     }
     //operation /
-    Matrix Matrix::operator/(double aScalar) const { return operatorMxA(&vdDivI, aScalar, *this);}
-    Matrix operator/(double aScalar, const Matrix &matrix) {return operatorAxM(&vdDivI, aScalar, matrix);}
+    Matrix Matrix::operator/(float aScalar) const { return operatorMxA(&vsDivI, aScalar, *this);}
+    Matrix operator/(float aScalar, const Matrix &matrix) {return operatorAxM(&vsDivI, aScalar, matrix);}
     Matrix Matrix::operator/(const Matrix &aMatrix) const {
         if (isScalar()){
             return getScalar()/aMatrix;
@@ -554,13 +554,13 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)/aMatrix.getScalar();
         }
-        return operatorMxM(vdDivI, vzDivI, *this, aMatrix);
+        return operatorMxM(vsDivI, vcDivI, *this, aMatrix);
     }
-    Matrix &operator/(double aScalar, Matrix &&matrix) {
-        return operatorAxM_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator/(float aScalar, Matrix &&matrix) {
+        return operatorAxM_RR(&vsDivI,aScalar, std::forward<Matrix&&>(matrix));
     }
-    Matrix &operator/(Matrix &&matrix,double aScalar) {
-        return operatorMxA_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
+    Matrix &operator/(Matrix &&matrix,float aScalar) {
+        return operatorMxA_RR(&vsDivI,aScalar, std::forward<Matrix&&>(matrix));
     }
     Matrix Matrix::operator/(Matrix &&aMatrix) const {
         if (isScalar()){
@@ -569,7 +569,7 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return (*this)/aMatrix.getScalar();
         }
-        return operatorMxM_RR(&vdDivI,&vzDivI,*this,std::forward<Matrix&&>(aMatrix));
+        return operatorMxM_RR(&vsDivI,&vcDivI,*this,std::forward<Matrix&&>(aMatrix));
     }
 
     Matrix operator/(Matrix &&aMatrix, Matrix &aOther) {
@@ -579,17 +579,17 @@ namespace Aurora {
         if (aMatrix.isScalar()){
             return aMatrix.getScalar()/aOther;
         }
-        return operatorMxM_RR(&vdDivI,&vzDivI,aOther,std::forward<Matrix&&>(aMatrix),true);
+        return operatorMxM_RR(&vsDivI,&vcDivI,aOther,std::forward<Matrix&&>(aMatrix),true);
     }
     //operator ^ (pow)
-    Matrix Matrix::operator^(int times) const { return operatorMxA(&vdPowI, times, *this);}
+    Matrix Matrix::operator^(int times) const { return operatorMxA(&vsPowI, times, *this);}
 
     Matrix operator^( Matrix &&matrix,int times) {
-        return operatorMxA(&vdPowI, times, std::forward<Matrix&&>(matrix));
+        return operatorMxA(&vsPowI, times, std::forward<Matrix&&>(matrix));
     }
 
-    double& Matrix::operator[](size_t index) { return getData()[index];}
-    double Matrix::operator[](size_t index) const { return getData()[index];}
+    float& Matrix::operator[](size_t index) { return getData()[index];}
+    float Matrix::operator[](size_t index) const { return getData()[index];}
 
     Matrix Matrix::block(int aDim,int aBeginIndex, int aEndIndex) const{
         if(aDim>2 ){
@@ -610,7 +610,7 @@ namespace Aurora {
         }
         int dimLength = std::abs(aEndIndex-aBeginIndex)+1;
         int dataSize = getDataSize()/getDimSize(aDim)*dimLength;
-        double * dataOutput = malloc(dataSize,isComplex());
+        float * dataOutput = malloc(dataSize,isComplex());
         int colStride = getDimSize(0);
         int sliceStride = getDimSize(0)*getDimSize(1);
         switch (aDim) {
@@ -621,14 +621,14 @@ namespace Aurora {
                 {
                     for (size_t j = 0; j < getDimSize(1); j++)
                     {
-                        cblas_dcopy(
+                        cblas_scopy(
                             dimLength,
                             getData() + (aBeginIndex + j * colStride +
                                 i * sliceStride)*getValueType(),
                             getValueType(), dataOutput + (colStride2 * j + i * sliceStride2)*getValueType(),
                             getValueType());
                         if(isComplex()){
-                            cblas_dcopy(
+                            cblas_scopy(
                             dimLength,
                             getData()+1 + (aBeginIndex + j * colStride +
                                 i * sliceStride)*getValueType(),
@@ -645,14 +645,14 @@ namespace Aurora {
                 int copySize = sliceStride2;
                 for (size_t i = 0; i < getDimSize(2); i++)
                 {
-                    cblas_dcopy(copySize,
+                    cblas_scopy(copySize,
                                 getData() + getValueType()*(aBeginIndex * colStride +
                                     i * sliceStride),
                                  getValueType(), dataOutput +  getValueType()*(i * copySize), 
                                  getValueType());
                     if (isComplex())
                     {
-                        cblas_dcopy(copySize,
+                        cblas_scopy(copySize,
                                 getData()+1 + getValueType()*(aBeginIndex * colStride +
                                     i * sliceStride),
                                  getValueType(), dataOutput+1 +  getValueType()*(i * copySize), 
@@ -663,7 +663,7 @@ namespace Aurora {
             }
             case 2:{
                 int copySize = dimLength*sliceStride;
-                cblas_dcopy(copySize*getValueType(), getData() + aBeginIndex * sliceStride*getValueType() ,1, dataOutput, 1);
+                cblas_scopy(copySize*getValueType(), getData() + aBeginIndex * sliceStride*getValueType() ,1, dataOutput, 1);
                 return Matrix::New(dataOutput,getDimSize(0),getDimSize(1),dimLength,getValueType());
             }
         }
@@ -671,7 +671,7 @@ namespace Aurora {
         
     }
 
-    bool Matrix::setBlockValue(int aDim,int aBeginIndex, int aEndIndex, double value) {
+    bool Matrix::setBlockValue(int aDim,int aBeginIndex, int aEndIndex, float value) {
         if(aDim>2 ){
             std::cerr<<"setBlockValue only support 1D-3D data!"<<std::endl;
             return false;
@@ -702,7 +702,7 @@ namespace Aurora {
                 {
                     for (size_t j = 0; j < getDimSize(1); j++)
                     {
-                        cblas_dcopy(
+                        cblas_scopy(
                             dimLength,
                             &value,0,
                             getData() + aBeginIndex + j * colStride +
@@ -718,7 +718,7 @@ namespace Aurora {
                 int copySize = sliceStride2;
                 for (size_t i = 0; i < getDimSize(2); i++)
                 {
-                    cblas_dcopy(copySize,
+                    cblas_scopy(copySize,
                                 &value, 0 ,
                                 getData() + aBeginIndex * colStride +
                                     i * sliceStride, 1);
@@ -727,7 +727,7 @@ namespace Aurora {
             }
             case 2:{
                 int copySize = dimLength*sliceStride;
-                cblas_dcopy(copySize, &value, 0 ,getData() + aBeginIndex * sliceStride ,1);
+                cblas_scopy(copySize, &value, 0 ,getData() + aBeginIndex * sliceStride ,1);
                 return true;
             }
         }
@@ -735,7 +735,7 @@ namespace Aurora {
         
     }
 
-    bool Matrix::setBlockComplexValue(int aDim,int aBeginIndex, int aEndIndex, std::complex<double> value) {
+    bool Matrix::setBlockComplexValue(int aDim,int aBeginIndex, int aEndIndex, std::complex<float> value) {
         if(aDim>2 ){
             std::cerr<<"setBlockValue only support 1D-3D data!"<<std::endl;
             return false;
@@ -766,15 +766,15 @@ namespace Aurora {
                 {
                     for (size_t j = 0; j < getDimSize(1); j++)
                     {
-                        double real = value.real();
-                        double imag = value.imag();
-                        cblas_dcopy(
+                        float real = value.real();
+                        float imag = value.imag();
+                        cblas_scopy(
                             dimLength,
                             &real,0,
                             getData() + aBeginIndex*2 + j * colStride*2 +
                                 i * sliceStride*2,
                             2);
-                        cblas_dcopy(
+                        cblas_scopy(
                             dimLength,
                             &imag,0,
                             getData() +1 + aBeginIndex*2 + j * colStride*2 +
@@ -790,13 +790,13 @@ namespace Aurora {
                 int copySize = sliceStride2;
                 for (size_t i = 0; i < getDimSize(2); i++)
                 {
-                    double real = value.real();
-                    double imag = value.imag();
-                    cblas_dcopy(copySize,
+                    float real = value.real();
+                    float imag = value.imag();
+                    cblas_scopy(copySize,
                                 &real, 0 ,
                                 getData() + aBeginIndex * colStride*2 +
                                     i * sliceStride*2,2);
-                    cblas_dcopy(copySize,
+                    cblas_scopy(copySize,
                                 &imag, 0 ,
                                 getData() +1 + aBeginIndex * colStride*2 +
                                     i * sliceStride*2, 2);
@@ -805,10 +805,10 @@ namespace Aurora {
             }
             case 2:{
                 int copySize = dimLength*sliceStride;
-                double real = value.real();
-                double imag = value.imag();
-                cblas_dcopy(copySize, &real, 0 ,getData() + aBeginIndex * sliceStride*2 ,2);
-                cblas_dcopy(copySize, &imag, 0 ,getData() +1 + aBeginIndex * sliceStride*2 ,2);
+                float real = value.real();
+                float imag = value.imag();
+                cblas_scopy(copySize, &real, 0 ,getData() + aBeginIndex * sliceStride*2 ,2);
+                cblas_scopy(copySize, &imag, 0 ,getData() +1 + aBeginIndex * sliceStride*2 ,2);
                 return true;
             }
         }
@@ -882,7 +882,7 @@ namespace Aurora {
                 {
                     for (size_t j = 0; j < getDimSize(1); j++)
                     {
-                        cblas_dcopy(dimLength,
+                        cblas_scopy(dimLength,
                                     src.getData() +
                                         j * colStride2 * getValueType() +
                                         i * sliceStride2 * getValueType(),
@@ -892,7 +892,7 @@ namespace Aurora {
                                         i * sliceStride * getValueType(),
                                     getValueType());
                         if (isComplex()){
-                            cblas_dcopy(dimLength ,
+                            cblas_scopy(dimLength ,
                                     src.getData() + 1 +
                                         j * colStride2 * getValueType() +
                                         i * sliceStride2 * getValueType(),
@@ -913,13 +913,13 @@ namespace Aurora {
                 int copySize = sliceStride2;
                 for (size_t i = 0; i < getDimSize(2); i++)
                 {
-                    cblas_dcopy(copySize,
+                    cblas_scopy(copySize,
                                 src.getData()+i * sliceStride2*getValueType(),getValueType(),
                                 getData() + aBeginIndex * colStride*getValueType() +
                                     i * sliceStride*getValueType(), getValueType());
                     if (isComplex())
                     {
-                        cblas_dcopy(copySize,
+                        cblas_scopy(copySize,
                                 src.getData()+1 +i * sliceStride2*getValueType(),getValueType(),
                                 getData()+1 + aBeginIndex * colStride*getValueType() +
                                     i * sliceStride*getValueType(), getValueType());
@@ -930,9 +930,9 @@ namespace Aurora {
             //copy slice
             case 2:{
                 int copySize = dimLength*sliceStride;
-                cblas_dcopy(copySize,src.getData(),getValueType(),getData() + aBeginIndex * sliceStride*getValueType() ,getValueType());
+                cblas_scopy(copySize,src.getData(),getValueType(),getData() + aBeginIndex * sliceStride*getValueType() ,getValueType());
                 if(isComplex()){
-                    cblas_dcopy(copySize,src.getData()+1,getValueType(),getData()+1 + aBeginIndex * sliceStride*getValueType() ,getValueType());
+                    cblas_scopy(copySize,src.getData()+1,getValueType(),getData()+1 + aBeginIndex * sliceStride*getValueType() ,getValueType());
                 }
                 return true;
             }
@@ -961,7 +961,7 @@ namespace Aurora {
                 ::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];
+                    float value = getData()[(k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i)*complexstep+1];
                     if(value<0)
                     {
                         mark = "";
@@ -1006,7 +1006,7 @@ namespace Aurora {
         int rowOffset = aRowIdx == $ ? 0 : aRowIdx;
         int colOffset = aColIdx == $ ? 0 : aColIdx;
         int sliceOffset = aSliceIdx == $ ? 0 : aSliceIdx;
-        double *startPointer = getData() + (rowStride * rowOffset
+        float *startPointer = getData() + (rowStride * rowOffset
                                             + colStride * colOffset
                                             + sliceStride * sliceOffset) * getValueType();
         int size1 = allDimIndex.empty()?1:getDimSize(allDimIndex[0]);
@@ -1031,27 +1031,27 @@ namespace Aurora {
         }
     }
 
-    Matrix Matrix::operator>(double aScalar) const {
+    Matrix Matrix::operator>(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() > aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator>(double aScalar, const Matrix &matrix) {
+    Matrix operator>(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar>v.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1074,36 +1074,36 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)>matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() > v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
-    Matrix Matrix::operator>=(double aScalar) const {
+    Matrix Matrix::operator>=(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() >= aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator>=(double aScalar, const Matrix &matrix) {
+    Matrix operator>=(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar >= v.array() ).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1126,36 +1126,36 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)>=matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() >= v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
-    Matrix Matrix::operator<(double aScalar) const {
+    Matrix Matrix::operator<(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() < aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator<(double aScalar, const Matrix &matrix) {
+    Matrix operator<(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar < v.array() ).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1178,36 +1178,36 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)<matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() < v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
-    Matrix Matrix::operator<=(double aScalar) const {
+    Matrix Matrix::operator<=(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() <= aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator<=(double aScalar, const Matrix &matrix) {
+    Matrix operator<=(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar <= v.array() ).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1230,36 +1230,36 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)<=matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() <= v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
-    Matrix Matrix::operator==(double aScalar) const {
+    Matrix Matrix::operator==(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() == aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator==(double aScalar, const Matrix &matrix) {
+    Matrix operator==(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar == v.array() ).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1282,36 +1282,36 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)==matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() == v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
     }
 
-    Matrix Matrix::operator!=(double aScalar) const {
+    Matrix Matrix::operator!=(float aScalar) const {
         if (isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        double *ret = malloc(getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        float *ret = malloc(getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, getDataSize());
         result.setConstant(0.0);
         result = (v.array() != aScalar).select(1.0, result);
         return New(ret, getDimSize(0), getDimSize(1), getDimSize(2));
     }
 
-    Matrix operator!=(double aScalar, const Matrix &matrix) {
+    Matrix operator!=(float aScalar, const Matrix &matrix) {
         if (matrix.isComplex()) {
             std::cerr << "Complex cann't compare!" << std::endl;
             return Matrix();
         }
-        Eigen::Map<Eigen::VectorXd> v(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (aScalar != v.array() ).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1334,10 +1334,10 @@ namespace Aurora {
         if(matrix.isScalar()){
             return (*this)!=matrix.getData()[0];
         }
-        Eigen::Map<Eigen::VectorXd> v(getData(), getDataSize());
-        Eigen::Map<Eigen::VectorXd> v2(matrix.getData(), matrix.getDataSize());
-        double *ret = malloc(matrix.getDataSize());
-        Eigen::Map<Eigen::VectorXd> result(ret, matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> v(getData(), getDataSize());
+        Eigen::Map<Eigen::VectorXf> v2(matrix.getData(), matrix.getDataSize());
+        float *ret = malloc(matrix.getDataSize());
+        Eigen::Map<Eigen::VectorXf> result(ret, matrix.getDataSize());
         result.setConstant(0.0);
         result = (v.array() != v2.array()).select(1.0, result);
         return Matrix::New(ret, matrix.getDimSize(0), matrix.getDimSize(1), matrix.getDimSize(2));
@@ -1345,20 +1345,20 @@ namespace Aurora {
 
     Matrix operator-(Matrix &&aMatrix) {
         int size = aMatrix.getDataSize()*aMatrix.getValueType();
-        double zero = 0.0;
-        vdSubI(size,&zero,0,aMatrix.getData(),1,aMatrix.getData(),1);
+        float zero = 0.0;
+        vsSubI(size,&zero,0,aMatrix.getData(),1,aMatrix.getData(),1);
         return aMatrix;
     }
 
     Matrix operator-(const Matrix &aMatrix) {
-        double *newBuffer = malloc(aMatrix.getDataSize(), aMatrix.getValueType()==Complex);
-        double zero = 0.0;
+        float *newBuffer = malloc(aMatrix.getDataSize(), aMatrix.getValueType()==Complex);
+        float zero = 0.0;
         if (aMatrix.isComplex()){
-            vdSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData(),2,newBuffer,2);
-            vdSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData()+1,2,newBuffer+1,2);
+            vsSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData(),2,newBuffer,2);
+            vsSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData()+1,2,newBuffer+1,2);
         }
         else{
-            vdSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData(),1,newBuffer,1);
+            vsSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData(),1,newBuffer,1);
         }
         
         return Matrix::New(newBuffer,aMatrix);
@@ -1391,7 +1391,7 @@ namespace Aurora {
                 if(getDimSize(1) == aOther.getDimSize(0))
                 {
                     auto result = deepCopy();
-                    cblas_dgemv(CblasColMajor,CblasTrans,aOther.getDimSize(0),aOther.getDimSize(1),1.0,
+                    cblas_sgemv(CblasColMajor,CblasTrans,aOther.getDimSize(0),aOther.getDimSize(1),1.0,
                     aOther.getData(),aOther.getDimSize(0),getData(),1,0,result.getData(),1);
                     return result;
                 }
@@ -1411,7 +1411,7 @@ namespace Aurora {
             // right size
             if (getDimSize(1) == aOther.getDimSize(0)) {
                 auto result = aOther.deepCopy();
-                cblas_dgemv(CblasColMajor, CblasNoTrans, getDimSize(0),
+                cblas_sgemv(CblasColMajor, CblasNoTrans, getDimSize(0),
                             getDimSize(1), 1.0, getData(),
                             getDimSize(0), aOther.getData(), 1, 0,
                             result.getData(), 1);
@@ -1428,11 +1428,11 @@ namespace Aurora {
             // fit size
             if (getDimSize(1) == aOther.getDimSize(0))
             {
-                double * output = malloc(getDimSize(0)*aOther.getDimSize(1));
+                float * output = malloc(getDimSize(0)*aOther.getDimSize(1));
                 int M =  getDimSize(0);
                 int N =  aOther.getDimSize(1);
                 int K =  getDimSize(1);
-                cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.0, getData(), M, aOther.getData(),K, 0, output, M);
+                cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.0, getData(), M, aOther.getData(),K, 0, output, M);
                 return Matrix::New(output,M,N);
             }
             std::cerr << "Matrix mul operation fail, can't do Matrix("
@@ -1443,7 +1443,7 @@ namespace Aurora {
         }
     }
 
-    Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2):
+    Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, float* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2):
     mSliceMode(aSliceMode),mData(aData),
     mSize(aSize),mSize2(aSize2),
     mStride(aStride),mStride2(aStride2),
@@ -1482,28 +1482,28 @@ namespace Aurora {
         switch (mSliceMode) {
             case 2:{
                 if (mType== Normal) {
-//                    cblas_dcopy_batch_strided(mSize, slice.mData, slice.mStride, slice.mStride2, mData, mStride,
+//                    cblas_scopy_batch_strided(mSize, slice.mData, slice.mStride, slice.mStride2, mData, mStride,
 //                                              mStride2, mSize2);
                     for (int i = 0; i < mSize2; ++i) {
-                        cblas_dcopy(mSize,slice.mData + i*slice.mStride2,slice.mStride,mData+i*mStride2,mStride);
+                        cblas_scopy(mSize,slice.mData + i*slice.mStride2,slice.mStride,mData+i*mStride2,mStride);
                     }
                 }
                 else {
-//                    cblas_zcopy_batch_strided(mSize,(std::complex<double>*)slice.mData,slice.mStride,slice.mStride2,
-//                            (std::complex<double>*)mData,mStride,mStride2,mSize2);
+//                    cblas_zcopy_batch_strided(mSize,(std::complex<float>*)slice.mData,slice.mStride,slice.mStride2,
+//                            (std::complex<float>*)mData,mStride,mStride2,mSize2);
                     for (int i = 0; i < mSize2; ++i) {
-                        cblas_zcopy(mSize, (std::complex<double> *) (slice.mData + i * slice.mStride2), slice.mStride,
-                                    (std::complex<double> *) (mData + i * mStride2), mStride);
+                        cblas_zcopy(mSize, (std::complex<float> *) (slice.mData + i * slice.mStride2), slice.mStride,
+                                    (std::complex<float> *) (mData + i * mStride2), mStride);
                     }
                 }
                 break;
             }
             case 1:{
                 if (mType== Normal){
-                    cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride);
+                    cblas_scopy(mSize,slice.mData,slice.mStride,mData,mStride);
                 }
                 else {
-                    cblas_zcopy(mSize,  (std::complex<double> *)slice.mData, slice.mStride,  (std::complex<double> *)mData, mStride);
+                    cblas_zcopy(mSize,  (std::complex<float> *)slice.mData, slice.mStride,  (std::complex<float> *)mData, mStride);
                 }
                 break;
             }
@@ -1571,23 +1571,23 @@ namespace Aurora {
             //matrix mode
             case 2:{
                 if (mType== Normal) {
-                    cblas_dcopy_batch_strided(mSize, matrix.getData(), 1, matrix.getDimSize(0), mData, mStride,
+                    cblas_scopy_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);
+                    cblas_zcopy_batch_strided(mSize,(std::complex<float>*)matrix.getData(),1,matrix.getDimSize(0),
+                                              (std::complex<float>*)mData,mStride,mStride2,mSize2);
                 }
                 break;
             }
             //vector mode
             case 1:{
                 if (mType== Normal){
-                    cblas_dcopy(mSize,matrix.getData(),1,mData,mStride);
+                    cblas_scopy(mSize,matrix.getData(),1,mData,mStride);
                 }
                 else {
-                    cblas_zcopy(mSize, (std::complex<double> *) matrix.getData(),1,
-                                (std::complex<double> *) mData, mStride);
+                    cblas_zcopy(mSize, (std::complex<float> *) matrix.getData(),1,
+                                (std::complex<float> *) mData, mStride);
                 }
                 break;
             }
@@ -1601,7 +1601,7 @@ namespace Aurora {
         return *this;
     }
 
-    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(double value) {
+    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(float value) {
         if(!mData){
             std::cerr <<"Assign value fail!Des data pointer is null!";
             return *this;
@@ -1622,7 +1622,7 @@ namespace Aurora {
         return *this;
     }
 
-    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(std::complex<double> value) {
+    Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(std::complex<float> value) {
         if(!mData){
             std::cerr <<"Assign value fail!Des data pointer is null!";
             return *this;
@@ -1645,28 +1645,28 @@ namespace Aurora {
     }
 
     Matrix Matrix::MatrixSlice::toMatrix() const {
-        double * data = (double *) malloc(mSize*(mSize2>0?mSize2:1) ,mType == Complex);
+        float * data = (float *) malloc(mSize*(mSize2>0?mSize2:1) ,mType == Complex);
 
         switch (mSliceMode) {
             case 2:{
                 if (mType== Normal) {
-                    cblas_dcopy_batch_strided(mSize,  mData, mStride,
+                    cblas_scopy_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,
+                    cblas_zcopy_batch_strided(mSize, (std::complex<float> *) mData, mStride, mStride2,
+                                              (std::complex<float> *) data, 1, mSize,
                                               mSize2);
                 }
                 break;
             }
             case 1:{
                 if (mType== Normal){
-                    cblas_dcopy(mSize,mData,mStride,data,1);
+                    cblas_scopy(mSize,mData,mStride,data,1);
                 }
                 else {
-                    cblas_zcopy(mSize, (std::complex<double> *) mData, mStride,
-                                (std::complex<double> *) data, 1);
+                    cblas_zcopy(mSize, (std::complex<float> *) mData, mStride,
+                                (std::complex<float> *) data, 1);
                 }
                 break;
             }
diff --git a/src/Matrix.h b/src/Matrix.h
index 25fdb7e..d671af9 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -20,15 +20,15 @@ namespace Aurora {
          */
         class MatrixSlice{
         public:
-            MatrixSlice(int aSize,int aStride, double* aData,ValueType aType = Normal,int SliceMode = 1,int aSize2 = 0, int aStride2 = 0);
+            MatrixSlice(int aSize,int aStride, float* aData,ValueType aType = Normal,int SliceMode = 1,int aSize2 = 0, int aStride2 = 0);
             MatrixSlice& operator=(const MatrixSlice& slice);
             MatrixSlice& operator=(const Matrix& matrix);
-            MatrixSlice& operator=(double value);
-            MatrixSlice& operator=(std::complex<double> value);
+            MatrixSlice& operator=(float value);
+            MatrixSlice& operator=(std::complex<float> value);
             Matrix toMatrix() const;
         private:
             int mSliceMode = 0;//0 scalar, 1 vector, 2 Matrix
-            double* mData;
+            float* mData;
             int mSize=0;
             int mSize2=0;
             int mStride=1;
@@ -36,16 +36,16 @@ namespace Aurora {
             ValueType mType;
             friend class Matrix;
         };
-        explicit Matrix(std::shared_ptr<double> aData = std::shared_ptr<double>(),
+        explicit Matrix(std::shared_ptr<float> aData = std::shared_ptr<float>(),
                         std::vector<int> aInfo = std::vector<int>(),
                         ValueType aValueType = Normal);
 
         explicit Matrix(const Matrix::MatrixSlice& slice);
 
         /**
-         * Create from a Raw data(double array).
-         * Use Raw data which create like new double[size]() as a data source
-         * and the share_ptr's deleter will be std::default_delete<double[]>
+         * Create from a Raw data(float array).
+         * Use Raw data which create like new float[size]() as a data source
+         * and the share_ptr's deleter will be std::default_delete<float[]>
          * @param data
          * @param rows
          * @param cols
@@ -53,10 +53,10 @@ namespace Aurora {
          * @param type
          * @return
          */
-        static Matrix fromRawData(double *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
+        static Matrix fromRawData(float *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
 
         /**
-         * Create from a Raw data(double array) with copy the data to a new mkl memory.
+         * Create from a Raw data(float array) with copy the data to a new mkl memory.
          * @param data
          * @param rows
          * @param cols
@@ -64,7 +64,7 @@ namespace Aurora {
          * @param type
          * @return
          */
-        static Matrix copyFromRawData(double *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
+        static Matrix copyFromRawData(float *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
 
         /**
          * New a mkl calculate based Matrix
@@ -76,7 +76,7 @@ namespace Aurora {
          * @param type
          * @return
          */
-        static Matrix New(double *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
+        static Matrix New(float *data, int rows, int cols = 1, int slices = 1, ValueType type = Normal);
 
         /**
          * New a mkl calculate based Matrix
@@ -85,7 +85,7 @@ namespace Aurora {
          * @param shapeMatrix
          * @return
          */
-        static Matrix New(double *data, const Matrix &shapeMatrix);
+        static Matrix New(float *data, const Matrix &shapeMatrix);
 
         /**
          * New a mkl calculate based Matrix
@@ -105,19 +105,19 @@ namespace Aurora {
         MatrixSlice operator()(int r, int c = $, int aSliceIdx = $) const;
 
         // add
-        Matrix operator+(double aScalar) const;
-        friend Matrix operator+(double aScalar, const Matrix &matrix);
-        friend Matrix& operator+(double aScalar, Matrix &&matrix);
-        friend Matrix& operator+(Matrix &&matrix,double aScalar);
+        Matrix operator+(float aScalar) const;
+        friend Matrix operator+(float aScalar, const Matrix &matrix);
+        friend Matrix& operator+(float aScalar, Matrix &&matrix);
+        friend Matrix& operator+(Matrix &&matrix,float aScalar);
         Matrix operator+(const Matrix &matrix) const;
         Matrix operator+(Matrix &&matrix) const;
         friend Matrix operator+(Matrix &&aMatrix,Matrix &aOther);
 
         // sub
-        Matrix operator-(double aScalar) const;
-        friend Matrix operator-(double aScalar, const Matrix &matrix);
-        friend Matrix& operator-(double aScalar, Matrix &&matrix);
-        friend Matrix& operator-(Matrix &&matrix,double aScalar);
+        Matrix operator-(float aScalar) const;
+        friend Matrix operator-(float aScalar, const Matrix &matrix);
+        friend Matrix& operator-(float aScalar, Matrix &&matrix);
+        friend Matrix& operator-(Matrix &&matrix,float aScalar);
         Matrix operator-(const Matrix &matrix) const;
         Matrix operator-(Matrix &&matrix) const;
         friend Matrix operator-(Matrix &&aMatrix,Matrix &aOther);
@@ -127,19 +127,19 @@ namespace Aurora {
         friend Matrix operator-(const Matrix &aMatrix);
 
         // mul
-        Matrix operator*(double aScalar) const;
-        friend Matrix operator*(double aScalar, const Matrix &matrix);
-        friend Matrix& operator*(double aScalar, Matrix &&matrix);
-        friend Matrix& operator*(Matrix &&matrix,double aScalar);
+        Matrix operator*(float aScalar) const;
+        friend Matrix operator*(float aScalar, const Matrix &matrix);
+        friend Matrix& operator*(float aScalar, Matrix &&matrix);
+        friend Matrix& operator*(Matrix &&matrix,float aScalar);
         Matrix operator*(const Matrix &matrix) const;
         Matrix operator*(Matrix &&matrix) const;
         friend Matrix operator*(Matrix &&aMatrix,Matrix &aOther);
 
         // div
-        Matrix operator/(double aScalar) const;
-        friend Matrix operator/(double aScalar, const Matrix &matrix);
-        friend Matrix& operator/(double aScalar, Matrix &&matrix);
-        friend Matrix& operator/(Matrix &&matrix,double aScalar);
+        Matrix operator/(float aScalar) const;
+        friend Matrix operator/(float aScalar, const Matrix &matrix);
+        friend Matrix& operator/(float aScalar, Matrix &&matrix);
+        friend Matrix& operator/(Matrix &&matrix,float aScalar);
         Matrix operator/(const Matrix &matrix) const;
         Matrix operator/(Matrix &&matrix) const;
         friend Matrix operator/(Matrix &&aMatrix, Matrix &aOther);
@@ -149,33 +149,33 @@ namespace Aurora {
         friend Matrix operator^(Matrix &&matrix,int times);
 
         //compare
-        Matrix operator>(double aScalar) const;
-        friend Matrix operator>(double aScalar, const Matrix &matrix);
+        Matrix operator>(float aScalar) const;
+        friend Matrix operator>(float aScalar, const Matrix &matrix);
         Matrix operator>(const Matrix &matrix) const;
 
-        Matrix operator<(double aScalar) const;
-        friend Matrix operator<(double aScalar, const Matrix &matrix);
+        Matrix operator<(float aScalar) const;
+        friend Matrix operator<(float aScalar, const Matrix &matrix);
         Matrix operator<(const Matrix &matrix) const;
 
-        Matrix operator>=(double aScalar) const;
-        friend Matrix operator>=(double aScalar, const Matrix &matrix);
+        Matrix operator>=(float aScalar) const;
+        friend Matrix operator>=(float aScalar, const Matrix &matrix);
         Matrix operator>=(const Matrix &matrix) const;
 
-        Matrix operator<=(double aScalar) const;
-        friend Matrix operator<=(double aScalar, const Matrix &matrix);
+        Matrix operator<=(float aScalar) const;
+        friend Matrix operator<=(float aScalar, const Matrix &matrix);
         Matrix operator<=(const Matrix &matrix) const;
 
-        Matrix operator==(double aScalar) const;
-        friend Matrix operator==(double aScalar, const Matrix &matrix);
+        Matrix operator==(float aScalar) const;
+        friend Matrix operator==(float aScalar, const Matrix &matrix);
         Matrix operator==(const Matrix &matrix) const;
 
-        Matrix operator!=(double aScalar) const;
-        friend Matrix operator!=(double aScalar, const Matrix &matrix);
+        Matrix operator!=(float aScalar) const;
+        friend Matrix operator!=(float aScalar, const Matrix &matrix);
         Matrix operator!=(const Matrix &matrix) const;
 
         // sub
-        double& operator[](size_t index);
-        double operator[](size_t index) const;
+        float& operator[](size_t index);
+        float operator[](size_t index) const;
 
         /**
          * 切块操作
@@ -187,8 +187,8 @@ namespace Aurora {
          */
         Matrix block(int aDim,int aBeginIndx, int aEndIndex) const;
 
-        bool setBlockValue(int aDim,int aBeginIndx, int aEndIndex,double value);
-        bool setBlockComplexValue(int aDim,int aBeginIndx, int aEndIndex,std::complex<double> value);
+        bool setBlockValue(int aDim,int aBeginIndx, int aEndIndex,float value);
+        bool setBlockComplexValue(int aDim,int aBeginIndx, int aEndIndex,std::complex<float> value);
 
 
         bool setBlock(int aDim,int aBeginIndx, int aEndIndex,const Matrix& src);
@@ -218,7 +218,7 @@ namespace Aurora {
 
         bool isScalar() const;
 
-        double getScalar() const;
+        float getScalar() const;
 
         /**
          * Get is the Matrix's data is empty or size is zero.
@@ -238,7 +238,7 @@ namespace Aurora {
          */
         int getDims() const;
 
-        double *getData() const;
+        float *getData() const;
 
         int getDimSize(int index) const;
 
@@ -257,8 +257,8 @@ namespace Aurora {
 
         /**
          * Get Value type as normal and complex,
-         * complex use std::complex<double>,
-         * it's contains two double value.
+         * complex use std::complex<float>,
+         * it's contains two float value.
          * @return
          */
         ValueType getValueType() const {
@@ -289,7 +289,7 @@ namespace Aurora {
 
     private:
         ValueType mValueType = Normal;
-        std::shared_ptr<double> mData;
+        std::shared_ptr<float> mData;
         std::vector<int> mInfo;
         bool mMKL_Allocated = false;
     };