UT and annotations for inverse_matrix and heapsort

source/module_base/mymath3.cpp

@@ -11,11 +11,13 @@ void heapAjust(double *r, int *ind, int s, int m)
     rc = r[s];
     ic = ind[s];
 
-    for (j = 2 * s;j <= m;j *= 2)
+    for (j = 2 * s; j <= m; j *= 2)
     {
-        if (j < m && (r[j] < r[j+1])) j++;
+        if (j < m && (r[j] < r[j + 1]))
+            j++;
 
-        if (!(rc < r[j])) break;
+        if (!(rc < r[j]))
+            break;
 
         r[s] = r[j];
 
@@ -32,24 +34,24 @@ void heapAjust(double *r, int *ind, int s, int m)
 
 void heapsort(const int n, double *r, int *ind)
 {
-    ModuleBase::timer::tick("mymath","heapsort");
+    ModuleBase::timer::tick("mymath", "heapsort");
     int i, ic;
     double rc;
 
     if (ind[0] == 0)
     {
-        for (i = 0;i < n;i++)
+        for (i = 0; i < n; i++)
         {
             ind[i] = i;
         }
     }
 
-    for (i = n / 2;i >= 0;i--)
+    for (i = n / 2; i >= 0; i--)
     {
         heapAjust(r, ind, i, n - 1);
     }
 
-    for (i= n - 1;i > 0;i--)
+    for (i = n - 1; i > 0; i--)
     {
         rc = r[0];
         r[0] = r[i];
@@ -59,7 +61,7 @@ void heapsort(const int n, double *r, int *ind)
         ind[i] = ic;
         heapAjust(r, ind, 0, i - 1);
     }
-    ModuleBase::timer::tick("mymath","heapsort");
+    ModuleBase::timer::tick("mymath", "heapsort");
     return;
 }
 
@@ -89,82 +91,83 @@ void hpsort(int n, double *ra, int *ind)
     int i, ir, j, k, iind;
     double rra;
 
-    if (ind[1] == 0)
+    if (ind[0] == 0)
     {
-        for (i = 1;i <= n;i++)
-            ind[i] = i;
+        for (i = 1; i <= n; i++)
+            ind[i - 1] = i;
     }
 
-    if (n < 2) return;  // nothing to order
+    if (n < 2)
+        return; // nothing to order
 
-    k  = n / 2 + 1;
+    k = n / 2;
 
-    ir = n;
+    ir = n - 1;
 
     while (true)
     {
-        if (k > 1)      // still in hiring phase
+        if (k > 0) // still in hiring phase
         {
-            k   = k - 1;
+            k = k - 1;
             rra = ra[k];
             iind = ind[k];
         }
-        else                 // in retirement-promotion phase.
+        else // in retirement-promotion phase.
         {
-            rra = ra[ir];      // clear a space at the end of the array
-            iind = ind[ir];    //
-            ra[ir] = ra[1];    // retire the top of the heap into it
-            ind[ir] = ind[1];  //
-            ir     = ir - 1;   // decrease the size of the corporation
+            rra = ra[ir]; // clear a space at the end of the array
+            iind = ind[ir]; //
+            ra[ir] = ra[0]; // retire the top of the heap into it
+            ind[ir] = ind[0]; //
+            ir = ir - 1; // decrease the size of the corporation
 
-            if (ir == 1)   // done with the last promotion
+            if (ir == 0) // done with the last promotion
             {
-                ra[1] = rra;    // the least competent worker at all //
-                ind[1] = iind;  //
+                ra[0] = rra; // the least competent worker at all //
+                ind[0] = iind; //
                 return;
             }
         }
 
-        i = k;               // wheter in hiring or promotion phase, we
+        i = k; // wheter in hiring or promotion phase, we
 
-        j = k + k;           // set up to place rra in its proper level
+        j = k + k + 1; // set up to place rra in its proper level
 
         while (j <= ir)
         {
             if (j < ir)
             {
-                if (ra[j] < ra[j+1])   // compare to better underling
+                if (ra[j] < ra[j + 1]) // compare to better underling
                 {
                     j = j + 1;
                 }
-                else if (ra[j] == ra[j+1])
+                else if (ra[j] == ra[j + 1])
                 {
-                    if (ind[j] < ind[j+1])
+                    if (ind[j] < ind[j + 1])
                         j = j + 1;
                 }
             }
 
-            if (rra < ra[j])   // demote rra
+            if (rra < ra[j]) // demote rra
             {
                 ra[i] = ra[j];
                 ind[i] = ind[j];
                 i = j;
-                j = j + j;
+                j = j + j + 1;
             }
             else if (rra == ra[j])
             {
-                if (iind < ind[j])   // demote rra
+                if (iind < ind[j]) // demote rra
                 {
                     ra[i] = ra[j];
                     ind[i] = ind[j];
                     i = j;
-                    j = j + j;
+                    j = j + j + 1;
                 }
                 else
-                    j = ir + 1;         // set j to terminate do-while loop
+                    j = ir + 1; // set j to terminate do-while loop
             }
-            else                   // this is the right place for rra
-                j = ir + 1;           // set j to terminate do-while loop
+            else // this is the right place for rra
+                j = ir + 1; // set j to terminate do-while loop
         }
 
         ra[i] = rra;
@@ -173,4 +176,4 @@ void hpsort(int n, double *ra, int *ind)
     }
 }
 
