Add Cuda matrix size function & unittest

src/Function3D.cpp

@@ -1,7 +1,12 @@
 #include <iostream>
+#include "CudaMatrix.h"
 #include "Function3D.h"
 #include "Function2D.h"
 #include "Function.h"
+#include "CudaMatrixPrivate.cuh"
+#include <cuda_runtime.h>
+
+
 
 //必须在Eigen之前
 #include "AuroraDefs.h"
@@ -61,14 +66,34 @@ Matrix Aurora::ones(int aRow, int aColumn, int aSlice) {
     size_t arraySize  = rowSize * colSize* sliceSize;
     float* data = malloc(arraySize);
     float one = 1.0;
-    cblas_scopy(arraySize,&one,0,data,1);
+    cblas_scopy(arraySize,&one,0,data,1.0f);
     return Matrix::New(data,rowSize,colSize,aSlice);
 }
 
+CudaMatrix Aurora::onesCuda(int aRow, int aColumn, int aSlice){
+    if (aRow == 0 || aColumn == 0)
+    {
+        std::cerr<<"ones function can create matrix with dim unit cont =0";
+        return CudaMatrix();
+    }
+    int rowSize =  aRow;
+    int colSize =  aColumn;
+    int sliceSize = aSlice == 0 ? 1 : aSlice;
+    size_t arraySize  = rowSize * colSize* sliceSize;
+    float* data = nullptr;
+    cudaMalloc((void**)&data,arraySize*sizeof(float));
+    ::thrustFill(data,data+arraySize,1.0f);
+    return CudaMatrix::fromRawData(data,rowSize,colSize,sliceSize);
+}
+
 Matrix Aurora::ones(int aSquareRow) {
     return Aurora::ones(aSquareRow, aSquareRow);
 }
 
