- move function read_hs() to diago_elpa_utils.h

- rename H-KPonts.dat to be H-KPoints.dat
- delete some redundant information output

Author: pxlxingliang

source/src_pdiag/test/CMakeLists.txt

@@ -4,7 +4,7 @@ AddTest(
   SOURCES diago_elpa_test.cpp
 )
 
-install(FILES H-KPonts.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+install(FILES H-KPoints.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 install(FILES H-GammaOnly.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 install(FILES S-KPoints.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 install(FILES S-GammaOnly.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})

source/src_pdiag/test/H-KPonts.dat renamed to source/src_pdiag/test/H-KPoints.dat

@@ -7,7 +7,8 @@
 
 #define DETAILINFO       true
 #define PASSTHRESHOLD    1e-4
-#define PRINT_NEW_MATRIX false // print out the H matrix after doing block-cycle distribution
+#define PRINT_NEW_MATRIX true // print out the H matrix after doing block-cycle distribution
+#define PRINT_HS    false
 
 /************************************************
  *  unit test of ELPA
@@ -42,10 +43,10 @@ template <class T> class ElpaPrepare
     int n = 10; // dimension of H/S matrix
     int nblk = 1; // dimension of sub-matrix
     int sparsity = 100; // sparsity of H matrix;
-    T *hmatrix = 0;
-    T *smatrix = 0;
-    double *e_lapack = 0; // eigenvalues solvered by LAPACK
-    double *e_elpa = 0;
+    T *hmatrix = nullptr;
+    T *smatrix = nullptr;
+    double *e_lapack = nullptr; // eigenvalues solved by LAPACK
+    double *e_elpa = nullptr;
     double lapack_time = 0.0, elpa_time = 0.0;
     int mypnum = 0;
 
@@ -54,66 +55,54 @@ template <class T> class ElpaPrepare
         LCAO_DIAGO_TEST::random_hs<T>(hmatrix, smatrix, n, sparsity);
     }
 
-    inline void readmatrix(std::ifstream &inf, T *&matrix, int &npw)
-    {
-        inf >> npw;
-        matrix = new T[npw * npw];
-        for (int i = 0; i < npw; i++)
-        {
-            for (int j = i; j < npw; j++)
-            {
-                inf >> matrix[i * npw + j];
-                if (i != j)
-                {
-                    matrix[j * npw + i] = matrix[i * npw + j];
-                }
-            }
-        }
-    }
-
-    void read_HS()
+    bool read_HS()
     {
+        bool readhfile = false;
+        bool readsfile = false;
         if (mypnum == 0)
         {
-            int hdim = 1;
-            int sdim = 1;
-
-            std::ifstream ifs;
+            std::vector<T> htmp,stmp;
             if (std::is_same<T, double>::value)
             {
-                ifs.open("H-GammaOnly.dat");
+                readhfile = LCAO_DIAGO_TEST::read_hs<std::vector<T>>(std::string("H-GammaOnly.dat"),htmp);
+                readsfile = LCAO_DIAGO_TEST::read_hs<std::vector<T>>(std::string("S-GammaOnly.dat"),stmp);
             }
             else if ((std::is_same<T, std::complex<double>>::value))
-                ifs.open("H-KPonts.dat");
-            readmatrix(ifs, hmatrix, hdim);
-            ifs.close();
-
-            if (std::is_same<T, double>::value)
             {
-                ifs.open("S-GammaOnly.dat");
+                readhfile = LCAO_DIAGO_TEST::read_hs<std::vector<T>>(std::string("H-KPoints.dat"),htmp);
+                readsfile = LCAO_DIAGO_TEST::read_hs<std::vector<T>>(std::string("S-KPoints.dat"),stmp);
             }
-            else if ((std::is_same<T, std::complex<double>>::value))
-                ifs.open("S-KPoints.dat");
-            readmatrix(ifs, smatrix, sdim);
-            ifs.close();
+            if (htmp.size() != stmp.size())
+            {
+                printf("Error: dimensions of H and S are not equal, %d, %d", htmp.size(), stmp.size());
+                readhfile = readsfile = false;
+            }
+
+            hmatrix = new T[htmp.size()];
+            smatrix = new T[stmp.size()];
 
-            if (hdim != sdim)
+            for(int i=0;i<htmp.size();i++)
             {
-                printf("Error: dimensions of H and S are not equal, %d, %d", hdim, sdim);
-                exit(1);
+                hmatrix[i] = htmp[i];
+                smatrix[i] = stmp[i];
             }
-            n = hdim;
+            
+            n = sqrt(stmp.size());
             nblk = 1;
         }
         MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
         MPI_Bcast(&nblk, 1, MPI_INT, 0, MPI_COMM_WORLD);
+        MPI_Bcast(&readhfile, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
+        MPI_Bcast(&readsfile, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
         e_lapack = new double[n];
         e_elpa = new double[n];
         for (int i = 0; i < n; i++)
         {
             e_lapack[i] = 0.0;
             e_elpa[i] = 1.0;
         }
+        if (readhfile && readsfile) return true;
+        else return false;
     }
 
     void new_matrix()
@@ -131,6 +120,7 @@ template <class T> class ElpaPrepare
 
     void print_hs()
     {
+        if(! PRINT_HS) {return;}
         if (mypnum == 0)
         {
             std::ofstream fp("hmatrix.dat");
@@ -207,16 +197,16 @@ class ElpaComplexDoubleTest : public ::testing::TestWithParam<ElpaPrepare<std::c
 TEST(VerifyElpaDiaoDouble, ReadHS)
 {
     ElpaPrepare<double> edp;
-    edp.read_HS();
-    // edp.print_hs();
+    ASSERT_TRUE(edp.read_HS());
+    edp.print_hs();
     edp.run_diago();
 }
 
 TEST(VerifyElpaDiaoComplexDouble, ReadHS)
 {
     ElpaPrepare<std::complex<double>> edp;
-    edp.read_HS();
-    // edp.print_hs();
+    ASSERT_TRUE(edp.read_HS());
+    edp.print_hs();
     edp.run_diago();
 }
 

source/src_pdiag/test/diago_elpa_test.cpp

@@ -136,6 +136,32 @@ template <class T> inline void print_matrix(std::ofstream &fp, T *matrix, int &n
     }
 }
 
+template <class T> bool read_hs(std::string fname, T &matrix)
+{
+    int ndim;
+    std::ifstream inf(fname);
+    if(! inf.is_open())
+    {
+        std::cout << "Error: open file " << fname << " failed, skip!" << std::endl;
+        return false;
+    }
+    inf >> ndim;
+    matrix.resize(ndim*ndim);
+    for (int i = 0; i < ndim; i++)
+    {
+        for (int j = i; j < ndim; j++)
+        {
+            inf >> matrix[i * ndim + j];
+            if (i != j)
+            {
+                matrix[j * ndim + i] = matrix[i * ndim + j];
+            }
+        }
+    }
+    inf.close();
+    return true;
+}
+
 int elpa_sethandle(elpa_t &handle,
                    int nFull,
                    int nev,
@@ -322,7 +348,6 @@ void lapack_diago(double *hmatrix, double *smatrix, double *e, int &nFull)
     }
 
     dsygv_(&itype, &jobz, &uplo, &nFull, a, &nFull, b, &nFull, e, ev, &lwork, &info);
-    std::cout << "lapack solver over: e[0]=" << e[0] << std::endl;
     if (info != 0)
     {
         std::cout << "ERROR: solvered by LAPACK error, info=" << info << std::endl;
@@ -364,12 +389,12 @@ void lapack_diago(std::complex<double> *hmatrix, std::complex<double> *smatrix,
     delete[] rwork;
 }
 
-// produce the random H/S matrix
 double conj(double &a)
 {
     return a;
 }
 
+// produce the random H/S matrix
 template <class T> void random_hs(T *hmatrix, T *smatrix, int &n, int &sparsity)
 {
     int min = 0;