Merge pull request #758 from hongriTianqi/develop

hongriTianqi · web-flow · commit ebf1afe35af9 · 2022-03-02T09:35:35.000+08:00
UT and annotations for inverse_matrix and heapsort
diff --git a/source/module_base/mymath3.cpp b/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
diff --git a/source/module_base/test/CMakeLists.txt b/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
+)
diff --git a/source/module_base/test/inverse_matrix_test.cpp b/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;
+	}
+}
diff --git a/source/module_base/test/mymath_test.cpp b/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]);
+	}
+}

Original file line number	Diff line number	Diff line change
`@@ -11,11 +11,13 @@ void heapAjust(double r, int ind, int s, int m)`
`11`	`11`	`rc = r[s];`
`12`	`12`	`ic = ind[s];`
`13`	`13`
`14`		`- for (j = 2 * s;j <= m;j *= 2)`
	`14`	`+ for (j = 2 * s; j <= m; j *= 2)`
`15`	`15`	`{`
`16`		`- if (j < m && (r[j] < r[j+1])) j++;`
	`16`	`+ if (j < m && (r[j] < r[j + 1]))`
	`17`	`+ j++;`
`17`	`18`
`18`		`- if (!(rc < r[j])) break;`
	`19`	`+ if (!(rc < r[j]))`
	`20`	`+ break;`
`19`	`21`
`20`	`22`	`r[s] = r[j];`
`21`	`23`
`@@ -32,24 +34,24 @@ void heapAjust(double r, int ind, int s, int m)`
`32`	`34`
`33`	`35`	`void heapsort(const int n, double r, int ind)`
`34`	`36`	`{`
`35`		`- ModuleBase::timer::tick("mymath","heapsort");`
	`37`	`+ ModuleBase::timer::tick("mymath", "heapsort");`
`36`	`38`	`int i, ic;`
`37`	`39`	`double rc;`
`38`	`40`
`39`	`41`	`if (ind[0] == 0)`
`40`	`42`	`{`
`41`		`- for (i = 0;i < n;i++)`
	`43`	`+ for (i = 0; i < n; i++)`
`42`	`44`	`{`
`43`	`45`	`ind[i] = i;`
`44`	`46`	`}`
`45`	`47`	`}`
`46`	`48`
`47`		`- for (i = n / 2;i >= 0;i--)`
	`49`	`+ for (i = n / 2; i >= 0; i--)`
`48`	`50`	`{`
`49`	`51`	`heapAjust(r, ind, i, n - 1);`
`50`	`52`	`}`
`51`	`53`
`52`		`- for (i= n - 1;i > 0;i--)`
	`54`	`+ for (i = n - 1; i > 0; i--)`
`53`	`55`	`{`
`54`	`56`	`rc = r[0];`
`55`	`57`	`r[0] = r[i];`
`@@ -59,7 +61,7 @@ void heapsort(const int n, double r, int ind)`
`59`	`61`	`ind[i] = ic;`
`60`	`62`	`heapAjust(r, ind, 0, i - 1);`
`61`	`63`	`}`
`62`		`- ModuleBase::timer::tick("mymath","heapsort");`
	`64`	`+ ModuleBase::timer::tick("mymath", "heapsort");`
`63`	`65`	`return;`
`64`	`66`	`}`
`65`	`67`
`@@ -89,82 +91,83 @@ void hpsort(int n, double ra, int ind)`
`89`	`91`	`int i, ir, j, k, iind;`
`90`	`92`	`double rra;`
`91`	`93`
`92`		`- if (ind[1] == 0)`
	`94`	`+ if (ind[0] == 0)`
`93`	`95`	`{`
`94`		`- for (i = 1;i <= n;i++)`
`95`		`- ind[i] = i;`
	`96`	`+ for (i = 1; i <= n; i++)`
	`97`	`+ ind[i - 1] = i;`
`96`	`98`	`}`
`97`	`99`
`98`		`- if (n < 2) return; // nothing to order`
	`100`	`+ if (n < 2)`
	`101`	`+ return; // nothing to order`
`99`	`102`
`100`		`- k = n / 2 + 1;`
	`103`	`+ k = n / 2;`
`101`	`104`
`102`		`- ir = n;`
	`105`	`+ ir = n - 1;`
`103`	`106`
`104`	`107`	`while (true)`
`105`	`108`	`{`
`106`		`- if (k > 1) // still in hiring phase`
	`109`	`+ if (k > 0) // still in hiring phase`