+CudaMatrix Aurora::onesCuda(int aSquareRow) {
+    return Aurora::onesCuda(aSquareRow, aSquareRow);
+}
+
 Matrix Aurora::zeros(int aRow, int aColumn, int aSlice) {
     if (aRow == 0 || aColumn == 0)
     {
@@ -81,18 +106,39 @@ Matrix Aurora::zeros(int aRow, int aColumn, int aSlice) {
     size_t arraySize  = rowSize * colSize* sliceSize;
     float* data = malloc(arraySize);
     float zero = 0.0;
-    cblas_scopy(arraySize,&zero,0,data,1);
+    cblas_scopy(arraySize,&zero,0,data,0.0f);
     return Matrix::New(data,rowSize,colSize,sliceSize);
 }
 
+CudaMatrix Aurora::zerosCuda(int aRow, int aColumn, int aSlice) {
+    if (aRow == 0 || aColumn == 0)
+    {
+        std::cerr<<"zeros function can create matrix with dim unit cont =0";
+        return CudaMatrix();
+    }
+    int rowSize =  aRow;
+    int colSize =  aColumn;
+    int sliceSize = aSlice == 0 ? 1 : aSlice;
+    size_t arraySize  = rowSize * colSize* sliceSize;
+    float* data = nullptr;
+    cudaMalloc((void**)&data,arraySize*sizeof(float));
+    ::thrustFill(data,data+arraySize,0.0f);
+    return CudaMatrix::fromRawData(data,rowSize,colSize,sliceSize);
+}
+
+
 Matrix Aurora::zeros(int aSquareRow) {
     return Aurora::zeros(aSquareRow, aSquareRow);
 }
 
+CudaMatrix Aurora::zerosCuda(int aSquareRow) {
+    return Aurora::zerosCuda(aSquareRow, aSquareRow);
+}
+
 Matrix Aurora::size(const Matrix &aMatrix)
 {
     if (aMatrix.isScalar()){
-        float * output = Aurora::malloc(1);
+        float * output =nullptr;
         output[0]=1;
         return Matrix::New(output,1,1,1);
     }
@@ -119,11 +165,50 @@ Matrix Aurora::size(const Matrix &aMatrix)
     }
 }
 
+CudaMatrix Aurora::size(const CudaMatrix &aMatrix){
+    float * output=nullptr;
+    if (aMatrix.isScalar()){
+        cudaMalloc((void**)&output,sizeof(float));
+        auto outMatrix = CudaMatrix::fromRawData(output,1,1,1);
+        outMatrix.setValue(0, 1);
+        return outMatrix;
+    }
+    else if (aMatrix.isVector()){
+        cudaMalloc((void**)&output,sizeof(float)*2);
+        auto outMatrix = CudaMatrix::fromRawData(output,2,1,1);
+        outMatrix.setValue(0, aMatrix.getDimSize(0));
+        outMatrix.setValue(1, aMatrix.getDimSize(1));
+        return outMatrix;
+    }
+    //3D
+    else if (aMatrix.getDimSize(2)>1){
+        cudaMalloc((void**)&output,sizeof(float)*3);
+        auto outMatrix = CudaMatrix::fromRawData(output,3,1,1);
+        outMatrix.setValue(0,aMatrix.getDimSize(0));
+        outMatrix.setValue(1,aMatrix.getDimSize(1));
+        outMatrix.setValue(2,aMatrix.getDimSize(2));
+        return outMatrix;
+    }
+    //2D matrix
+    else{
+        cudaMalloc((void**)&output,sizeof(float)*2);
+        auto outMatrix = CudaMatrix::fromRawData(output,2,1,1);
+        outMatrix.setValue(0,aMatrix.getDimSize(0));
+        outMatrix.setValue(1,aMatrix.getDimSize(1));
+        return outMatrix;
+    }
+}
+
 int Aurora::size(const Matrix &aMatrix,int dims)
 {
     return aMatrix.getDimSize(dims-1);
 }
 
+int Aurora::size(const CudaMatrix &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, float aExtrapval)
 {
     std::vector<Matrix> zTemps;

src/Function3D.h

@@ -4,6 +4,7 @@
 
 #include "Matrix.h"
 #include "Function1D.h"
+#include "CudaMatrix.h"
 
 namespace Aurora {
 
@@ -16,6 +17,8 @@ namespace Aurora {
      */
     Matrix ones(int aRow, int aColumn, int aSlice = 0);
 
+    CudaMatrix onesCuda(int aRow, int aColumn, int aSlice = 0);
+
     /**
      * 创建全部为1的方阵
      * @param aSquareRow
@@ -23,6 +26,8 @@ namespace Aurora {
      */
     Matrix ones(int aSquareRow);
 
+    CudaMatrix onesCuda(int aSquareRow);
+
     /**
      * 创建全部为0的数组，矩阵
      * @param aRow 行数，必须大于0
@@ -32,18 +37,25 @@ namespace Aurora {
      */
     Matrix zeros(int aRow, int aColumn, int aSlice = 0);
 
+    CudaMatrix zerosCuda(int aRow, int aColumn, int aSlice = 0);
+
     /**
      * 创建全部为0的方阵
      * @param aSquareRow
      * @return 全部为0的方阵
      */
     Matrix zeros(int aSquareRow);
+    CudaMatrix zerosCuda(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, 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);
+    CudaMatrix size(const CudaMatrix &aMatrix);
+
     int size(const Matrix &aMatrix,int dims);
+    int size(const CudaMatrix &aMatrix,int dims);
+
 };
 
 

test/Function3D_Test.cpp

@@ -4,6 +4,8 @@
 #include "TestUtility.h"
 
 #include "Matrix.h"
+#include "CudaMatrix.h"
+
 #include "Function.h"
 #include "Function1D.h"
 #include "Function2D.h"
@@ -93,50 +95,98 @@ TEST_F(Function3D_Test, interpn) {
     EXPECT_FLOAT_AE(result.getData()[2],94.7908);
 }
 
-TEST_F(Function3D_Test, zerosAndones){
+TEST_F(Function3D_Test, zerosAndonesAndsizesMKL) {
-
+    // mklversion
-    Aurora::Matrix zerosM = Aurora::zeros( 3, 4,5);
+    Aurora::Matrix zerosM = Aurora::zeros(3, 4, 5);
-    EXPECT_EQ(60,zerosM.getDataSize());
+    EXPECT_EQ(60, zerosM.getDataSize());
+    auto size = Aurora::size(zerosM);
+    EXPECT_EQ(3,size[0]);
+    EXPECT_EQ(4,size[1]);
+    EXPECT_EQ(5,size[2]);
+    EXPECT_EQ(3,Aurora:: size(zerosM,1));
+    EXPECT_EQ(4,Aurora:: size(zerosM,2));
+    EXPECT_EQ(5,Aurora:: size(zerosM,3));
     for (int i = 0; i < zerosM.getDataSize(); ++i) {
-        EXPECT_EQ(0,zerosM.getData()[i])<<" error at index："<<i;
+        EXPECT_EQ(0, zerosM.getData()[i]) << " error at index:" << i;
     }
-    zerosM = Aurora::zeros( 3, 4);
+    zerosM = Aurora::zeros(3, 4);
-    EXPECT_EQ(12,zerosM.getDataSize());
+    EXPECT_EQ(12, zerosM.getDataSize());
     for (int i = 0; i < zerosM.getDataSize(); ++i) {
-        EXPECT_EQ(0,zerosM.getData()[i])<<" error at index："<<i;
+      EXPECT_EQ(0, zerosM.getData()[i]) << " error at index:" << i;
     }
-    zerosM = Aurora::zeros( 3, 1);
+    zerosM = Aurora::zeros(3, 1);
-    EXPECT_EQ(3,zerosM.getDataSize());
+    EXPECT_EQ(3, zerosM.getDataSize());
     for (int i = 0; i < zerosM.getDataSize(); ++i) {
-        EXPECT_EQ(0,zerosM.getData()[i])<<" error at index："<<i;
+      EXPECT_EQ(0, zerosM.getData()[i]) << " error at index:" << i;
     }
-    zerosM = Aurora::zeros( 20);
+    zerosM = Aurora::zeros(20);
-    EXPECT_EQ(400,zerosM.getDataSize());
+    EXPECT_EQ(400, zerosM.getDataSize());
     for (int i = 0; i < zerosM.getDataSize(); ++i) {
-        EXPECT_EQ(0,zerosM.getData()[i])<<" error at index："<<i;
+      EXPECT_EQ(0, zerosM.getData()[i]) << " error at index:" << i;
     }
-    Aurora::Matrix onesM = Aurora::ones( 9, 9,9);
+    Aurora::Matrix onesM = Aurora::ones(9, 9, 9);
-    EXPECT_EQ(729,onesM.getDataSize());
+    EXPECT_EQ(729, onesM.getDataSize());
     for (int i = 0; i < onesM.getDataSize(); ++i) {
-        EXPECT_FLOAT_EQ(1.0,onesM.getData()[i])<<" error at index："<<i;
+      EXPECT_FLOAT_EQ(1.0, onesM.getData()[i]) << " error at index:" << i;
     }
 
-    onesM = Aurora::ones( 9, 9);
+    onesM = Aurora::ones(9, 9);
-    EXPECT_EQ(81,onesM.getDataSize());
+    EXPECT_EQ(81, onesM.getDataSize());
     for (int i = 0; i < onesM.getDataSize(); ++i) {
-        EXPECT_FLOAT_EQ(1.0,onesM.getData()[i])<<" error at index："<<i;
+      EXPECT_FLOAT_EQ(1.0, onesM.getData()[i]) << " error at index:" << i;
     }
 
-    onesM = Aurora::ones( 9, 1);
+    onesM = Aurora::ones(9, 1);
-    EXPECT_EQ(9,onesM.getDataSize());
+    EXPECT_EQ(9, onesM.getDataSize());
     for (int i = 0; i < onesM.getDataSize(); ++i) {
-        EXPECT_FLOAT_EQ(1.0,onesM.getData()[i])<<" error at index："<<i;
+      EXPECT_FLOAT_EQ(1.0, onesM.getData()[i]) << " error at index:" << i;
     }
 
-    onesM = Aurora::ones( 9);
+    onesM = Aurora::ones(9);
-    EXPECT_EQ(81,onesM.getDataSize());
+    EXPECT_EQ(81, onesM.getDataSize());
     for (int i = 0; i < onesM.getDataSize(); ++i) {
-        EXPECT_FLOAT_EQ(1.0,onesM.getData()[i])<<" error at index："<<i;
+      EXPECT_FLOAT_EQ(1.0, onesM.getData()[i]) << " error at index:" << i;
+    }
+}
+
+TEST_F(Function3D_Test, zerosAndonesAndsizeCuda) {
+    Aurora::CudaMatrix zerosM = Aurora::zerosCuda(3, 4, 5);
+    auto size = Aurora::size(zerosM);
+    EXPECT_EQ(3,size.getValue(0));
+    EXPECT_EQ(4,size.getValue(1));
+    EXPECT_EQ(5,size.getValue(2));
+    EXPECT_EQ(3,Aurora:: size(zerosM,1));
+    EXPECT_EQ(4,Aurora:: size(zerosM,2));
+    EXPECT_EQ(5,Aurora:: size(zerosM,3));
+    EXPECT_EQ(60, zerosM.getDataSize());
+    for (int i = 0; i < zerosM.getDataSize(); ++i) {
+        EXPECT_EQ(0, zerosM.getValue(i)) << " error at index:" << i;
+    }
+    zerosM = Aurora::zerosCuda(3, 4);
+    EXPECT_EQ(12, zerosM.getDataSize());
+    for (int i = 0; i < zerosM.getDataSize(); ++i) {
+      EXPECT_EQ(0, zerosM.getValue(i)) << " error at index:" << i;
+    }
+    zerosM = Aurora::zerosCuda(20);
+    EXPECT_EQ(400, zerosM.getDataSize());
+    for (int i = 0; i < zerosM.getDataSize(); ++i) {
+      EXPECT_EQ(0, zerosM.getValue(i)) << " error at index:" << i;
+    }
+    auto onesM = Aurora::onesCuda(9, 9, 9);
+    EXPECT_EQ(729, onesM.getDataSize());
+    for (int i = 0; i < onesM.getDataSize(); ++i) {
+      ASSERT_FLOAT_EQ(1.0, onesM.getValue(i)) << " error at index:" << i;
     }
 
+    onesM = Aurora::onesCuda(9, 3);
+    EXPECT_EQ(27, onesM.getDataSize());
+    for (int i = 0; i < onesM.getDataSize(); ++i) {
+      EXPECT_FLOAT_EQ(1.0, onesM.getValue(i)) << " error at index:" << i;
+    }
+
+    onesM = Aurora::onesCuda(9);
+    EXPECT_EQ(81, onesM.getDataSize());
+    for (int i = 0; i < onesM.getDataSize(); ++i) {
+      EXPECT_FLOAT_EQ(1.0, onesM.getValue(i)) << " error at index:" << i;
+    }
 }