Merge pull request #1172 from caic99/math_info

feature: log math lib info

Author: dyzheng

CMakeLists.txt

@@ -20,6 +20,7 @@ option(USE_OPENMP " Enable OpenMP in abacus." ON)
 option(ENABLE_ASAN "Enable AddressSanitizer" OFF)
 option(BUILD_TESTING "Build ABACUS unit tests" OFF)
 option(GENERATE_TEST_REPORTS "Enable test report generation" OFF)
+option(INFO "Enable gathering of math library information" OFF)
 
 set(ABACUS_BIN_NAME abacus)
 list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/modules)
@@ -269,6 +270,12 @@ add_compile_definitions(
   TEST_EXX_RADIAL=1
 )
 
+if(INFO)
+  message(STATUS "Will gather math lib info.")
+  add_compile_definitions(GATHER_INFO)
+  # modifications on blas_connector and lapack_connector
+endif()
+
 IF (BUILD_TESTING)
   set(CMAKE_CXX_STANDARD 14) # Required in orbital
   include(CTest)

source/module_base/CMakeLists.txt

@@ -7,6 +7,7 @@ add_library(
     export.cpp
     integral.cpp
     inverse_matrix.cpp
+    gather_math_lib_info.cpp
     global_file.cpp
     global_function.cpp
     global_function_ddotreal.cpp
@@ -16,7 +17,7 @@ add_library(
     math_polyint.cpp
     math_sphbes.cpp
     math_ylmreal.cpp
-	math_bspline.cpp
+    math_bspline.cpp
     math_chebyshev.cpp
     mathzone.cpp
     mathzone_add1.cpp

source/module_base/blas_connector.h

@@ -18,13 +18,13 @@ extern "C"
 	void daxpy_(const int *N, const double *alpha, const double *X, const int *incX, double *Y, const int *incY);
 	void caxpy_(const int *N, const std::complex<float> *alpha, const std::complex<float> *X, const int *incX, std::complex<float> *Y, const int *incY);
 	void zaxpy_(const int *N, const std::complex<double> *alpha, const std::complex<double> *X, const int *incX, std::complex<double> *Y, const int *incY);
-	
+
 	void dcopy_(long const *n, const double *a, int const *incx, double *b, int const *incy);
-	void zcopy_(long const *n, const std::complex<double> *a, int const *incx, std::complex<double> *b, int const *incy); 
+	void zcopy_(long const *n, const std::complex<double> *a, int const *incx, std::complex<double> *b, int const *incy);
 
 	//reason for passing results as argument instead of returning it:
 	//see https://www.numbercrunch.de/blog/2014/07/lost-in-translation/
-	void zdotc_(std::complex<double> *result, const int *n, const std::complex<double> *zx, 
+	void zdotc_(std::complex<double> *result, const int *n, const std::complex<double> *zx,
 		const int *incx, const std::complex<double> *zy, const int *incy);
 	// Peize Lin add ?dot 2017-10-27, to compute d=x*y
 	float sdot_(const int *N, const float *X, const int *incX, const float *Y, const int *incY);
@@ -36,36 +36,36 @@ extern "C"
 	double dznrm2_( const int *n, const std::complex<double> *X, const int *incX );
 
 	// level 2: matrix-std::vector operations, O(n^2) data and O(n^2) work.
-	void dgemv_(const char *transa, const int *m, const int *n, const double *alpha,  const double *a,  
+	void dgemv_(const char *transa, const int *m, const int *n, const double *alpha,  const double *a,
 		const int *lda, const double *x, const int *incx, const double *beta, double *y, const int *incy);
-		
+
 	void zgemv_(const char *trans, const int *m, const int *n, const std::complex<double> *alpha,
 			const std::complex<double> *a, const int *lda, const std::complex<double> *x, const int *incx,
 			const std::complex<double> *beta, std::complex<double> *y, const int *incy);
 
-	void dsymv_(const char *uplo, const int *n, 
-		const double *alpha, const double *a, const int *lda, 
-		const double *x, const int *incx, 
-		const double *beta, double *y, const int *incy);			
+	void dsymv_(const char *uplo, const int *n,
+		const double *alpha, const double *a, const int *lda,
+		const double *x, const int *incx,
+		const double *beta, double *y, const int *incy);
 
     // A := alpha x * y.T + A
 	void dger_(int *m, int *n, double *alpha, double *x, int *incx, double *y, int *incy, double *a, int *lda);
 	void zgerc_(int *m, int *n, std::complex<double> *alpha,std::complex<double> *x, int *incx, std::complex<double> *y, int *incy,std::complex<double> *a, int *lda);
 
 	// level 3: matrix-matrix operations, O(n^2) data and O(n^3) work.
-	
+
 	// Peize Lin add ?gemm 2017-10-27, to compute C = a * A.? * B.? + b * C
 	// A is general
 	void sgemm_(const char *transa, const char *transb, const int *m, const int *n, const int *k,
-		const float *alpha, const float *a, const int *lda, const float *b, const int *ldb, 
+		const float *alpha, const float *a, const int *lda, const float *b, const int *ldb,
 		const float *beta, float *c, const int *ldc);
 	void dgemm_(const char *transa, const char *transb, const int *m, const int *n, const int *k,
-		const double *alpha, const double *a, const int *lda, const double *b, const int *ldb, 
+		const double *alpha, const double *a, const int *lda, const double *b, const int *ldb,
 		const double *beta, double *c, const int *ldc);
 	void zgemm_(const char *transa, const char *transb, const int *m, const int *n, const int *k,
-		const std::complex<double> *alpha, const std::complex<double> *a, const int *lda, const std::complex<double> *b, const int *ldb, 
+		const std::complex<double> *alpha, const std::complex<double> *a, const int *lda, const std::complex<double> *b, const int *ldb,
 		const std::complex<double> *beta, std::complex<double> *c, const int *ldc);
-	
+
 	//a is symmetric
 	void dsymm_(const char *side, const char *uplo, const int *m, const int *n,
 		const double *alpha, const double *a, const int *lda, const double *b, const int *ldb,
@@ -91,50 +91,50 @@ class BlasConnector
 
 	// Peize Lin add 2016-08-04
 	// y=a*x+y
-	static inline 
+	static inline
 	void axpy( const int n, const float alpha, const float *X, const int incX, float *Y, const int incY)
 	{
 		saxpy_(&n, &alpha, X, &incX, Y, &incY);
-	}	
-	static inline 
+	}
+	static inline
 	void axpy( const int n, const double alpha, const double *X, const int incX, double *Y, const int incY)
 	{
 		daxpy_(&n, &alpha, X, &incX, Y, &incY);
-	}	
-	static inline 
+	}
+	static inline
 	void axpy( const int n, const std::complex<float> alpha, const std::complex<float> *X, const int incX, std::complex<float> *Y, const int incY)
 	{
 		caxpy_(&n, &alpha, X, &incX, Y, &incY);
-	}	
-	static inline 
+	}
+	static inline
 	void axpy( const int n, const std::complex<double> alpha, const std::complex<double> *X, const int incX, std::complex<double> *Y, const int incY)
 	{
 		zaxpy_(&n, &alpha, X, &incX, Y, &incY);
-	}	
-	
+	}
+
 	// Peize Lin add 2016-08-04
 	// x=a*x
-	static inline 
+	static inline
 	void scal( const int n,  const float alpha, float *X, const int incX)
 	{
 		sscal_(&n, &alpha, X, &incX);
-	}	
-	static inline 
+	}
+	static inline
 	void scal( const int n, const double alpha, double *X, const int incX)
 	{
 		dscal_(&n, &alpha, X, &incX);
-	}	
-	static inline 
+	}
+	static inline
 	void scal( const int n, const std::complex<float> alpha, std::complex<float> *X, const int incX)
 	{
 		cscal_(&n, &alpha, X, &incX);
-	}	
-	static inline 
+	}
+	static inline
 	void scal( const int n, const std::complex<double> alpha, std::complex<double> *X, const int incX)
 	{
 		zscal_(&n, &alpha, X, &incX);
-	}	
-	
+	}
+
 	// Peize Lin add 2017-10-27
 	// d=x*y
 	static inline