-}
+} // namespace ModuleBase

source/module_base/test/CMakeLists.txt

@@ -86,3 +86,12 @@ AddTest(
   TARGET base_math_bspline
   SOURCES math_bspline_test.cpp ../math_bspline.cpp
 )
+AddTest(
+  TARGET base_inverse_matrix
+  LIBS ${math_libs}
+  SOURCES inverse_matrix_test.cpp ../inverse_matrix.cpp ../complexmatrix.cpp ../matrix.cpp ../timer.cpp
+)
+AddTest(
+  TARGET base_mymath
+  SOURCES mymath_test.cpp ../mymath3.cpp ../timer.cpp
+)

source/module_base/test/inverse_matrix_test.cpp

@@ -0,0 +1,61 @@
+#include "../inverse_matrix.h"
+#include "gtest/gtest.h"
+#include <iostream>
+
+/************************************************
+ *  unit test of class Inverse_Matrix_Complex
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - Inverse
+ *     - use Inverse_Matrix_Complex to inverse a Hermite matrix
+ *     - functions: init and using_zheev
+ *
+ */
+
+//a mock function of WARNING_QUIT, to avoid the uncorrected call by matrix.cpp at line 37.
+namespace ModuleBase
+{
+    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+}
+
+TEST(InverseMatrixComplexTest,Inverse)
+{
+	int dim = 10;
+	ModuleBase::ComplexMatrix B(dim,dim);
+	ModuleBase::ComplexMatrix C(dim,dim);
+	ModuleBase::ComplexMatrix D(dim,dim);
+	double a;
+	double b;
+	double c;
+	// construct a Hermite matrix
+	for(int j=0;j<dim;j++)
+	{
+		for(int i=0;i<=j;i++)
+		{
+			if (i==j)
+			{
+				c = std::rand();
+				B(i,j) = std::complex<double> (c,0.0);
+			}
+			else
+			{
+				a = std::rand();
+				b = std::rand();
+				B(i,j) = std::complex<double> (a,b);
+				B(j,i) = conj(B(i,j));
+			}
+		}
+	}
+	ModuleBase::Inverse_Matrix_Complex IMC;
+	IMC.init(dim);
+	IMC.using_zheev(B,C);
+	D = B * C;
+	for(int i=0;i<dim;i++)
+	{
+		EXPECT_NEAR(D(i,i).real(),1.0,1e-14);
+		EXPECT_NEAR(D(i,i).imag(),0.0,1e-14);
+		//std::cout << D(i,i).real() << " " << D(i,i).imag()  << std::endl;
+	}
+}

source/module_base/test/mymath_test.cpp

@@ -0,0 +1,66 @@
+#include "../mymath.h"
+#include "gtest/gtest.h"
+#include <iostream>
+#include <algorithm>
+
+/************************************************
+ *  unit test of heap sort functions
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - Heapsort
+ *     - heapsort function in mymath3.cpp
+ *   - Hpsort
+ *     - hpsort function in mymath3.cpp
+ *
+ */
+
+class MymathTest : public testing::Test
+{
+protected:
+	int number;
+};
+
+TEST_F(MymathTest,Heapsort)
+{
+	number = 5;
+	double rr[number];
+	int index[number];
+	rr[0] = 10;
+	rr[1] = 9;
+	rr[2] = 8;
+	rr[3] = 7;
+	rr[4] = 6;
+	index[0] = 0;
+	//for (int i=0;i<number;i++)
+	//	std::cout << i << " " << rr[i] << std::endl;
+	ModuleBase::heapsort(number,rr,index);
+	for (int i=0;i<number-1;i++)
+	{
+		//std::cout << i << " " << rr[i] << std::endl;
+		EXPECT_LT(rr[i],rr[i+1]);
+	}
+}
+
+TEST_F(MymathTest,Hpsort)
+{
+	number = 5;
+	double rr[number];
+	int index[number];
+	rr[0] = 10;
+	rr[1] = 9;
+	rr[2] = 8;
+	rr[3] = 7;
+	rr[4] = 6;
+	index[0] = 0;
+	//for (int i=0;i<number;i++)
+	//	std::cout << i << " " << rr[i] << std::endl;
+	ModuleBase::hpsort(number,rr,index);
+	//std::sort(rr,rr+number);
+	for (int i=0;i<number-1;i++)
+	{
+		//std::cout << i << " " << rr[i] << std::endl;
+		EXPECT_LT(rr[i],rr[i+1]);
+	}
+}