`107`	`110`	`{`
`108`		`- k = k - 1;`
	`111`	`+ k = k - 1;`
`109`	`112`	`rra = ra[k];`
`110`	`113`	`iind = ind[k];`
`111`	`114`	`}`
`112`		`- else // in retirement-promotion phase.`
	`115`	`+ else // in retirement-promotion phase.`
`113`	`116`	`{`
`114`		`- rra = ra[ir]; // clear a space at the end of the array`
`115`		`- iind = ind[ir]; //`
`116`		`- ra[ir] = ra[1]; // retire the top of the heap into it`
`117`		`- ind[ir] = ind[1]; //`
`118`		`- ir = ir - 1; // decrease the size of the corporation`
	`117`	`+ rra = ra[ir]; // clear a space at the end of the array`
	`118`	`+ iind = ind[ir]; //`
	`119`	`+ ra[ir] = ra[0]; // retire the top of the heap into it`
	`120`	`+ ind[ir] = ind[0]; //`
	`121`	`+ ir = ir - 1; // decrease the size of the corporation`
`119`	`122`
`120`		`- if (ir == 1) // done with the last promotion`
	`123`	`+ if (ir == 0) // done with the last promotion`
`121`	`124`	`{`
`122`		`- ra[1] = rra; // the least competent worker at all //`
`123`		`- ind[1] = iind; //`
	`125`	`+ ra[0] = rra; // the least competent worker at all //`
	`126`	`+ ind[0] = iind; //`
`124`	`127`	`return;`
`125`	`128`	`}`
`126`	`129`	`}`
`127`	`130`
`128`		`- i = k; // wheter in hiring or promotion phase, we`
	`131`	`+ i = k; // wheter in hiring or promotion phase, we`
`129`	`132`
`130`		`- j = k + k; // set up to place rra in its proper level`
	`133`	`+ j = k + k + 1; // set up to place rra in its proper level`
`131`	`134`
`132`	`135`	`while (j <= ir)`
`133`	`136`	`{`
`134`	`137`	`if (j < ir)`
`135`	`138`	`{`
`136`		`- if (ra[j] < ra[j+1]) // compare to better underling`
	`139`	`+ if (ra[j] < ra[j + 1]) // compare to better underling`
`137`	`140`	`{`
`138`	`141`	`j = j + 1;`
`139`	`142`	`}`
`140`		`- else if (ra[j] == ra[j+1])`
	`143`	`+ else if (ra[j] == ra[j + 1])`
`141`	`144`	`{`
`142`		`- if (ind[j] < ind[j+1])`
	`145`	`+ if (ind[j] < ind[j + 1])`
`143`	`146`	`j = j + 1;`
`144`	`147`	`}`
`145`	`148`	`}`
`146`	`149`
`147`		`- if (rra < ra[j]) // demote rra`
	`150`	`+ if (rra < ra[j]) // demote rra`
`148`	`151`	`{`
`149`	`152`	`ra[i] = ra[j];`
`150`	`153`	`ind[i] = ind[j];`
`151`	`154`	`i = j;`
`152`		`- j = j + j;`
	`155`	`+ j = j + j + 1;`
`153`	`156`	`}`
`154`	`157`	`else if (rra == ra[j])`
`155`	`158`	`{`
`156`		`- if (iind < ind[j]) // demote rra`
	`159`	`+ if (iind < ind[j]) // demote rra`
`157`	`160`	`{`
`158`	`161`	`ra[i] = ra[j];`
`159`	`162`	`ind[i] = ind[j];`
`160`	`163`	`i = j;`
`161`		`- j = j + j;`
	`164`	`+ j = j + j + 1;`
`162`	`165`	`}`
`163`	`166`	`else`
`164`		`- j = ir + 1; // set j to terminate do-while loop`
	`167`	`+ j = ir + 1; // set j to terminate do-while loop`
`165`	`168`	`}`
`166`		`- else // this is the right place for rra`
`167`		`- j = ir + 1; // set j to terminate do-while loop`
	`169`	`+ else // this is the right place for rra`
	`170`	`+ j = ir + 1; // set j to terminate do-while loop`
`168`	`171`	`}`
`169`	`172`
`170`	`173`	`ra[i] = rra;`
`@@ -173,4 +176,4 @@ void hpsort(int n, double ra, int ind)`
`173`	`176`	`}`
`174`	`177`	`}`
`175`	`178`
`176`		`-}`
	`179`	`+} // namespace ModuleBase`