1048 lines
41 KiB
C++
1048 lines
41 KiB
C++
#include "Matrix.h"
|
|
|
|
#include <string>
|
|
#include <cstring>
|
|
#include <iostream>
|
|
|
|
#include "AuroraDefs.h"
|
|
|
|
#include <Eigen/Core>
|
|
#include <Eigen/Dense>
|
|
|
|
#include "Function.h"
|
|
|
|
namespace Aurora{
|
|
typedef void(*CalcFuncD)(const MKL_INT n, const double a[], const MKL_INT inca, const double b[],
|
|
const MKL_INT incb, double r[], const MKL_INT incr);
|
|
|
|
typedef void(*CalcFuncZ)(const MKL_INT n, const MKL_Complex16 a[], const MKL_INT inca, const MKL_Complex16 b[],
|
|
const MKL_INT incb, MKL_Complex16 r[], const MKL_INT incr);
|
|
|
|
inline Matrix operatorMxA(CalcFuncD aFunc, double aScalar, const Matrix &aMatrix) {
|
|
double *output = malloc(aMatrix.getDataSize(), aMatrix.getValueType() == Complex);
|
|
//复数时，+和-需要特别操作，只影响实部
|
|
if (aMatrix.getValueType() == Complex && (aFunc == &vdAddI || aFunc == &vdSubI)) {
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData() , 2, &aScalar, 0, output ,
|
|
2);
|
|
double zero = 0.0;
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData()+1 , 2, &zero, 0, output+1 ,
|
|
2);
|
|
} else{
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0, output, 1);
|
|
}
|
|
return Matrix::New(output, aMatrix.getDimSize(0), aMatrix.getDimSize(1), aMatrix.getDimSize(2),
|
|
aMatrix.getValueType());
|
|
}
|
|
|
|
inline Matrix &operatorMxA_RR(CalcFuncD aFunc, double aScalar, Aurora::Matrix &&aMatrix) {
|
|
if (aMatrix.getValueType() == Complex) {
|
|
//针对实部操作
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData(), 2, &aScalar, 0,
|
|
aMatrix.getData(),
|
|
2);
|
|
//乘法除法需特别操作，影响虚部
|
|
if (aFunc == &vdDivI || aFunc == &vdMulI){
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData() + 1, 2, &aScalar, 0,
|
|
aMatrix.getData() + 1, 2);
|
|
}
|
|
} else {
|
|
aFunc(aMatrix.getDataSize(), aMatrix.getData(), 1, &aScalar, 0,
|
|
aMatrix.getData(),
|
|
1);
|
|
}
|
|
return aMatrix;
|
|
}
|
|
|
|
|
|
inline void V_MxM_CN_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* xC,double* yN,double *output, int DimsStride) {
|
|
aFuncD(size, xC, DimsStride * 2, yN, 1, output, 2);
|
|
if (aFuncD == &vdDivI || aFuncD == &vdMulI){
|
|
aFuncD(size, xC + 1, DimsStride * 2, yN, 1, output + 1, 2);
|
|
}
|
|
else{
|
|
double zero = 0.0;
|
|
aFuncD(size, xC+1 , DimsStride*2, &zero, 0, output + 1, 2);
|
|
}
|
|
}
|
|
|
|
inline double* _MxM_CN_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* xC,double* yN, int dimsStride)
|
|
{
|
|
double *output = malloc(size, true);
|
|
V_MxM_CN_Calc(aFuncD, size, xC, yN, output, dimsStride);
|
|
return output;
|
|
}
|
|
|
|
inline void V_MxM_NC_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* xN,double* yC,double *output, int DimsStride) {
|
|
aFuncD(size, xN, DimsStride, yC, 2, output, 2);
|
|
if (aFuncD == &vdDivI || aFuncD == &vdMulI){
|
|
aFuncD(size, xN , DimsStride, yC+ 1, 2, output + 1, 2);
|
|
}
|
|
else{
|
|
double zero = 0.0;
|
|
aFuncD(size, yC+ 1, 2, &zero, 0, output + 1, 2);
|
|
}
|
|
}
|
|
|
|
inline double* _MxM_NC_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* xN,double* yC, int dimsStride)
|
|
{
|
|
double *output = malloc(size, true);
|
|
V_MxM_NC_Calc(aFuncD, size, xN, yC, output, dimsStride);
|
|
return output;
|
|
}
|
|
|
|
inline void V_MxM_NN_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* x,double* y,double* output, int DimsStride) {
|
|
aFuncD(size, x, DimsStride, y, 1, output,1);
|
|
}
|
|
|
|
inline double* _MxM_NN_Calc(
|
|
CalcFuncD aFuncD,
|
|
const int size, double* x,double* y, int DimsStride) {
|
|
double *output = malloc(size);
|
|
V_MxM_NN_Calc(aFuncD, size, x, y, output, DimsStride);
|
|
return output;
|
|
}
|
|
|
|
inline void V_MxM_CC_Calc(
|
|
CalcFuncZ aFuncZ, const int size, double* x,double* y,double* output,
|
|
int DimsStride) {
|
|
aFuncZ(size, (std::complex<double> *) x, DimsStride,
|
|
(std::complex<double> *) y, 1, (std::complex<double> *) output, 1);
|
|
}
|
|
|
|
inline double* _MxM_CC_Calc(
|
|
CalcFuncZ aFuncZ, const int size, double* x,double* y,
|
|
int DimsStride) {
|
|
double *output = malloc(size, true);
|
|
V_MxM_CC_Calc(aFuncZ, size, x, y, output, DimsStride);
|
|
return output;
|
|
}
|
|
|
|
inline Matrix operatorMxM(CalcFuncD aFuncD, CalcFuncZ aFuncZ, const Matrix &aMatrix,
|
|
const Matrix &aOther) {
|
|
// 2v2,1v1,3v3
|
|
if (aMatrix.compareShape(aOther)) {
|
|
int DimsStride = 1;
|
|
double *data = nullptr;
|
|
if (aMatrix.getValueType() != aOther.getValueType()) {
|
|
if (aMatrix.getValueType() == Normal) {
|
|
data = _MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
|
|
DimsStride);
|
|
return Matrix::New(data, aOther);
|
|
} else {
|
|
data = _MxM_CN_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),
|
|
DimsStride);
|
|
return Matrix::New(data, aMatrix);
|
|
}
|
|
|
|
} else if (aMatrix.getValueType() == Normal) {
|
|
data = _MxM_NN_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), DimsStride);
|
|
return Matrix::New(data, aMatrix);
|
|
} else {
|
|
data = _MxM_CC_Calc(aFuncZ, aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(), DimsStride);
|
|
return Matrix::New(data, aMatrix);
|
|
}
|
|
}
|
|
//0v3, 0v2
|
|
else if (aMatrix.getDataSize()==1){
|
|
if (aMatrix.getValueType() ==Normal)return operatorMxA(aFuncD,aMatrix.getData()[0],aOther);
|
|
else{
|
|
std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
//3v0, 2v0
|
|
else if (aOther.getDataSize()==1){
|
|
if (aOther.getValueType() ==Normal)return operatorMxA(aFuncD,aOther.getData()[0],aMatrix);
|
|
else{
|
|
std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
//other
|
|
else {
|
|
std::cerr<<"M * M Shape error!If you want do vector * Matrix, use repmat to replicate the vector"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
|
|
inline Matrix operatorMxM_RR(CalcFuncD aFuncD, CalcFuncZ aFuncZ, const Aurora::Matrix &aMatrix,
|
|
Aurora::Matrix &&aOther,bool oppside = false) {
|
|
if (aMatrix.compareShape(aOther)) {
|
|
int DimsStride = 1;
|
|
double* X = oppside?aOther.getData():aMatrix.getData();
|
|
double* 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)
|
|
:_MxM_CN_Calc(aFuncD,aMatrix.getDataSize(), aMatrix.getData(), aOther.getData(),DimsStride);
|
|
return Matrix::New(output, aMatrix);
|
|
}
|
|
//aOther is a complex matrix, use aOther as output
|
|
if (oppside){
|
|
V_MxM_CN_Calc(aFuncD, aMatrix.getDataSize(), aOther.getData(),aMatrix.getData(),aOther.getData(),DimsStride);
|
|
}
|
|
else
|
|
{
|
|
V_MxM_NC_Calc(aFuncD, aMatrix.getDataSize(), aMatrix.getData(),aOther.getData(),aOther.getData(),DimsStride);
|
|
}
|
|
return aOther;
|
|
} else if (aMatrix.getValueType() == Normal) {
|
|
V_MxM_NN_Calc(aFuncD, aMatrix.getDataSize(), X, Y, aOther.getData(),
|
|
DimsStride);
|
|
return aOther;
|
|
} else {
|
|
V_MxM_CC_Calc(aFuncZ, aMatrix.getDataSize(), X, Y, aOther.getData(),
|
|
DimsStride);
|
|
return aOther;
|
|
}
|
|
}
|
|
//0v3, 0v2
|
|
else if (aMatrix.getDataSize()==1){
|
|
if (aMatrix.getValueType() ==Normal){
|
|
return operatorMxA(aFuncD,aMatrix.getData()[0],std::forward<Aurora::Matrix &&>(aOther));
|
|
}
|
|
else{
|
|
std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
//3v0, 2v0
|
|
else if (aOther.getDataSize()==1){
|
|
if (aOther.getValueType() ==Normal)return operatorMxA(aFuncD,aOther.getData()[0],aMatrix);
|
|
else{
|
|
std::cerr<<"M * M fail, Complex scalar * not support now!"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
//other
|
|
else {
|
|
std::cerr<<"M * M Shape error!If you want do vector * Matrix, use repmat to replicate the vector"<<std::endl;
|
|
return Matrix();
|
|
}
|
|
}
|
|
};
|
|
|
|
|
|
namespace Aurora {
|
|
Matrix::Matrix(std::shared_ptr<double> aData, std::vector<int> aInfo, ValueType aValueType)
|
|
: mValueType(aValueType)
|
|
, mData(aData)
|
|
, mInfo(aInfo)
|
|
{
|
|
size_t infoSize = mInfo.size();
|
|
for(;infoSize<3;++infoSize)
|
|
{
|
|
mInfo.push_back(1);
|
|
}
|
|
}
|
|
|
|
Matrix::Matrix(const Matrix::MatrixSlice& slice) {
|
|
auto temp = slice.toMatrix();
|
|
this->mData = temp.mData;
|
|
this->mInfo = temp.mInfo;
|
|
this->mValueType = temp.mValueType;
|
|
}
|
|
|
|
bool Matrix::isNull() const {
|
|
return !mData || mInfo.empty();
|
|
}
|
|
|
|
bool Matrix::isNan() const
|
|
{
|
|
for(size_t i=0; i<getDataSize(); ++i)
|
|
{
|
|
if(mData.get()[i] == mData.get()[i])
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool Matrix::isScalar() const {
|
|
return (getDimSize(0) == 1 &&
|
|
getDimSize(1) == 1 &&
|
|
getDimSize(2) < 2);
|
|
}
|
|
|
|
double Matrix::getScalar() const {
|
|
if (isNull()) return 0.0;
|
|
if (isNull()) return 0.0;
|
|
return getData()[0];
|
|
}
|
|
|
|
bool Matrix::isVector() const{
|
|
if (getDimSize(2)>1) return false;
|
|
if (isScalar()) return false;
|
|
return getDimSize(0) == 1 ||
|
|
getDimSize(1) == 1;
|
|
}
|
|
|
|
bool Matrix::isMKLAllocated() const
|
|
{
|
|
return mMKL_Allocated;
|
|
}
|
|
|
|
int Matrix::getDims() const {
|
|
if(mInfo[2] > 1)
|
|
{
|
|
return 3;
|
|
}
|
|
return 2;
|
|
}
|
|
|
|
double *Matrix::getData() const {
|
|
return mData.get();
|
|
}
|
|
|
|
int Matrix::getDimSize(int aIndex) const {
|
|
if (aIndex >= 0 && aIndex < 3) {
|
|
return mInfo.at(aIndex);
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
|
|
size_t Matrix::getDataSize() const {
|
|
if (!mData.get())return 0;
|
|
size_t ret = 1;
|
|
for (auto v: mInfo) {
|
|
ret *= v;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
void Matrix::forceReshape(int rows, int columns, int slices)
|
|
{
|
|
mInfo = {rows,columns,slices};
|
|
}
|
|
|
|
bool Matrix::compareShape(const Matrix &other) const {
|
|
if (mInfo[2] == 1 && other.mInfo[2] == 1) {
|
|
if (mInfo[0]==1 && other.mInfo[1] == 1 && mInfo[1] == other.mInfo[0]) return true;
|
|
if (mInfo[1]==1 && other.mInfo[0] == 1 && mInfo[0] == other.mInfo[1]) return true;
|
|
}
|
|
for (int i = 0; i < mInfo.size(); ++i) {
|
|
if (mInfo[i] != other.mInfo[i]) return false;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
Matrix Matrix::New(double *data, int rows, int cols, int slices, ValueType type) {
|
|
if (!data) return Matrix();
|
|
std::vector<int> vector(3);
|
|
vector[0]=rows;
|
|
vector[1] = (cols > 0?cols:1);
|
|
vector[2] = (slices > 0?slices:1);
|
|
Matrix ret({data, free}, vector);
|
|
if (type != Normal)ret.setValueType(type);
|
|
ret.mMKL_Allocated = true;
|
|
return ret;
|
|
}
|
|
|
|
Matrix Matrix::New(double *data, const Matrix &shapeMatrix) {
|
|
return New(data,
|
|
shapeMatrix.getDimSize(0),
|
|
shapeMatrix.getDimSize(1),
|
|
shapeMatrix.getDimSize(2),
|
|
shapeMatrix.getValueType());
|
|
}
|
|
|
|
Matrix Matrix::New(const Matrix &shapeMatrix) {
|
|
double *newBuffer = malloc(shapeMatrix.getDataSize(), shapeMatrix.getValueType());
|
|
return New(newBuffer, shapeMatrix);
|
|
}
|
|
|
|
Matrix Matrix::fromRawData(double *data, int rows, int cols, int slices, ValueType type) {
|
|
if (!data) return Matrix();
|
|
std::vector<int> vector;
|
|
vector.push_back(rows);
|
|
if (cols > 0)vector.push_back(cols);
|
|
if (slices > 0)vector.push_back(slices);
|
|
Matrix ret({data, std::default_delete<double[]>()}, vector);
|
|
if (type != Normal)ret.setValueType(type);
|
|
return ret;
|
|
}
|
|
|
|
Matrix Matrix::copyFromRawData(double *data, int rows, int cols, int slices, ValueType type) {
|
|
int colsize = cols>0?cols:1;
|
|
int slicesize = slices>0?slices:1;
|
|
int size = rows*colsize*slicesize;
|
|
double *newBuffer = Aurora::malloc(size, type);
|
|
cblas_dcopy(size*type,data,1,newBuffer,1);
|
|
return New(newBuffer,rows,cols,slices,type);
|
|
}
|
|
|
|
Matrix Matrix::deepCopy() const {
|
|
double *newBuffer = malloc(getDataSize(), getValueType()==Complex);
|
|
cblas_dcopy(getDataSize() * getValueType(),getData(),1,newBuffer,1);
|
|
return New(newBuffer,
|
|
getDimSize(0),
|
|
getDimSize(1),
|
|
getDimSize(2),
|
|
getValueType());
|
|
}
|
|
//operation +
|
|
Matrix Matrix::operator+(double aScalar) const { return operatorMxA(&vdAddI, aScalar, *this);}
|
|
Matrix operator+(double aScalar, const Matrix &matrix) {return matrix + aScalar;}
|
|
Matrix Matrix::operator+(const Matrix &matrix) const {return operatorMxM(vdAddI, vzAddI, *this, matrix);}
|
|
Matrix &operator+(double aScalar, Matrix &&matrix) {
|
|
return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix &operator+(Matrix &&matrix,double aScalar) {
|
|
return operatorMxA_RR(&vdAddI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix Matrix::operator+(Matrix &&aMatrix) const {
|
|
return operatorMxM_RR(&vdAddI,&vzAddI,*this,std::forward<Matrix&&>(aMatrix));
|
|
}
|
|
|
|
Matrix operator+(Matrix &&aMatrix, Matrix &aOther) {
|
|
return operatorMxM_RR(&vdAddI,&vzAddI,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 matrix - aScalar;}
|
|
Matrix Matrix::operator-(const Matrix &matrix) const {return operatorMxM(vdSubI, vzSubI, *this, matrix);}
|
|
Matrix &operator-(double aScalar, Matrix &&matrix) {
|
|
return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix &operator-(Matrix &&matrix,double aScalar) {
|
|
return operatorMxA_RR(&vdSubI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix Matrix::operator-(Matrix &&aMatrix) const {
|
|
return operatorMxM_RR(&vdSubI,&vzSubI,*this,std::forward<Matrix&&>(aMatrix));
|
|
}
|
|
|
|
Matrix operator-(Matrix &&aMatrix, Matrix &aOther) {
|
|
return operatorMxM_RR(&vdSubI,&vzSubI,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*(const Matrix &matrix) const {return operatorMxM(vdMulI, vzMulI, *this, matrix);}
|
|
Matrix &operator*(double aScalar, Matrix &&matrix) {
|
|
return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix &operator*(Matrix &&matrix,double aScalar) {
|
|
return operatorMxA_RR(&vdMulI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix Matrix::operator*(Matrix &&aMatrix) const {
|
|
return operatorMxM_RR(&vdMulI,&vzMulI,*this,std::forward<Matrix&&>(aMatrix));
|
|
}
|
|
|
|
Matrix operator*(Matrix &&aMatrix, Matrix &aOther) {
|
|
return operatorMxM_RR(&vdMulI,&vzMulI,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 matrix / aScalar;}
|
|
Matrix Matrix::operator/(const Matrix &matrix) const {return operatorMxM(vdDivI, vzDivI, *this, matrix);}
|
|
Matrix &operator/(double aScalar, Matrix &&matrix) {
|
|
return operatorMxA_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix &operator/(Matrix &&matrix,double aScalar) {
|
|
return operatorMxA_RR(&vdDivI,aScalar, std::forward<Matrix&&>(matrix));
|
|
}
|
|
Matrix Matrix::operator/(Matrix &&aMatrix) const {
|
|
return operatorMxM_RR(&vdDivI,&vzDivI,*this,std::forward<Matrix&&>(aMatrix));
|
|
}
|
|
|
|
Matrix operator/(Matrix &&aMatrix, Matrix &aOther) {
|
|
return operatorMxM_RR(&vdDivI,&vzDivI,aOther,std::forward<Matrix&&>(aMatrix),true);
|
|
}
|
|
//operator ^ (pow)
|
|
Matrix Matrix::operator^(int times) const { return operatorMxA(&vdPowI, times, *this);}
|
|
|
|
Matrix operator^( Matrix &&matrix,int times) {
|
|
return operatorMxA(&vdPowI, times, std::forward<Matrix&&>(matrix));
|
|
}
|
|
void Matrix::printf() {
|
|
if(isNull())
|
|
{
|
|
::printf("null matrix\n");
|
|
return;
|
|
}
|
|
int k_count = getDimSize(2);
|
|
int j_count = getDimSize(1);
|
|
int complexstep = 1;
|
|
const char* mark = "+";
|
|
if(mValueType == Complex)
|
|
{
|
|
complexstep = 2;
|
|
}
|
|
for (int k = 0; k <k_count; ++k) {
|
|
::printf("slice %d:\r\n[",k);
|
|
for (int i = 0; i < getDimSize(0); ++i) {
|
|
::printf("[");
|
|
for (int j = 0; j < j_count; ++j) {
|
|
::printf("%f2",getData()[(k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i)*complexstep]);
|
|
if(mValueType == Complex)
|
|
{
|
|
double value = getData()[(k*getDimSize(1)*getDimSize(0)+j*getDimSize(0)+i)*complexstep+1];
|
|
if(value<0)
|
|
{
|
|
mark = "";
|
|
}
|
|
::printf("%s%f2i",mark,value);
|
|
mark = "+";
|
|
}
|
|
::printf(", ");
|
|
}
|
|
::printf("]\r\n");
|
|
}
|
|
::printf("]\r\n");
|
|
}
|
|
}
|
|
|
|
void Matrix::printfShape() {
|
|
std::cerr << "Matrix shape:(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
}
|
|
|
|
Matrix::MatrixSlice Matrix::operator()(int aRowIdx, int aColIdx, int aSliceIdx) const {
|
|
std::vector<int> dims = {aRowIdx, aColIdx, aSliceIdx};
|
|
std::vector<int> allDimIndex;
|
|
int mode = 0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
if (dims[j]==$ && getDims()>j){
|
|
++mode;
|
|
allDimIndex.push_back(j);
|
|
}
|
|
}
|
|
int rowStride = 1;
|
|
int colStride = getDimSize(0);
|
|
int sliceStride = getDimSize(0)*getDimSize(1);
|
|
int strides[3] = {rowStride, colStride, sliceStride};
|
|
int rowOffset = aRowIdx == $ ? 0 : aRowIdx;
|
|
int colOffset = aColIdx == $ ? 0 : aColIdx;
|
|
int sliceOffset = aSliceIdx == $ ? 0 : aSliceIdx;
|
|
double *startPointer = getData() + (rowStride * rowOffset
|
|
+ colStride * colOffset
|
|
+ sliceStride * sliceOffset) * getValueType();
|
|
int size1 = allDimIndex.empty()?1:getDimSize(allDimIndex[0]);
|
|
int stride1 = allDimIndex.empty()?1:strides[allDimIndex[0]];
|
|
switch (mode) {
|
|
//matrix mode
|
|
case 2:{
|
|
int size2 = getDimSize(allDimIndex[1]);
|
|
int stride2 = strides[allDimIndex[1]];
|
|
return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode, size2, stride2);
|
|
}
|
|
//vector mode & default
|
|
case 1:
|
|
{
|
|
return Matrix::MatrixSlice(size1, stride1, startPointer, getValueType(), mode);
|
|
}
|
|
//scalar mode or ALL $
|
|
case 0:
|
|
default: {
|
|
return Matrix::MatrixSlice(1 , 1, startPointer,getValueType(), 0);
|
|
}
|
|
}
|
|
}
|
|
|
|
Matrix Matrix::operator>(double 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());
|
|
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) {
|
|
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());
|
|
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));
|
|
}
|
|
|
|
Matrix Matrix::operator>(const Matrix &matrix) const {
|
|
if (isComplex() || matrix.isComplex()) {
|
|
std::cerr << "Complex cann't compare!" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if (!compareShape(matrix) && !isScalar() && !matrix.isScalar()) {
|
|
std::cerr << "Matrix not equal, matrix 1(" << matrix.getDimSize(0) << "," << matrix.getDimSize(1) << ","
|
|
<< matrix.getDimSize(2) << "), matrix 2(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if(isScalar()){
|
|
return getData()[0]>matrix;
|
|
}
|
|
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());
|
|
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 {
|
|
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());
|
|
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) {
|
|
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());
|
|
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));
|
|
}
|
|
|
|
Matrix Matrix::operator>=(const Matrix &matrix) const {
|
|
if (isComplex() || matrix.isComplex()) {
|
|
std::cerr << "Complex cann't compare!" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if (!compareShape(matrix) && !isScalar() && !matrix.isScalar()) {
|
|
std::cerr << "Matrix not equal, matrix 1(" << matrix.getDimSize(0) << "," << matrix.getDimSize(1) << ","
|
|
<< matrix.getDimSize(2) << "), matrix 2(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if(isScalar()){
|
|
return getData()[0]<=matrix;
|
|
}
|
|
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());
|
|
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 {
|
|
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());
|
|
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) {
|
|
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());
|
|
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));
|
|
}
|
|
|
|
Matrix Matrix::operator<(const Matrix &matrix) const {
|
|
if (isComplex() || matrix.isComplex()) {
|
|
std::cerr << "Complex cann't compare!" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if (!compareShape(matrix) && !isScalar() && !matrix.isScalar()) {
|
|
std::cerr << "Matrix not equal, matrix 1(" << matrix.getDimSize(0) << "," << matrix.getDimSize(1) << ","
|
|
<< matrix.getDimSize(2) << "), matrix 2(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if(isScalar()){
|
|
return getData()[0]<matrix;
|
|
}
|
|
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());
|
|
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 {
|
|
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());
|
|
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) {
|
|
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());
|
|
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));
|
|
}
|
|
|
|
Matrix Matrix::operator<=(const Matrix &matrix) const {
|
|
if (isComplex() || matrix.isComplex()) {
|
|
std::cerr << "Complex cann't compare!" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if (!compareShape(matrix) && !isScalar() && !matrix.isScalar()) {
|
|
std::cerr << "Matrix not equal, matrix 1(" << matrix.getDimSize(0) << "," << matrix.getDimSize(1) << ","
|
|
<< matrix.getDimSize(2) << "), matrix 2(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if(isScalar()){
|
|
return getData()[0]<=matrix;
|
|
}
|
|
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());
|
|
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 {
|
|
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());
|
|
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) {
|
|
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());
|
|
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));
|
|
}
|
|
|
|
Matrix Matrix::operator==(const Matrix &matrix) const {
|
|
if (isComplex() || matrix.isComplex()) {
|
|
std::cerr << "Complex cann't compare!" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if (!compareShape(matrix) && !isScalar() && !matrix.isScalar()) {
|
|
std::cerr << "Matrix not equal, matrix 1(" << matrix.getDimSize(0) << "," << matrix.getDimSize(1) << ","
|
|
<< matrix.getDimSize(2) << "), matrix 2(" << getDimSize(0) << "," << getDimSize(1) << ","
|
|
<< getDimSize(2) << ")" << std::endl;
|
|
return Matrix();
|
|
}
|
|
if(isScalar()){
|
|
return getData()[0]==matrix;
|
|
}
|
|
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());
|
|
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 operator-(Matrix &&aMatrix) {
|
|
int size = aMatrix.getDataSize()*aMatrix.getValueType();
|
|
double zero = 0.0;
|
|
vdSubI(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;
|
|
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);
|
|
}
|
|
else{
|
|
vdSubI( aMatrix.getDataSize(),&zero,0,aMatrix.getData(),1,newBuffer,1);
|
|
}
|
|
|
|
return Matrix::New(newBuffer,aMatrix);
|
|
}
|
|
|
|
Matrix::MatrixSlice::MatrixSlice(int aSize,int aStride, double* aData, ValueType aType, int aSliceMode,int aSize2, int aStride2):
|
|
mSliceMode(aSliceMode),mData(aData),
|
|
mSize(aSize),mSize2(aSize2),
|
|
mStride(aStride),mStride2(aStride2),
|
|
mType(aType)
|
|
{
|
|
|
|
}
|
|
|
|
Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix::MatrixSlice &slice) {
|
|
if (this==&slice) return *this;
|
|
if(!mData){
|
|
std::cerr <<"Assign value fail!Des data pointer is null!";
|
|
return *this;
|
|
}
|
|
if(!slice.mData){
|
|
std::cerr <<"Assign value fail!Src data pointer is null!"<<std::endl;
|
|
return *this;
|
|
}
|
|
if (slice.mSliceMode!=mSliceMode) {
|
|
std::cerr << "Assign value fail!Src slice(dims count:" << slice.mSliceMode
|
|
<< "), not match of des(dims count:" << mSliceMode << ")!" << std::endl;
|
|
return *this;
|
|
}
|
|
if (slice.mSize!=mSize) {
|
|
std::cerr <<"Assign value fail!Src slice(dim 1 size:"<< slice.mSize <<"), not match of des(dim 1 size:"<<mSize<<")!"<<std::endl;
|
|
return *this;
|
|
}
|
|
if (slice.mSize2!=mSize2) {
|
|
std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< slice.mSize2 <<"), not match of des(dim 2 size:"<<mSize2<<")!"<<std::endl;
|
|
return *this;
|
|
}
|
|
if (slice.mType!=mType) {
|
|
std::cerr <<"Assign value fail!Src slice(value type:"<< slice.mType <<"), not match of des(value type:"<<mType<<")!"<<std::endl;
|
|
return *this;
|
|
}
|
|
switch (mSliceMode) {
|
|
case 2:{
|
|
if (mType== Normal) {
|
|
// cblas_dcopy_batch_strided(mSize, slice.mData, slice.mStride, slice.mStride2, mData, mStride,
|
|
// mStride2, mSize2);
|
|
for (int i = 0; i < mSize2; ++i) {
|
|
cblas_dcopy(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);
|
|
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);
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
case 1:{
|
|
if (mType== Normal){
|
|
cblas_dcopy(mSize,slice.mData,slice.mStride,mData,mStride);
|
|
}
|
|
else {
|
|
cblas_dcopy(mSize*2, slice.mData, slice.mStride, mData, mStride);
|
|
}
|
|
break;
|
|
}
|
|
case 0:
|
|
default:{
|
|
mData[0] = slice.mData[0];
|
|
if (mType != Normal)mData[1] = slice.mData[1];
|
|
}
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(const Matrix &matrix) {
|
|
if(!mData){
|
|
std::cerr <<"Assign value fail!Des data pointer is null!";
|
|
return *this;
|
|
}
|
|
if(!matrix.getData()){
|
|
std::cerr <<"Assign value fail!Src data pointer is null!";
|
|
return *this;
|
|
}
|
|
if (matrix.getDims()!=mSliceMode) {
|
|
std::cerr <<"Assign value fail!Src matrix(dims count:"<< matrix.getDims() <<"), not match of des(dims count:"<<mSliceMode<<")!";
|
|
return *this;
|
|
}
|
|
if (matrix.getDimSize(0)!=mSize) {
|
|
std::cerr <<"Assign value fail!Src matrix(dim 1 size:"<< matrix.getDimSize(0)<<"), not match of des(dim 1 size:"<<mSize<<")!";
|
|
return *this;
|
|
}
|
|
if (matrix.getDimSize(1)!=mSize2) {
|
|
std::cerr <<"Assign value fail!Src slice(dim 2 size:"<< matrix.getDimSize(1) <<"), not match of des(dim 2 size:"<<mSize2<<")!";
|
|
return *this;
|
|
}
|
|
if (matrix.getValueType()!=mType) {
|
|
std::cerr <<"Assign value fail!Src slice(value type:"<< matrix.getValueType() <<"), not match of des(value type:"<<mType<<")!";
|
|
return *this;
|
|
}
|
|
switch (mSliceMode) {
|
|
case 2:{
|
|
if (mType== Normal) {
|
|
cblas_dcopy_batch_strided(mSize, matrix.getData(), 1, matrix.getDimSize(0), mData, mStride,
|
|
mStride2, mSize2);
|
|
}
|
|
else {
|
|
cblas_zcopy_batch_strided(mSize,(std::complex<double>*)matrix.getData(),1,matrix.getDimSize(0),
|
|
(std::complex<double>*)mData,mStride,mStride2,mSize2);
|
|
}
|
|
break;
|
|
}
|
|
case 1:{
|
|
if (mType== Normal){
|
|
cblas_dcopy(mSize,matrix.getData(),1,mData,mStride);
|
|
}
|
|
else {
|
|
cblas_zcopy(mSize, (std::complex<double> *) matrix.getData(),1,
|
|
(std::complex<double> *) mData, mStride);
|
|
}
|
|
break;
|
|
}
|
|
case 0:
|
|
default:{
|
|
mData[0] = matrix.getData()[0];
|
|
if (mType != Normal)mData[1] = matrix.getData()[1];
|
|
}
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(double value) {
|
|
if(!mData){
|
|
std::cerr <<"Assign value fail!Des data pointer is null!";
|
|
return *this;
|
|
}
|
|
if (mSliceMode!=0) {
|
|
std::cerr <<"Assign value fail!Des slicemode is"<<mSliceMode<<", not scalar mode!";
|
|
return *this;
|
|
}
|
|
if (mSize!=1) {
|
|
std::cerr <<"Assign value fail!Des size:"<<mSize<<", not scalar mode!";
|
|
return *this;
|
|
}
|
|
if (mType!=Normal) {
|
|
std::cerr <<"Assign value fail!Des type is complex!";
|
|
return *this;
|
|
}
|
|
mData[0]=value;
|
|
return *this;
|
|
}
|
|
|
|
Matrix::MatrixSlice &Matrix::MatrixSlice::operator=(std::complex<double> value) {
|
|
if(!mData){
|
|
std::cerr <<"Assign value fail!Des data pointer is null!";
|
|
return *this;
|
|
}
|
|
if (mSliceMode!=0) {
|
|
std::cerr <<"Assign value fail!Des slicemode is"<<mSliceMode<<", not scalar mode!";
|
|
return *this;
|
|
}
|
|
if (mSize!=1) {
|
|
std::cerr <<"Assign value fail!Des size:"<<mSize<<", not scalar mode!";
|
|
return *this;
|
|
}
|
|
if (mType!=Complex) {
|
|
std::cerr <<"Assign value fail!Des type is not complex!";
|
|
return *this;
|
|
}
|
|
mData[0]=value.real();
|
|
mData[1]=value.imag();
|
|
return *this;
|
|
}
|
|
|
|
Matrix Matrix::MatrixSlice::toMatrix() const {
|
|
double * data = (double *) malloc(mSize*(mSize2>0?mSize2:1) ,mType == Complex);
|
|
|
|
switch (mSliceMode) {
|
|
case 2:{
|
|
if (mType== Normal) {
|
|
cblas_dcopy_batch_strided(mSize, mData, mStride,
|
|
mStride2,data, 1, mSize, mSize2);
|
|
}
|
|
else {
|
|
cblas_zcopy_batch_strided(mSize, (std::complex<double> *) mData, mStride, mStride2,
|
|
(std::complex<double> *) data, 1, mSize,
|
|
mSize2);
|
|
}
|
|
break;
|
|
}
|
|
case 1:{
|
|
if (mType== Normal){
|
|
cblas_dcopy(mSize,mData,mStride,data,1);
|
|
}
|
|
else {
|
|
cblas_zcopy(mSize, (std::complex<double> *) mData, mStride,
|
|
(std::complex<double> *) data, 1);
|
|
}
|
|
break;
|
|
}
|
|
case 0:
|
|
default:{
|
|
data[0]= mData[0];
|
|
if (mType != Normal) data[1] = mData[1];
|
|
}
|
|
}
|
|
|
|
return Matrix::New(data,mSize,mSize2,1,mType);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|