sunwen/Aurora
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add vecnorm and vecnorm's unittest.

Browse Source 
Add norm function with complex support.
This commit is contained in:
sunwen
2023-04-26 11:17:15 +08:00
parent d1389df8a4
commit 7693ce9f9e
3 changed files with 98 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
src/Function1D.cpp
View File

@@ -1,4 +1,5 @@
#include "Function1D.h"
#include "Function2D.h"
#include "Function.h"
//必须在Eigen之前
@@ -420,12 +421,16 @@ Matrix Aurora::conj(const Matrix& aMatrix)
double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
{
if(aMatrix.isComplex() || aMatrix.isNull())
if(aMatrix.isNull())
{
return NAN;
}
size_t size = aMatrix.getDataSize();
if(aMatrix.isComplex())
{
size*=2;
}
int column = aMatrix.getDimSize(1);
int row = aMatrix.getDimSize(0);
if (aNormMethod == NormMethod::Norm1)
@@ -433,7 +438,7 @@ double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
double value = 0;
for(int i=0; i<column; ++i)
{
double temp = cblas_dasum(row, aMatrix($,i,$).toMatrix().getData(), 1);
double temp = Aurora::sum(abs(aMatrix($,i,$).toMatrix())).getData()[0];
if(temp > value)
{
value = temp;
@@ -450,9 +455,18 @@ double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
//columns > 1
if(aMatrix.getDimSize(1) > 1)
{
Eigen::Map<Eigen::MatrixXd> eMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
Eigen::JacobiSVD<Eigen::MatrixXd> svd(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
return svd.singularValues()(0);
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);
}
else
{
Eigen::Map<Eigen::MatrixXd> eMatrix(aMatrix.getData(), row, column);
Eigen::JacobiSVD<Eigen::MatrixXd> svd(eMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
return svd.singularValues()(0);
}
}
else
{
@@ -507,3 +521,20 @@ Matrix Aurora::horzcat(const Matrix& aMatrix1, const Matrix& aMatrix2)
return Matrix::New(resultData, row, column1+column2, 1, aMatrix1.getValueType());
}
Matrix Aurora::vecnorm(const Matrix& aMatrix, NormMethod aNormMethod, int aDim)
{
//only surpport aDim = 1 for now.
if(aDim != 1 || aNormMethod == NormMethod::NormF || aMatrix.isNull())
{
return Matrix();
}
int column = aMatrix.getDimSize(1);
double* resultData = Aurora::malloc(column);
for(int i=0; i<column; ++i)
{
resultData[i] = norm(aMatrix($,i,$).toMatrix(), aNormMethod);
}
return Matrix::New(resultData,column);
}

2
src/Function1D.h
View File

@@ -70,6 +70,8 @@ namespace Aurora {
Matrix horzcat(const Matrix& aMatrix1, const Matrix& aMatrix2);
Matrix vecnorm(const Matrix& aMatrix, NormMethod aNormMethod, int aDim);
/**
* 多项式计算
* @brief 例如p[1 0 1],x[3 2 5],代表对多项式 y = x^2 + 1 求(x=3, x=2, x=5)时所有的y