@@ -149,32 +149,32 @@ class BlasConnector
 	}
 
 	// Peize Lin add 2017-10-27, fix bug trans 2019-01-17
-	// C = a * A.? * B.? + b * C 
+	// C = a * A.? * B.? + b * C
 	static inline
 	void gemm(const char transa, const char transb, const int m, const int n, const int k,
-		const float alpha, const float *a, const int lda, const float *b, const int ldb, 
+		const float alpha, const float *a, const int lda, const float *b, const int ldb,
 		const float beta, float *c, const int ldc)
 	{
 		sgemm_(&transb, &transa, &n, &m, &k,
-			&alpha, b, &ldb, a, &lda, 
+			&alpha, b, &ldb, a, &lda,
 			&beta, c, &ldc);
 	}
 	static inline
 	void gemm(const char transa, const char transb, const int m, const int n, const int k,
-		const double alpha, const double *a, const int lda, const double *b, const int ldb, 
+		const double alpha, const double *a, const int lda, const double *b, const int ldb,
 		const double beta, double *c, const int ldc)
 	{
 		dgemm_(&transb, &transa, &n, &m, &k,
-			&alpha, b, &ldb, a, &lda, 
+			&alpha, b, &ldb, a, &lda,
 			&beta, c, &ldc);
 	}
 	static inline
 	void gemm(const char transa, const char transb, const int m, const int n, const int k,
-		const std::complex<double> alpha, const std::complex<double> *a, const int lda, const std::complex<double> *b, const int ldb, 
+		const std::complex<double> alpha, const std::complex<double> *a, const int lda, const std::complex<double> *b, const int ldb,
 		const std::complex<double> beta, std::complex<double> *c, const int ldc)
 	{
 		zgemm_(&transb, &transa, &n, &m, &k,
-			&alpha, b, &ldb, a, &lda, 
+			&alpha, b, &ldb, a, &lda,
 			&beta, c, &ldc);
 	}
 
