Matrix class update.

src/Matrix.cpp

@@ -292,34 +292,67 @@ namespace Aurora {
 
 
     void Matrix::printf() {
-        ::printf("[");
+        int k_count = getDimSize(2)==0?1:getDimSize(2);
-        for (int i = 0; i < mInfo[0]; ++i) {
+        int j_count = getDimSize(1)==0?1:getDimSize(1);
+        for (int k = 0; k <k_count; ++k) {
+        ::printf("slice %d:\r\n[",k);
+        for (int i = 0; i < getDimSize(0); ++i) {
             ::printf("[");
-            for (int j = 0; j < mInfo[1]; ++j) {
+            for (int j = 0; j < j_count; ++j) {
-                ::printf("%f2, ",getData()[getDimSize(0)*j+i]);
+                ::printf("%f2, ",getData()[k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i]);
             }
             ::printf("]\r\n");
         }
         ::printf("]\r\n");
+        }
     }
 
-    Matrix::MatrixSlice Matrix::operator()(int aRowIdx, int aColIdx, int s) const {
+    Matrix::MatrixSlice Matrix::operator()(int aRowIdx, int aColIdx, int aSliceIdx) const {
-
+        std::vector<int> dims = {aRowIdx, aColIdx, aSliceIdx};
-        if (aRowIdx == ALL_DIM) {
+        std::vector<int> allDimIndex;
-            return Matrix::MatrixSlice((int) (aColIdx > 0 ? getDimSize(0) : getDataSize()), 1,
+        int mode = 0;
-                                       getData() + getDimSize(0) * (aColIdx > 0 ? aColIdx : 0));
+        for (int j = 0; j < 3; ++j) {
+            if (dims[j]==$ && this->getDims()>j){
+                ++mode;
+                allDimIndex.push_back(j);
+            }
         }
-        else if (aColIdx == ALL_DIM){
+        int rowStride = 1;
-            return Matrix::MatrixSlice((int) (aRowIdx > 0 ? getDimSize(1) : getDataSize()), getDimSize(0),
+        int colStride =  getDimSize(0);
-                                       getData() +(aRowIdx > 0 ? aRowIdx : 0));
+        int sliceStride =  getDimSize(0)*getDimSize(1);
+        int strides[3] = {rowStride, colStride, sliceStride};
+        int rowOffset = aRowIdx == $ ? 0 : aRowIdx;
+        int colOffset = aColIdx == $ ? 0 : aColIdx;
+        int sliceOffset = aSliceIdx == $ ? 0 : aSliceIdx;
+        double *startPointer = getData() + (rowStride * rowOffset
+                                            + colStride * colOffset
+                                            + sliceStride * sliceOffset) * getValueType();
+        int size1 = getDimSize(allDimIndex[0]);
+        int stride1 = strides[allDimIndex[0]];
+        switch (mode) {
+            //matrix mode
+            case 2:{
+                int size2 = getDimSize(allDimIndex[1]);
+                int stride2 = strides[allDimIndex[1]];
+                return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode, size2, stride2);
+            }
+            //vector mode & default
+            case 1:
+            {
+                return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode);
+            }
+            //scalar mode or ALL $
+            case 0:
+            default: {
+                return Matrix::MatrixSlice(1 , 1, startPointer,getValueType(), mode);
+            }
         }
-        return Matrix::MatrixSlice(1, 1, getData()+ getDimSize(0) * aColIdx +aRowIdx );
     }
 
-    Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType):
+    Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2):
-    mData(aData),
+    mSliceMode(aSliceMode),mData(aData),
-    mSize(aSize),
+    mSize(aSize),mSize2(aSize2),
-    mStride(aStride),
+    mStride(aStride),mStride2(aStride2),
     mType(aType)
     {
 
@@ -327,26 +360,150 @@ namespace Aurora {
 
     Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix::MatrixSlice &slice) {
         if (this==&slice) return *this;
-        if (mSize == slice.mSize && mType == slice.mType){
+        if(!mData){
-            if (mType== Normal)cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride);
+            std::cerr <<"Assign value fail!Des data pointer is null!";
-            else cblas_zcopy(mSize,(std::complex<double>*)slice.mData,slice.mStride,(std::complex<double>*)mData,mStride);
+            return *this;
+        }
+        if(!slice.mData){
+            std::cerr <<"Assign value fail!Src data pointer is null!";
+            return *this;
+        }
+        if (slice.mSliceMode!=mSliceMode) {
+            std::cerr <<"Assign value fail!Src slice(dims count:"<< slice.mSliceMode <<"), not match of des(dims count:"<<mSliceMode<<")!";
+            return *this;
+        }
+        if (slice.mSize!=mSize) {
+            std::cerr <<"Assign value fail!Src slice(dim 1 size:"<< slice.mSize <<"), not match of des(dim 1 size:"<<mSize<<")!";
+            return *this;
+        }
+        if (slice.mSize2!=mSize2) {
+            std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< slice.mSize2 <<"), not match of des(dim 2 size:"<<mSize2<<")!";
+            return *this;
+        }
+        if (slice.mType!=mType) {
+            std::cerr <<"Assign value fail!Src slice(value type:"<< slice.mType <<"), not match of des(value type:"<<mType<<")!";
+            return *this;
+        }
+        switch (mSliceMode) {
+            case 2:{
+                if (mType== Normal) {
+                    cblas_dcopy_batch_strided(mSize, slice.mData, slice.mStride, slice.mStride2, mData, mStride,
+                                              mStride2, mSize2);
+                }
+                else {
+                    cblas_zcopy_batch_strided(mSize,(std::complex<double>*)slice.mData,slice.mStride,slice.mStride2,
+                            (std::complex<double>*)mData,mStride,mStride2,mSize2);
+                }
+                break;
+            }
+            case 1:{
+                if (mType== Normal){
+                    cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride);
+                }
+                else {
+                    cblas_zcopy(mSize, (std::complex<double> *) slice.mData, slice.mStride,
+                                (std::complex<double> *) mData, mStride);
+                }
+                break;
+            }
+            case 0:
+            default:{
+                mData[0] = slice.mData[0];
+                if (mType != Normal)mData[1] = slice.mData[1];
+            }
         }
         return *this;
     }
 
     Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix &matrix) {
-        if (mSize == matrix.getDataSize() && mType == matrix.getValueType()){
+        if(!mData){
-            if (mType== Normal)cblas_dcopy(mSize,matrix.getData(),1,mData,mStride);
+            std::cerr <<"Assign value fail!Des data pointer is null!";
-            else cblas_zcopy(mSize,(std::complex<double>*)matrix.getDataSize(),1,(std::complex<double>*)mData,mStride);
+            return *this;
+        }
+        if(!matrix.getData()){
+            std::cerr <<"Assign value fail!Src data pointer is null!";
+            return *this;
+        }
+        if (matrix.getDims()!=mSliceMode) {
+            std::cerr <<"Assign value fail!Src matrix(dims count:"<< matrix.getDims() <<"), not match of des(dims count:"<<mSliceMode<<")!";
+            return *this;
+        }
+        if (matrix.getDimSize(0)!=mSize) {
+            std::cerr <<"Assign value fail!Src matrix(dim 1 size:"<< matrix.getDimSize(0)<<"), not match of des(dim 1 size:"<<mSize<<")!";
+            return *this;
+        }
+        if (matrix.getDimSize(1)!=mSize2) {
+            std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< matrix.getDimSize(1) <<"), not match of des(dim 2 size:"<<mSize2<<")!";
+            return *this;
+        }
+        if (matrix.getValueType()!=mType) {
+            std::cerr <<"Assign value fail!Src slice(value type:"<< matrix.getValueType() <<"), not match of des(value type:"<<mType<<")!";
+            return *this;
+        }
+        switch (mSliceMode) {
+            case 2:{
+                if (mType== Normal) {
+                    cblas_dcopy_batch_strided(mSize, matrix.getData(), 1, matrix.getDimSize(0), mData, mStride,
+                                              mStride2, mSize2);
+                }
+                else {
+                    cblas_zcopy_batch_strided(mSize,(std::complex<double>*)matrix.getData(),1,matrix.getDimSize(0),
+                                              (std::complex<double>*)mData,mStride,mStride2,mSize2);
+                }
+                break;
+            }
+            case 1:{
+                if (mType== Normal){
+                    cblas_dcopy(mSize,matrix.getData(),1,mData,mStride);
+                }
+                else {
+                    cblas_zcopy(mSize, (std::complex<double> *) matrix.getData(),1,
+                                (std::complex<double> *) mData, mStride);
+                }
+                break;
+            }
+            case 0:
+            default:{
+                mData[0] = matrix.getData()[0];
+                if (mType != Normal)mData[1] = matrix.getData()[1];
+            }
         }
         return *this;
     }
 
     Matrix Matrix::MatrixSlice::toMatrix() {
-        auto data = malloc(mSize ,mType == Complex);
+        auto data = (double *) mkl_malloc(mSize*(mSize2>0?mSize2:1) * sizeof(double)*mType, 64);
-        if (mType== Normal)cblas_dcopy(mSize,mData,mStride,data,1);
+        auto matrix = Matrix::New(data,mSize,mSize2,0,mType);
-        else cblas_zcopy(mSize,(std::complex<double>*)mData,mStride,(std::complex<double>*)data,1);
+        switch (mSliceMode) {
-        return Matrix::New(data,mSize,0,0,mType);
+            case 2:{
+                if (mType== Normal) {
+                    cblas_dcopy_batch_strided(mSize,  mData, mStride,
+                                              mStride2,matrix.getData(), 1, matrix.getDimSize(0), mSize2);
+                }
+                else {
+                    cblas_zcopy_batch_strided(mSize, (std::complex<double> *) mData, mStride, mStride2,
+                                              (std::complex<double> *) matrix.getData(), 1, matrix.getDimSize(0),
+                                              mSize2);
+                }
+                break;
+            }
+            case 1:{
+                if (mType== Normal){
+                    cblas_dcopy(mSize,mData,mStride,matrix.getData(),1);
+                }
+                else {
+                    cblas_zcopy(mSize, (std::complex<double> *) mData, mStride,
+                                (std::complex<double> *) matrix.getData(), 1);
+                }
+                break;
+            }
+            case 0:
+            default:{
+                matrix.getData()[0]= mData[0];
+                if (mType != Normal) matrix.getData()[1] = mData[1];
+            }
+        }
+        return matrix;
     }
 }
 

src/Matrix.h

@@ -4,12 +4,13 @@
 #include <memory>
 #include <vector>
 
-#define ALL_DIM -1
+
 namespace Aurora {
     enum ValueType{
         Normal=1,
         Complex
     };
+    const int $ = -1;
     
     class Matrix {
     public:
@@ -28,14 +29,17 @@ namespace Aurora {
          */
         class MatrixSlice{
         public:
-            MatrixSlice(int aSize,int aStride, double* aData,ValueType aType = Normal);
+            MatrixSlice(int aSize,int aStride, double* aData,ValueType aType = Normal,int SliceMode = 1,int aSize2 = 0, int aStride2 = 0);
             MatrixSlice& operator=(const MatrixSlice& slice);
             MatrixSlice& operator=(const Matrix& matrix);
             Matrix toMatrix();
         private:
+            int mSliceMode = 0;//0 scalar, 1 vector, 2 Matrix
             double* mData;
             int mSize;
+            int mSize2;
             int mStride;
+            int mStride2;
             ValueType mType;
         };
 
@@ -50,7 +54,7 @@ namespace Aurora {
         Matrix deepCopy() const;
 
         //切片，暂时不支持三维
-        MatrixSlice operator()(int r, int c = ALL_DIM, int s = ALL_DIM) const;
+        MatrixSlice operator()(int r, int c = $, int aSliceIdx = $) const;
 
         // add
         Matrix operator+(double aScalar) const;

src/main.cxx

@@ -8,28 +8,79 @@
 
 
 #include "Matrix.h"
+#include "Function1D.h"
 
 #include "spdlog/spdlog.h"
 #include "spdlog/sinks/stdout_color_sinks.h"
 #include "spdlog/sinks/basic_file_sink.h"
 
 int main() {
-    double * dataA =new double[9];
+    {
-    double * dataB =new double[9];
+        double * dataA =new double[8];
-    for (int i = 0; i < 9; ++i) {
+        double * dataB =new double[8];
-        dataA[i]=(double)(i+1);
+        double * dataC =new double[4];
-        dataB[i]=(double)(i+2);
+        for (int i = 0; i < 8; ++i) {
+            dataA[i]=(double)(i-3);
+            dataB[i]=(double)(i+2);
+            dataC[i/2]=(double)(i+9);
+        }
+        Aurora::Matrix A = Aurora::Matrix::New(dataA,2,2,2);
+        printf("A:\r\n");
+        A.printf();
+        Aurora::Matrix B = Aurora::Matrix::New(dataB,2,2,2);
+        printf("B:\r\n");
+        B.printf();
+        printf("A slice 1 = B slice 0\r\n");
+        A(Aurora::$,Aurora::$,1) = B(Aurora::$,Aurora::$,0);
+        printf("B:\r\n");
+        B.printf();
+        printf("New A:\r\n");
+        A.printf();
+        printf("A Row slice 1 = B Row slice 0\r\n");
+        A(1,Aurora::$,Aurora::$) = B(0,Aurora::$,Aurora::$);
+        printf("B:\r\n");
+        B.printf();
+        printf("New A:\r\n");
+        A.printf();
+
+        printf("A Col slice 1 = B Col slice 0\r\n");
+        A(Aurora::$,1,Aurora::$) = B(Aurora::$,0,Aurora::$);
+        printf("B:\r\n");
+        B.printf();
+        printf("New A:\r\n");
+        A.printf();
+
+            printf("New A col slice 1 toMatrix:\r\n");
+            auto Ds = A(Aurora::$,1,Aurora::$);
+            auto D = Ds.toMatrix();
+            D.printf();
+
+
+        Aurora::Matrix C = Aurora::Matrix::New(dataC,2,2);
+        printf("C:\r\n");
+        C.printf();
+        A(Aurora::$,Aurora::$,0) = C;
+        printf("New A:\r\n");
+        A.printf();
+        return 0;
+
     }
-    Aurora::Matrix A = Aurora::Matrix::New(dataA,3,3);
+    {
-    A.printf();
+        double * dataA =new double[9];
-    Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3);
+        double * dataB =new double[9];
-    printf("B:\r\n");
+        for (int i = 0; i < 9; ++i) {
-    B.printf();
+            dataA[i]=(double)(i-3);
-    printf("A*B:\r\n");
+            dataB[i]=(double)(i+2);
-    (A*B).printf();
+        }
-    printf("A*B+A:\r\n");
+        Aurora::Matrix A = Aurora::Matrix::New(dataA,3,3);
-    (A*B+A).printf();
+        printf("A:\r\n");
-    A = (A*B+A)+3.;
+        A.printf();
-    A.printf();
+        Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3);
+        printf("B:\r\n");
+        B.printf();
+        printf("sign(A):\r\n");
+        sign(A*B).printf();
+    }
+
     return 0;
 }

test/FunctionTester.cpp

@@ -1,10 +1,7 @@
-//
-// Created by Krad on 2023/4/7.
-//
-
 #include <gtest/gtest.h>
 #include "Matrix.h"
 #include "Function.h"
+#include "Function1D.h"
 
 
 class FunctionTester:public ::testing::Test{
@@ -24,6 +21,43 @@ double fourDecimalRound(double src){
     return round(src*10000.0)/10000.0;
 }
 
+TEST_F(FunctionTester, matrixSlice) {
+    double * dataA =new double[8];
+    double * dataB =new double[8];
+    for (int i = 0; i < 8; ++i) {
+        dataA[i]=(double)(i-3);
+        dataB[i]=(double)(i+2);
+    }
+    Aurora::Matrix A = Aurora::Matrix::New(dataA,2,2,2);
+    printf("A:\r\n");
+    A.printf();
+    Aurora::Matrix B = Aurora::Matrix::New(dataB,2,2,2);
+    printf("B:\r\n");
+    B.printf();
+    A(Aurora::$,Aurora::$,1) = B(Aurora::$,Aurora::$,0);
+    printf("New A:\r\n");
+    A.printf();
+    printf("New B:\r\n");
+    B.printf();
+}
+
+TEST_F(FunctionTester, sign) {
+    double * dataA =new double[9];
+    double * dataB =new double[9];
+    for (int i = 0; i < 9; ++i) {
+        dataA[i]=(double)(i-3);
+        dataB[i]=(double)(i+2);
+    }
+    Aurora::Matrix A = Aurora::Matrix::New(dataA,3,3);
+    printf("A:\r\n");
+    A.printf();
+    Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3);
+    printf("B:\r\n");
+    B.printf();
+    printf("sign(A):\r\n");
+    sign(A*B).printf();
+}
+
 TEST_F(FunctionTester, matrix) {
     double * dataA =new double[9];
     double * dataB =new double[9];