Add the USE_KML option and change abstol/orfac for KML

Author: ZhouXY-PKU

CMakeLists.txt

@@ -17,6 +17,7 @@ option(USE_CUDA_MPI "Enable CUDA-aware MPI" OFF)
 option(USE_CUDA_ON_DCU "Enable CUDA on DCU" OFF)
 option(USE_ROCM "Enable ROCm" OFF)
 option(USE_DSP "Enable DSP" OFF)
+option(USE_KML "Enable Kunpeng Math Library" OFF)
 option(USE_SW "Enable SW Architecture" OFF)
 
 option(USE_ABACUS_LIBM "Build libmath from source to speed up" OFF)
@@ -337,6 +338,50 @@ if (USE_DSP)
   target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_device.a)
   target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_host.a)
 endif()
+
+
+if(USE_KML)
+    add_compile_definitions(__KML)
+    message(STATUS "Huawei KML support enabled. Defining __KML.")
+#    TODO: Create FindKML.cmake
+#    if(NOT DEFINED KML_ROOT)
+#        if(DEFINED ENV{KML_ROOT})
+#            set(KML_ROOT $ENV{KML_ROOT})
+#        else()
+#            message(WARNING "KML_ROOT is not set. Trying default system paths for KML.")
+#        endif()
+#    endif()
+#
+#    find_library(KML_BLAS_LIB NAMES kblas PATHS ${KML_ROOT}/lib ${KML_ROOT}/lib64 NO_DEFAULT_PATH)
+#    find_library(KML_LAPACK_LIB NAMES klapack_full PATHS ${KML_ROOT}/lib ${KML_ROOT}/lib64 NO_DEFAULT_PATH)
+#    find_library(KML_SCALAPACK_LIB NAMES kscalapack_full PATHS ${KML_ROOT}/lib ${KML_ROOT}/lib64 NO_DEFAULT_PATH)
+#    find_library(KML_FFTW_LIB NAMES fftw3 PATHS ${KML_ROOT}/lib ${KML_ROOT}/lib64 NO_DEFAULT_PATH)
+#
+#    set(KML_LIBS_FOUND TRUE)
+#    foreach(LIB_VAR KML_BLAS_LIB KML_LAPACK_LIB KML_SCALAPACK_LIB KML_FFTW_LIB)
+#        if(NOT ${LIB_VAR})
+#            message(WARNING "${LIB_VAR} not found in KML_ROOT! Please check your KML installation.")
+#            set(KML_LIBS_FOUND FALSE)
+#        endif()
+#    endforeach()
+#
+#    if(KML_LIBS_FOUND)
+#        target_link_libraries(abacus PUBLIC
+#            ${KML_BLAS_LIB}
+#            ${KML_LAPACK_LIB}
+#            ${KML_SCALAPACK_LIB}
+#            ${KML_FFTW_LIB}
+#        )
+#        message(STATUS "Huawei KML libraries found and linked successfully.")
+#    else()
+#        message(FATAL_ERROR "Failed to find all required KML libraries. Aborting.")
+#    endif()
+#
+#    set(BLAS_libraries ${KML_BLAS_LIB})
+#    set(LAPACK_libraries ${KML_LAPACK_LIB})
+endif(USE_KML)
+
+
 if (USE_SW)
   add_compile_definitions(__SW)
   set(SW ON)

source/source_base/module_external/lapack_connector.h

@@ -30,6 +30,23 @@
 #include "../complexmatrix.h"
 #include "../global_function.h"
 
+#include <limits>
+// =========================================================
+// Tolerances for LAPACK/ScaLAPACK eigenvalue routines
+// =========================================================
+// Huawei KML strictly validates input parameters.
+// It rejects abstol=0 and orfac=-1, which are standard
+// defaults in open-source LAPACK to trigger internal logic.
+// We must explicitly pass the mathematically equivalent defaults.
+#ifdef __KML
+    constexpr double LAPACK_ABSTOL = 2*std::numeric_limits<double>::min();          // 2*PDLAMCH('S')
+    constexpr double LAPACK_ORFAC  = 2.0e-12;                                       // Default value
+#else
+    constexpr double LAPACK_ABSTOL = 0.0;
+    constexpr double LAPACK_ORFAC  = -1.0;
+#endif
+// =========================================================
+
 //Naming convention of lapack subroutines : ammxxx, where
 //"a" specifies the data type:
 // - s stands for float

source/source_base/module_external/scalapack_connector.h

@@ -1,6 +1,23 @@
 #ifndef SCALAPACK_CONNECTOR_H
 #define SCALAPACK_CONNECTOR_H
 
+#include <limits>
+// =========================================================
+// Tolerances for LAPACK/ScaLAPACK eigenvalue routines
+// =========================================================
+// Huawei KML strictly validates input parameters.
+// It rejects abstol=0 and orfac=-1, which are standard 
+// defaults in open-source ScaLAPACK to trigger internal logic.
+// We must explicitly pass the mathematically equivalent defaults.
+#ifdef __KML
+    constexpr double SCALAPACK_ABSTOL = 2*std::numeric_limits<double>::min();          // 2*PDLAMCH('S')
+    constexpr double SCALAPACK_ORFAC  = 2.0e-12;                                       // Default value
+#else
+    constexpr double SCALAPACK_ABSTOL = 0.0;
+    constexpr double SCALAPACK_ORFAC  = -1.0;
+#endif
+// =========================================================
+
 #ifdef __MPI
 
 #include <complex>

source/source_hsolver/diago_lapack.cpp

@@ -98,7 +98,7 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0;
     double vl = 0, vu = 0;
-    const double abstol = 0, orfac = -1;
+    const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC;
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(PARAM.globalv.nlocal, 0);
@@ -206,7 +206,7 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     const char jobz = 'V', range = 'I', uplo = 'U';
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0;
-    const double abstol = 0, orfac = -1;
+    const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC;
 
     const double vl = 0, vu = 0;
     std::vector<std::complex<double>> work(1, 0);

source/source_hsolver/diago_pxxxgvx.cpp

@@ -502,8 +502,8 @@ void pxxxgvx_diag(const int* const desc,
         int lrwork = -1;
         int liwork = -1;
         int info = 0;
-        const typename GetTypeReal<T>::type abstol = 0;
-        const typename GetTypeReal<T>::type orfac = -1;
+        const typename GetTypeReal<T>::type abstol = SCALAPACK_ABSTOL;
+        const typename GetTypeReal<T>::type orfac = SCALAPACK_ORFAC;
         const typename GetTypeReal<T>::type vl = 0;
         const typename GetTypeReal<T>::type vu = 0;
         std::vector<T> work(1, 0);
@@ -667,4 +667,4 @@ template void pxxxgvx_diag(const int* const desc,
 
 #endif
 
-} // namespace hsolver
+} // namespace hsolver

source/source_hsolver/diago_scalapack.cpp

@@ -95,7 +95,7 @@ namespace hsolver
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0;
     double vl = 0, vu = 0;
-    const double abstol = 0, orfac = -1;
+    const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC;
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(PARAM.globalv.nlocal, 0);
@@ -219,7 +219,7 @@ namespace hsolver
     const char jobz = 'V', range = 'I', uplo = 'U';
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0;
-    const double abstol = 0, orfac = -1;
+    const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC;
     //Note: pzhegvx_ has a bug
     //      We must give vl,vu a value, although we do not use range 'V'
     //      We must give rwork at least a memory of sizeof(double) * 3
@@ -445,4 +445,4 @@ namespace hsolver
     }
 }
 
-} // namespace hsolver
+} // namespace hsolver