@@ -196,7 +196,7 @@ class BlasConnector
 		return dznrm2_( &n, X, &incX );
 	}
 
-	// copies a into b	
+	// copies a into b
 	static inline
 	void copy(const long n, const double *a, const int incx, double *b, const int incy)
 	{
@@ -206,8 +206,52 @@ class BlasConnector
 	void copy(const long n, const std::complex<double> *a, const int incx, std::complex<double> *b, const int incy)
 	{
 		zcopy_(&n, a, &incx, b, &incy);
-	}	
-
+	}
 };
 
-#endif
+// If GATHER_INFO is defined, the original function is replaced with a "i" suffix,
+// preventing changes on the original code.
+// The real function call is at gather_math_lib_info.cpp
+#ifdef GATHER_INFO
+
+#define zgemm_ zgemm_i
+void zgemm_i(const char *transa,
+             const char *transb,
+             const int *m,
+             const int *n,
+             const int *k,
+             const std::complex<double> *alpha,
+             const std::complex<double> *a,
+             const int *lda,
+             const std::complex<double> *b,
+             const int *ldb,
+             const std::complex<double> *beta,
+             std::complex<double> *c,
+             const int *ldc);
+
+#define zaxpy_  zaxpy_i
+void zaxpy_i(const int *N,
+            const std::complex<double> *alpha,
+            const std::complex<double> *X,
+            const int *incX,
+            std::complex<double> *Y,
+            const int *incY);
+
+/*
+#define zgemv_ zgemv_i
+
+void zgemv_i(const char *trans,
+             const int *m,
+             const int *n,
+             const std::complex<double> *alpha,
+             const std::complex<double> *a,
+             const int *lda,
+             const std::complex<double> *x,
+             const int *incx,
+             const std::complex<double> *beta,
+             std::complex<double> *y,
+             const int *incy);
+*/
+
+#endif // GATHER_INFO
+#endif // BLAS_CONNECTOR_H

source/module_base/gather_math_lib_info.cpp

@@ -0,0 +1,78 @@
+// This file defines the math lib wrapper for output information before executing computations.
+
+/*
+When INFO is defined in cmake configure, a macro of GATHER_INFO will be defined.
+This macro will be used to output information before executing computations.
+Results will output to OUT/math_info.log, see ModuleBase::Global_File::make_dir_out .
+*/
+
+#include <ostream>
+#undef GATHER_INFO
+#include "module_base/blas_connector.h"
+#include "module_base/lapack_connector.h"
+
+#include <iostream>
+
+void zgemm_i(const char *transa,
+             const char *transb,
+             const int *m,
+             const int *n,
+             const int *k,
+             const std::complex<double> *alpha,
+             const std::complex<double> *a,
+             const int *lda,
+             const std::complex<double> *b,
+             const int *ldb,
+             const std::complex<double> *beta,
+             std::complex<double> *c,
+             const int *ldc)
+{
+    GlobalV::ofs_info << std::defaultfloat << "zgemm " << *transa << " " << *transb << " " << *m << " " << *n << " "
+                      << *k << " " << *alpha << " " << *lda << " " << *ldb << " " << *beta << " " << *ldc << std::endl;
+    zgemm_(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
+}
+
+void zaxpy_i(const int *N,
+             const std::complex<double> *alpha,
+             const std::complex<double> *X,
+             const int *incX,
+             std::complex<double> *Y,
+             const int *incY)
+{
+    // std::cout << "zaxpy " << *N << std::endl;
+    // alpha is a coefficient
+    // incX, incY is always 1
+    zaxpy_(N, alpha, X, incX, Y, incY);
+}
+
+void zhegvx_i(const int *itype,
+              const char *jobz,
+              const char *range,
+              const char *uplo,
+              const int *n,
+              std::complex<double> *a,
+              const int *lda,
+              std::complex<double> *b,
+              const int *ldb,
+              const double *vl,
+              const double *vu,
+              const int *il,
+              const int *iu,
+              const double *abstol,
+              const int *m,
+              double *w,
+              std::complex<double> *z,
+              const int *ldz,
+              std::complex<double> *work,
+              const int *lwork,
+              double *rwork,
+              int *iwork,
+              int *ifail,
+              int *info)
+{
+    GlobalV::ofs_info << std::defaultfloat << "zhegvx " << *itype << " " << *jobz << " " << *range << " " << *uplo
+                      << " " << *n << " " << *lda << " " << *ldb << " " << *vl << " " << *vu << " " << *il << " " << *iu
+                      << " " << *abstol << " " << *m << " " << *lwork << " " << *info << std::endl;
+    zhegvx_(itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork,
+            iwork, ifail, info);
+}