From c54471ef6ab108a6df7015438dd8bf124790115b Mon Sep 17 00:00:00 2001
From: sunwen <sunwen@realcan.cn>
Date: Tue, 25 Apr 2023 11:24:24 +0800
Subject: [PATCH] Add norm and norm's unittest.

---
 src/Function1D.cpp       | 47 ++++++++++++++++++++++++++++++++++++++++
 src/Function1D.h         |  7 ++++++
 test/Function1D_Test.cpp | 20 +++++++++++++++++
 3 files changed, 74 insertions(+)

diff --git a/src/Function1D.cpp b/src/Function1D.cpp
index e06e33c..1699233 100644
--- a/src/Function1D.cpp
+++ b/src/Function1D.cpp
@@ -11,6 +11,8 @@
 #include <Eigen/Core>
 #include <Eigen/Eigen>
 #include <Eigen/Dense>
+#include <Eigen/SVD>
+#include <mkl_lapack.h>
 
 using namespace Aurora;
 
@@ -415,3 +417,48 @@ Matrix Aurora::conj(const Matrix& aMatrix)
     vzConj(size,(MKL_Complex16*)aMatrix.getData(), (MKL_Complex16*)data);
     return Matrix::New(data, aMatrix);
 }
+
+double Aurora::norm(const Matrix& aMatrix, NormMethod aNormMethod)
+{
+    if(aMatrix.isComplex() || aMatrix.isNull())
+    {
+        return NAN;
+    }
+
+    size_t size = aMatrix.getDataSize();
+    int column = aMatrix.getDimSize(1);
+    int row = aMatrix.getDimSize(0);
+    if (aNormMethod == NormMethod::Norm1)
+    {
+        double value = 0;
+        for(int i=0; i<column; ++i)
+        {
+            double temp = cblas_dasum(row, aMatrix($,i,$).toMatrix().getData(), 1);
+            if(temp > value)
+            {
+                value = temp;
+            }
+        }
+        return value;
+    }
+    else if(aNormMethod == NormMethod::NormF)
+    {
+        return cblas_dnrm2(size, aMatrix.getData(), 1);
+    }
+    else if(aNormMethod == NormMethod::Norm2)
+    {
+        //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);
+        }
+        else
+        {
+            return cblas_dnrm2(size, aMatrix.getData(), 1);
+        }
+    }
+
+    return 0;
+}
diff --git a/src/Function1D.h b/src/Function1D.h
index 076dc80..814cb5d 100644
--- a/src/Function1D.h
+++ b/src/Function1D.h
@@ -10,6 +10,11 @@ namespace Aurora {
         Spline=0,Linear
     };
 
+    enum NormMethod
+    {
+        Norm1=1,Norm2,NormF
+    };
+
     Matrix complex(const Matrix& matrix);
 
     Matrix real(const Matrix& matrix);
@@ -59,6 +64,8 @@ namespace Aurora {
 
     Matrix conj(const Matrix& aMatrix);
 
+    double norm(const Matrix& aMatrix, NormMethod aNormMethod);
+
     /**
      * 多项式计算
      * @brief 例如p[1 0 1],x[3 2 5],代表对多项式 y = x^2 + 1 求(x=3, x=2, x=5)时所有的y
diff --git a/test/Function1D_Test.cpp b/test/Function1D_Test.cpp
index 85e2d1d..34242a0 100644
--- a/test/Function1D_Test.cpp
+++ b/test/Function1D_Test.cpp
@@ -288,3 +288,23 @@ TEST_F(Function1D_Test, conj) {
     EXPECT_DOUBLE_AE(result.getData()[5],2);
     EXPECT_DOUBLE_AE((double)result.getValueType(), (double)Aurora::Normal);
 }
+
+TEST_F(Function1D_Test, norm) {
+    double *data = new double[3]{1,2,-3};
+    auto matrix = Aurora::Matrix::fromRawData(data, 3);
+    auto result = Aurora::norm(matrix,Aurora::NormMethod::Norm1);
+    EXPECT_DOUBLE_AE(result,6);
+    result = Aurora::norm(matrix,Aurora::NormMethod::Norm2);
+    EXPECT_DOUBLE_AE(result,3.74166);
+    result = Aurora::norm(matrix,Aurora::NormMethod::NormF);
+    EXPECT_DOUBLE_AE(result,3.74166);
+
+    data = new double[8]{1,2,-3,6,7,9,22.3,-8.6};
+    matrix = Aurora::Matrix::fromRawData(data, 4,2);
+    result = Aurora::norm(matrix,Aurora::NormMethod::Norm1);
+    EXPECT_DOUBLE_AE(result,46.9);
+    result = Aurora::norm(matrix,Aurora::NormMethod::Norm2);
+    EXPECT_DOUBLE_AE(result,26.7284);
+    result = Aurora::norm(matrix,Aurora::NormMethod::NormF);
+    EXPECT_DOUBLE_AE(result,27.4089);
+}