UT and annotations for math_bspline

source/module_base/math_bspline.cpp

@@ -32,18 +32,12 @@ namespace ModuleBase
         }
     }
 
-    void Bspline::cleanp()
-    {
-        delete[] bezier;
-        bezier = NULL;
-    }
-
     double Bspline::bezier_ele(int n)
     {
         return this->bezier[n];
     }
 
-    void Bspline::getbslpine(double x)
+    void Bspline::getbspline(double x)
     {
         bezier[0] = 1.0;
         for(int k = 1 ; k <= norder ; ++k)

source/module_base/math_bspline.h

@@ -4,57 +4,54 @@
 namespace ModuleBase
 {
 
-
-//
-//DESCRIPTION:
-//   A class to treat Cardinal B-spline interpolation.
-//   qianrui created 2021-09-14
-//MATH:
-//   Only uniform nodes are considered: xm-x[m-1]=Dx(>= 0) for control node: X={x0,x1,...,xm};
-//   Any function p(x) can be written by
-//   p(x)=\sum_i{ci*M_ik(x)} (k->infinity), 
-//   where M_ik is the i-th k-order Cardinal B-spline base function 
-//   and ci is undetermined coefficient.
-//      M_i0 =  H(x-xi)-H(x-x[i+1]), H(x): step function
-//                  x-xi                    x[i+k+1]-x
-//      M_ik(x)=  ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x)  ( xi <= x <= x[i+1] )
-//                x[i+k]-xi               x[i+k+1]-x[i+1]
-//   For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x)
-//   If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx):
-//               x+n*Dx-xi               xi+(k-n+1)*Dx-x
-//      Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1]
-//                 k*Dx                        k*Dx
-//USAGE：
-//   ModuleBase::Bspline bp;
-//   bp.init(10,0.7,2); //Dx = 0.7, xi = 2
-//   bp.getbslpine(0.5); //x = 0.5
-//   cout<<bp.bspline_ele(3)<<endl; //print M_ik[xi+3*Dx+x]: M_i[10](4.6) 
-//
+/**
+ * @brief A class to treat Cardinal B-spline interpolation.
+ *
+ * @author qianrui created 2021-09-14
+ * @details see: J. Chem. Phys. 103, 8577 (1995).
+ * Math:
+ * Only uniform nodes are considered: xm-x[m-1]=Dx(>= 0) for control node: X={x0,x1,...,xm};
+ * Any function p(x) can be written by
+ * p(x)=\sum_i{ci*M_ik(x)} (k->infinity),
+ * where M_ik is the i-th k-order Cardinal B-spline base function
+ * and ci is undetermined coefficient.
+ * M_i0 =  H(x-xi)-H(x-x[i+1]), H(x): step function
+ *                  x-xi                    x[i+k+1]-x
+ *      M_ik(x)=  ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x)  ( xi <= x <= x[i+1] )
+ *                x[i+k]-xi               x[i+k+1]-x[i+1]
+ *   For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x)
+ *   If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx):
+ *               x+n*Dx-xi               xi+(k-n+1)*Dx-x
+ *      Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1]
+ *                 k*Dx                        k*Dx
+ * USAGE：
+ *   ModuleBase::Bspline bp;
+ *   bp.init(10,0.7,2); //Dx = 0.7, xi = 2
+ *   bp.getbslpine(0.5); //x = 0.5
+ *   cout<<bp.bezier_ele(3)<<endl; //print M_ik[xi+3*Dx+x]: M_i[10](4.6)
+ *
+ */
 class Bspline
 {
-    private:
-        int norder; // the order of bezier base; norder >= 0
-        double Dx; //Dx: the interval of control node
-        double xi; // xi: the starting point
-        double *bezier; //bezier[n] = Bk[n]
+  private:
+    int norder; // the order of bezier base; norder >= 0
+    double Dx; // Dx: the interval of control node
+    double xi; // xi: the starting point
+    double *bezier; // bezier[n] = Bk[n]
 
-    public:
-        Bspline();
-        ~Bspline();
+  public:
+    Bspline();
+    ~Bspline();
 
-        //Init norder, Dx, xi
-        void init(int norderin, double Dxin, double xiin);
+    void init(int norderin, double Dxin, double xiin);
 
-        //delete[] bezier
-        void cleanp();
+    // Get the result of i-th bezier base functions for different input x+xi+n*Dx.
+    // x should be in [0,Dx]
+    // n-th result is stored in bezier[n];
+    void getbspline(double x);
 
-        //Get the result of i-th bezier base functions for different input x+xi+n*Dx.
-        //x should be in [0,Dx]
-        //n-th result is stored in bezier[n];
-        void getbslpine(double x);
-
-        //get the element of bezier
-        double bezier_ele(int n);
+    // get the element of bezier
+    double bezier_ele(int n);
 };
-}
+} // namespace ModuleBase
 #endif

