From 3b7c918cf1fe0e3f3934b32c4301bcef6436b1f4 Mon Sep 17 00:00:00 2001 From: Krad Date: Wed, 19 Apr 2023 15:54:52 +0800 Subject: [PATCH] Matrix class update. --- src/Matrix.cpp | 211 +++++++++++++++++++++++++++++++++++----- src/Matrix.h | 10 +- src/main.cxx | 83 +++++++++++++--- test/FunctionTester.cpp | 42 +++++++- 4 files changed, 296 insertions(+), 50 deletions(-) diff --git a/src/Matrix.cpp b/src/Matrix.cpp index 57dd81b..ffc788b 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -292,34 +292,67 @@ namespace Aurora { void Matrix::printf() { - ::printf("["); - for (int i = 0; i < mInfo[0]; ++i) { + int k_count = getDimSize(2)==0?1:getDimSize(2); + int j_count = getDimSize(1)==0?1:getDimSize(1); + for (int k = 0; k 0 ? getDimSize(0) : getDataSize()), 1, - getData() + getDimSize(0) * (aColIdx > 0 ? aColIdx : 0)); + Matrix::MatrixSlice Matrix::operator()(int aRowIdx, int aColIdx, int aSliceIdx) const { + std::vector dims = {aRowIdx, aColIdx, aSliceIdx}; + std::vector allDimIndex; + int mode = 0; + for (int j = 0; j < 3; ++j) { + if (dims[j]==$ && this->getDims()>j){ + ++mode; + allDimIndex.push_back(j); + } } - else if (aColIdx == ALL_DIM){ - return Matrix::MatrixSlice((int) (aRowIdx > 0 ? getDimSize(1) : getDataSize()), getDimSize(0), - getData() +(aRowIdx > 0 ? aRowIdx : 0)); + int rowStride = 1; + int colStride = getDimSize(0); + int sliceStride = getDimSize(0)*getDimSize(1); + int strides[3] = {rowStride, colStride, sliceStride}; + int rowOffset = aRowIdx == $ ? 0 : aRowIdx; + int colOffset = aColIdx == $ ? 0 : aColIdx; + int sliceOffset = aSliceIdx == $ ? 0 : aSliceIdx; + double *startPointer = getData() + (rowStride * rowOffset + + colStride * colOffset + + sliceStride * sliceOffset) * getValueType(); + int size1 = 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): - mData(aData), - mSize(aSize), - mStride(aStride), + Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2): + mSliceMode(aSliceMode),mData(aData), + mSize(aSize),mSize2(aSize2), + mStride(aStride),mStride2(aStride2), mType(aType) { @@ -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 (mType== Normal)cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride); - else cblas_zcopy(mSize,(std::complex*)slice.mData,slice.mStride,(std::complex*)mData,mStride); + if(!mData){ + std::cerr <<"Assign value fail!Des data pointer is null!"; + return *this; + } + if(!slice.mData){ + std::cerr <<"Assign value fail!Src data pointer is null!"; + return *this; + } + if (slice.mSliceMode!=mSliceMode) { + std::cerr <<"Assign value fail!Src slice(dims count:"<< slice.mSliceMode <<"), not match of des(dims count:"<*)slice.mData,slice.mStride,slice.mStride2, + (std::complex*)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 *) slice.mData, slice.mStride, + (std::complex *) 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 (mType== Normal)cblas_dcopy(mSize,matrix.getData(),1,mData,mStride); - else cblas_zcopy(mSize,(std::complex*)matrix.getDataSize(),1,(std::complex*)mData,mStride); + if(!mData){ + std::cerr <<"Assign value fail!Des data pointer is null!"; + return *this; + } + if(!matrix.getData()){ + std::cerr <<"Assign value fail!Src data pointer is null!"; + return *this; + } + if (matrix.getDims()!=mSliceMode) { + std::cerr <<"Assign value fail!Src matrix(dims count:"<< matrix.getDims() <<"), not match of des(dims count:"<*)matrix.getData(),1,matrix.getDimSize(0), + (std::complex*)mData,mStride,mStride2,mSize2); + } + break; + } + case 1:{ + if (mType== Normal){ + cblas_dcopy(mSize,matrix.getData(),1,mData,mStride); + } + else { + cblas_zcopy(mSize, (std::complex *) matrix.getData(),1, + (std::complex *) 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); - if (mType== Normal)cblas_dcopy(mSize,mData,mStride,data,1); - else cblas_zcopy(mSize,(std::complex*)mData,mStride,(std::complex*)data,1); - return Matrix::New(data,mSize,0,0,mType); + auto data = (double *) mkl_malloc(mSize*(mSize2>0?mSize2:1) * sizeof(double)*mType, 64); + auto matrix = Matrix::New(data,mSize,mSize2,0,mType); + switch (mSliceMode) { + 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 *) mData, mStride, mStride2, + (std::complex *) 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 *) mData, mStride, + (std::complex *) matrix.getData(), 1); + } + break; + } + case 0: + default:{ + matrix.getData()[0]= mData[0]; + if (mType != Normal) matrix.getData()[1] = mData[1]; + } + } + return matrix; } } diff --git a/src/Matrix.h b/src/Matrix.h index 885b04c..203df62 100644 --- a/src/Matrix.h +++ b/src/Matrix.h @@ -4,12 +4,13 @@ #include #include -#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; diff --git a/src/main.cxx b/src/main.cxx index c34cc70..b039de5 100644 --- a/src/main.cxx +++ b/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]; - for (int i = 0; i < 9; ++i) { - dataA[i]=(double)(i+1); - dataB[i]=(double)(i+2); + { + double * dataA =new double[8]; + double * dataB =new double[8]; + double * dataC =new double[4]; + 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(); - Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3); - printf("B:\r\n"); - B.printf(); - printf("A*B:\r\n"); - (A*B).printf(); - printf("A*B+A:\r\n"); - (A*B+A).printf(); - A = (A*B+A)+3.; - A.printf(); + { + 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(); + } + return 0; } diff --git a/test/FunctionTester.cpp b/test/FunctionTester.cpp index 952653e..d4525c7 100644 --- a/test/FunctionTester.cpp +++ b/test/FunctionTester.cpp @@ -1,10 +1,7 @@ -// -// Created by Krad on 2023/4/7. -// - #include #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];