<Test> Add unit tests for plain, pulay, and broyden mixing.

source/module_base/CMakeLists.txt

@@ -75,6 +75,7 @@ if(BUILD_TESTING)
   add_subdirectory(test)
   add_subdirectory(test_parallel)
   add_subdirectory(kernels/test)
+  add_subdirectory(module_mixing/test)
   if (USE_ABACUS_LIBM)
     add_subdirectory(libm/test)
   endif()

source/module_base/module_mixing/pulay_mixing.h

@@ -134,13 +134,9 @@ class Pulay_Mixing : public Mixing
                 "One Mixing object can only bind one Mixing_Data object to calculate coefficients");
         const int length = mdata.length;
         FPTYPE* FP_F = static_cast<FPTYPE*>(F);
-        // std::cout << 1 << std::endl;
-        std::cout << "ndim_use    " << mdata.ndim_use << std::endl;
-        std::cout << "ndim_tot    " << mdata.ndim_tot << std::endl;
 
         if (mdata.ndim_use > 1)
         {
-		    std::cout << "beta" << std::endl;
             const int ndim_use = mdata.ndim_use;
             ModuleBase::matrix beta_tmp(ndim_use, ndim_use);
             //beta(i, j) = <F_i, F_j>
@@ -155,17 +151,12 @@ class Pulay_Mixing : public Mixing
                     }
                     FPTYPE* Fj = FP_F + j * length;
                     beta(i, j) = beta_tmp(i, j) = inner_dot(Fi, Fj);
-				    std::cout << beta(i,j) << "\t";
                     if (j != i)
                     {
                         beta(j, i) = beta_tmp(j, i) = beta_tmp(i, j);
                     }
                 }
-			    std::cout << "\n";
             }
-			std::cout << "\n";
-
-        // std::cout << 2 << std::endl;
 
             double* work = new double[ndim_use];
             int* iwork = new int[ndim_use];
@@ -185,17 +176,6 @@ class Pulay_Mixing : public Mixing
                 }
             }
 
-		    std::cout << "beta^-1" << std::endl;
-            for (int i = 0; i < ndim_use; ++i)
-            {
-                for (int j = i; j < ndim_use; ++j)
-                {
-				    std::cout << beta_tmp(i,j) << "\t";
-                }
-			    std::cout << "\n";
-            }
-			std::cout << "\n";
-
             // coef{i} = \sum_j beta{ij} / \sum_k \sum_j beta{kj} 
             double sum_beta = 0.;
             for (int i = 0; i < ndim_use; ++i)
@@ -205,9 +185,6 @@ class Pulay_Mixing : public Mixing
                     sum_beta += beta_tmp(j, i);
                 }
             }
-            std::cout << "sum_beta    " << sum_beta << std::endl;
-
-    		std::cout << "alpha" << std::endl;
             for (int i = 0; i < ndim_use; ++i)
             {
                 coef[i] = 0.;
@@ -216,13 +193,9 @@ class Pulay_Mixing : public Mixing
                     coef[i] += beta_tmp(i,j);
                 }
                 coef[i] /= sum_beta;
-                std::cout << coef[i] << "\t";
             }
-			std::cout << "\n";
             delete[] work;
             delete[] iwork;