source/module_base/test/CMakeLists.txt

@@ -82,3 +82,7 @@ AddTest(
   LIBS ${math_libs}
   SOURCES math_ylmreal_test.cpp ../math_ylmreal.cpp ../realarray.cpp ../timer.cpp ../matrix.cpp 
 )
+AddTest(
+  TARGET base_math_bspline
+  SOURCES math_bspline_test.cpp ../math_bspline.cpp
+)

source/module_base/test/math_bspline_test.cpp

@@ -0,0 +1,62 @@
+#include "../math_bspline.h"
+#include "gtest/gtest.h"
+
+/************************************************
+ *  unit test of class Bspline
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - Init
+ *     - norder must be even
+ *     - norder mush be positive
+ *   - Properties
+ *     - \sum_i M_n(u+i) = 1 (i=0,1,2,...n)
+ *
+ */
+
+class MathBsplineTest : public testing::Test
+{
+protected:
+	ModuleBase::Bspline bp;
+	int norder;
+};
+
+TEST_F(MathBsplineTest,Init)
+{
+	EXPECT_DEATH(
+		{
+			norder = 3; // norder must be even
+			bp.init(norder,0.05,0);
+		},""
+	);
+	EXPECT_DEATH(
+		{
+			norder = 0; // norder must be positive
+			bp.init(norder,0.05,0);
+		},""
+	);
+}
+
+// summation over n is unity
+TEST_F(MathBsplineTest,Properties)
+{
+	int by = 2;
+	for (norder=2;norder<=20;norder=norder+by)
+	{
+		bp.init(norder,1.0,0);
+		bp.getbspline(0.2);
+		double sum=0.0;
+		//std::cout << "\n" << "norder : "<< norder<<std::endl;
+		for (int i=0;i<=norder;i++)
+		{
+			//std::cout << i << " " << bp.bezier_ele(i) << std::endl;
+			sum += bp.bezier_ele(i);
+		}
+		//std::cout<<"sum "<< sum<< std::endl;
+		EXPECT_NEAR(sum,1.0,1.e-15);
+		bp.getbspline(0.0); // must be set to 0.0 to test M_n(0)
+		EXPECT_NEAR(bp.bezier_ele(0),0,1.e-16); 
+		EXPECT_NEAR(bp.bezier_ele(norder+1),0,1.e-16); 
+	}
+}

source/src_pw/bspline_sf.cpp

@@ -46,9 +46,9 @@ void PW_Basis::bspline_sf(const int norder)
             bsx.init(norder, 1, 0);
             bsy.init(norder, 1, 0);
             bsz.init(norder, 1, 0);
-            bsx.getbslpine(dx);
-            bsy.getbslpine(dy);
-            bsz.getbslpine(dz);
+            bsx.getbspline(dx);
+            bsy.getbspline(dy);
+            bsz.getbspline(dz);
 
             for(int iz = 0 ; iz <= norder ; ++iz)
             {
@@ -109,7 +109,7 @@ void PW_Basis:: bsplinecoef(complex<double> *b1, complex<double> *b2, complex<do
     const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;
     ModuleBase::Bspline bsp;
     bsp.init(norder, 1, 0);
-    bsp.getbslpine(1.0);
+    bsp.getbspline(1.0);
     for(int ix = 0 ; ix < this->nx ; ++ix)
     {
         complex<double> fracx=0;
@@ -137,4 +137,4 @@ void PW_Basis:: bsplinecoef(complex<double> *b1, complex<double> *b2, complex<do
         }
         b3[iz] = exp(ci_tpi*double(norder*iz)/double(nz))/fracz;
     }
-}
+}