-        // std::cout << 3 << std::endl;
-
         }
         else
         {

source/module_base/module_mixing/test/CMakeLists.txt

@@ -0,0 +1,16 @@
+remove_definitions(-D__MPI)
+AddTest(
+  TARGET base_broyden_mixing
+  LIBS base device ${math_libs}
+  SOURCES broyden_mixing_test.cpp
+)
+AddTest(
+  TARGET base_pulay_mixing
+  LIBS base device ${math_libs}
+  SOURCES pulay_mixing_test.cpp
+)
+AddTest(
+  TARGET base_plain_mixing
+  LIBS ${math_libs}
+  SOURCES plain_mixing_test.cpp ../mixing_data.cpp
+)

source/module_base/module_mixing/test/broyden_mixing_test.cpp

@@ -0,0 +1,97 @@
+#include "gtest/gtest.h"
+#include "../broyden_mixing.h"
+
+#define DOUBLETHRESHOLD 1e-8
+
+class Broyden_Mixing_Test : public testing::Test
+{
+protected:
+    const double mixing_beta = 0.8;
+    const int mixing_ndim = 7;
+    Base_Mixing::Mixing_Data xdata;
+    Base_Mixing::Broyden_Mixing* broyden = nullptr;
+
+    double thr = 1e-5;
+    int niter = 0;
+    int maxiter = 50;
+    double* x_out = nullptr;
+
+    void SetUp()
+    {
+        this->broyden = new Base_Mixing::Broyden_Mixing(this->mixing_ndim, this->mixing_beta);
+    }
+    void TearDown()
+    {
+        delete this->broyden;
+        delete[] this->x_out;
+    }
+
+    void solve()
+    {
+        double* x_in = new double[3];
+        this->x_out = new double[3];
+        double* delta_x = new double[3];
+        for (int i = 0; i < 3; ++i)
+            x_in[i] = x_out[i] = 0.;
+
+        auto screen = std::bind(&Broyden_Mixing_Test::Kerker_mock, this, std::placeholders::_1);
+        auto inner_dot = std::bind(&Broyden_Mixing_Test::inner_dot_mock, this, std::placeholders::_1, std::placeholders::_2);
+        this->xdata.resize(this->mixing_ndim, 3, sizeof(double));
+
+        double residual = 10.;
+
+        while (niter < maxiter)
+        {
+            this->iteration(x_in, x_out);
+            niter++;
+
+            for (int i = 0; i < 3; ++i)
+            {
+                delta_x[i] = x_out[i] - x_in[i];
+            }
+            residual = this->inner_dot_mock(delta_x, delta_x);
+            if (residual <= thr)
+            {
+                break;
+            }
+
+            this->broyden->push_data(this->xdata, x_in, x_out, screen, true);
+
+            this->broyden->cal_coef(this->xdata, inner_dot);
+
+            this->broyden->mix_data(this->xdata, x_in);
+        }
+
+        delete[] x_in;
+        delete[] delta_x;
+    }
+
+    void iteration(double* x_in, double* x_out)
+    {
+        x_out[0] = (3. * x_in[1] - 2. * x_in[2] + 20.) / 8.;
+        x_out[1] = (-4. * x_out[0] + 1. * x_in[2] + 33.) / 11.;
+        x_out[2] = (-6. * x_out[0] - 3. * x_out[1] + 36.) / 12.;
+    }
+
+    void Kerker_mock(double *drho){}
+
+    double inner_dot_mock(double* x1, double* x2)
+    {
+        double xnorm = 0.0;
+        for (int ir = 0; ir < 3; ++ir)
+        {
+            xnorm += x1[ir] * x2[ir];
+        }
+        return xnorm;
+    }
+};
+
+
+TEST_F(Broyden_Mixing_Test, Solve_LinearEq)
+{
+    solve();
+    EXPECT_NEAR(x_out[0], 3., DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[1], 2., DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[2], 1., DOUBLETHRESHOLD);
+    ASSERT_EQ(niter, 5);
+}

source/module_base/module_mixing/test/plain_mixing_test.cpp

@@ -0,0 +1,96 @@
+#include "gtest/gtest.h"
+#include "../plain_mixing.h"
+
+#define DOUBLETHRESHOLD 1e-8
+
+class Plain_Mixing_Test : public testing::Test
+{
+protected:
+    const double mixing_beta = 0.8;
+    Base_Mixing::Mixing_Data xdata;
+    Base_Mixing::Plain_Mixing* plain = nullptr;
+
+    double thr = 1e-5;
+    int niter = 0;
+    int maxiter = 50;
+    double* x_out = nullptr;
+
+    void SetUp()
+    {
+        this->plain = new Base_Mixing::Plain_Mixing(this->mixing_beta);
+    }
+    void TearDown()
+    {
+        delete this->plain;
+        delete[] this->x_out;
+    }
+
+    void solve()
+    {
+        double* x_in = new double[3];
+        this->x_out = new double[3];
+        double* delta_x = new double[3];
+        for (int i = 0; i < 3; ++i)
+            x_in[i] = x_out[i] = 0.;
+
+        auto screen = std::bind(&Plain_Mixing_Test::Kerker_mock, this, std::placeholders::_1);
+        auto inner_dot = std::bind(&Plain_Mixing_Test::inner_dot_mock, this, std::placeholders::_1, std::placeholders::_2);
+        this->xdata.resize(1, 3, sizeof(double));
+
+        double residual = 10.;
+
+        while (niter < maxiter)
+        {
+            this->iteration(x_in, x_out);
+            niter++;
+
+            for (int i = 0; i < 3; ++i)
+            {
+                delta_x[i] = x_out[i] - x_in[i];
+            }
+            residual = this->inner_dot_mock(delta_x, delta_x);
+            if (residual <= thr)
+            {
+                break;
+            }
+
+            this->plain->push_data(this->xdata, x_in, x_out, screen, true);
+
+            this->plain->cal_coef(this->xdata, inner_dot);
+
+            this->plain->mix_data(this->xdata, x_in);
+        }
+
+        delete[] x_in;
+        delete[] delta_x;
+    }
+
+    void iteration(double* x_in, double* x_out)
+    {
+        x_out[0] = (3. * x_in[1] - 2. * x_in[2] + 20.) / 8.;
+        x_out[1] = (-4. * x_out[0] + 1. * x_in[2] + 33.) / 11.;
+        x_out[2] = (-6. * x_out[0] - 3. * x_out[1] + 36.) / 12.;
+    }
+
+    void Kerker_mock(double *drho){}
+
+    double inner_dot_mock(double* x1, double* x2)
+    {
+        double xnorm = 0.0;
+        for (int ir = 0; ir < 3; ++ir)
+        {
+            xnorm += x1[ir] * x2[ir];
+        }
+        return xnorm;
+    }
+};
+
+
+TEST_F(Plain_Mixing_Test, Solve_LinearEq)
+{
+    solve();
+    EXPECT_NEAR(x_out[0], 3.0000363451961291616, DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[1], 1.999993538038743468, DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[2], 0.99998344289224938564, DOUBLETHRESHOLD);
+    ASSERT_EQ(niter, 7);
+}

source/module_base/module_mixing/test/pulay_mixing_test.cpp

@@ -0,0 +1,97 @@
+#include "gtest/gtest.h"
+#include "../pulay_mixing.h"
+
+#define DOUBLETHRESHOLD 1e-8
+
+class Pulay_Mixing_Test : public testing::Test
+{
+protected:
+    const double mixing_beta = 0.8;
+    const int mixing_ndim = 7;
+    Base_Mixing::Mixing_Data xdata;
+    Base_Mixing::Pulay_Mixing* pulay = nullptr;
+
+    double thr = 1e-5;
+    int niter = 0;
+    int maxiter = 50;
+    double* x_out = nullptr;
+
+    void SetUp()
+    {
+        this->pulay = new Base_Mixing::Pulay_Mixing(this->mixing_ndim, this->mixing_beta);
+    }
+    void TearDown()
+    {
+        delete this->pulay;
+        delete[] this->x_out;
+    }
+
+    void solve()
+    {
+        double* x_in = new double[3];
+        this->x_out = new double[3];
+        double* delta_x = new double[3];
+        for (int i = 0; i < 3; ++i)
+            x_in[i] = x_out[i] = 0.;
+
+        auto screen = std::bind(&Pulay_Mixing_Test::Kerker_mock, this, std::placeholders::_1);
+        auto inner_dot = std::bind(&Pulay_Mixing_Test::inner_dot_mock, this, std::placeholders::_1, std::placeholders::_2);
+        this->xdata.resize(this->mixing_ndim, 3, sizeof(double));
+
+        double residual = 10.;
+
+        while (niter < maxiter)
+        {
+            this->iteration(x_in, x_out);
+            niter++;
+
+            for (int i = 0; i < 3; ++i)
+            {
+                delta_x[i] = x_out[i] - x_in[i];
+            }
+            residual = this->inner_dot_mock(delta_x, delta_x);
+            if (residual <= thr)
+            {
+                break;
+            }
+
+            this->pulay->push_data(this->xdata, x_in, x_out, screen, true);
+
+            this->pulay->cal_coef(this->xdata, inner_dot);
+
+            this->pulay->mix_data(this->xdata, x_in);
+        }
+
+        delete[] x_in;
+        delete[] delta_x;
+    }
+
+    void iteration(double* x_in, double* x_out)
+    {
+        x_out[0] = (3. * x_in[1] - 2. * x_in[2] + 20.) / 8.;
+        x_out[1] = (-4. * x_out[0] + 1. * x_in[2] + 33.) / 11.;
+        x_out[2] = (-6. * x_out[0] - 3. * x_out[1] + 36.) / 12.;
+    }
+
+    void Kerker_mock(double *drho){}
+
+    double inner_dot_mock(double* x1, double* x2)
+    {
+        double xnorm = 0.0;
+        for (int ir = 0; ir < 3; ++ir)
+        {
+            xnorm += x1[ir] * x2[ir];
+        }
+        return xnorm;
+    }
+};
+
+
+TEST_F(Pulay_Mixing_Test, Solve_LinearEq)
+{
+    solve();
+    EXPECT_NEAR(x_out[0], 3., DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[1], 2., DOUBLETHRESHOLD);
+    EXPECT_NEAR(x_out[2], 1., DOUBLETHRESHOLD);
+    ASSERT_EQ(niter, 5);
+}