Diagonalizer gpu.v2

CMake/Modules/FindPython.cmake

@@ -3,11 +3,11 @@ UNSET(PYTHON_LIBRARIES CACHE)
 FIND_PATH(
 	PYTHON_INCLUDES
 	NAMES Python.h
-	PATH_SUFFIXES python3.0 python3.1 python3.2 python3.3 python3.4 python3.5 python3.6 python3.7 python3.8 python3.9 python3.0m python3.1m python3.2m python3.3m python3.4m python3.5m python3.6m python3.7m python3.8m python3.9m python3.0d python3.1d python3.2d python3.3d python3.4d python3.5d python3.6d python3.7d python3.8d python3.9d python3.0u python3.1u python3.2u python3.3u python3.4u python3.5u python3.6u python3.7u python3.8u python3.9u
+	PATH_SUFFIXES python3.0 python3.1 python3.2 python3.3 python3.4 python3.5 python3.6 python3.7 python3.8 python3.9 python3.10 python3.0m python3.1m python3.2m python3.3m python3.4m python3.5m python3.6m python3.7m python3.8m python3.9m python3.0d python3.1d python3.2d python3.3d python3.4d python3.5d python3.6d python3.7d python3.8d python3.9d python3.0u python3.1u python3.2u python3.3u python3.4u python3.5u python3.6u python3.7u python3.8u python3.9u
 )
 FIND_LIBRARY(
 	PYTHON_LIBRARIES
-	NAMES python3.0 python3.1 python3.2 python3.3 python3.4 python3.5 python3.6 python3.7 python3.8 python3.9 python3.0m python3.1m python3.2m python3.3m python3.4m python3.5m python3.6m python3.7m python3.8m python3.9m python3.0d python3.1d python3.2d python3.3d python3.4d python3.5d python3.6d python3.7d python3.8d python3.9d python3.0u python3.1u python3.2u python3.3u python3.4u python3.5u python3.6u python3.7u python3.8u python3.9u
+	NAMES python3.0 python3.1 python3.2 python3.3 python3.4 python3.5 python3.6 python3.7 python3.8 python3.9 python3.10 python3.0m python3.1m python3.2m python3.3m python3.4m python3.5m python3.6m python3.7m python3.8m python3.9m python3.0d python3.1d python3.2d python3.3d python3.4d python3.5d python3.6d python3.7d python3.8d python3.9d python3.0u python3.1u python3.2u python3.3u python3.4u python3.5u python3.6u python3.7u python3.8u python3.9u
 	PATH_SUFFIXES lib lib32 lib64
 )
 IF(("${PYTHON_INCLUDES}" MATCHES "PYTHON_INCLUDES-NOTFOUND") OR ("${PYTHON_LIBRARIES}" MATCHES "PYTHON_LIBRARIES-NOTFOUND"))

CMake/TBTKConfigTemplate.cmake

@@ -64,6 +64,7 @@ IF(CUDA_FOUND)
 	TBTK_MESSAGE("[X] CUDA")
 	LIST(APPEND TBTK_LIBRARIES ${CUDA_LIBRARIES})
 	LIST(APPEND TBTK_LIBRARIES ${CUDA_cusparse_LIBRARY})
+	LIST(APPEND TBTK_LIBRARIES ${CUDA_cusolver_LIBRARY})
 	ADD_DEFINITIONS(-DTBTK_CUDA_ENABLED)
 ELSE(CUDA_FOUND)
 	TBTK_MESSAGE("[ ] CUDA")

CMakeLists.txt

@@ -306,7 +306,7 @@ ENDIF(TBTK_ENABLE_ALL_OPTIONS)
 MESSAGE("================================================================================")
 
 #C++ flags
-SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -O3")
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -Wall -O3")
 IF("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang")
 	SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-undefined-var-template")
 ENDIF("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang")

Lib/include/TBTK/Solver/Diagonalizer.h

@@ -92,6 +92,12 @@ class Diagonalizer : public Solver, public Communicator{
 	 *  self-consistent callculation. */
 	void setMaxIterations(int maxIterations);
 
+	/** Set if you want to use GPU acceleration provided by CUDA routines.
+	 *
+	 *  @param useGPUAcceleration Turns GPU acceleration for the Diagonalizer 
+	 *  solver off or on. */
+	void setUseGPUAcceleration(bool useGPUAcceleration);
+
 	/** Run calculations. Diagonalizes ones if no self-consistency callback
 	 *  have been set, or otherwise multiple times until self-consistencey
 	 *  or maximum number of iterations has been reached. */
@@ -132,7 +138,7 @@ class Diagonalizer : public Solver, public Communicator{
 	 *  @param state The state number, ordered in accending order.
 	 *
 	 *  @return The eigenvalue for the given state. */
-	const double getEigenValue(int state);
+	double getEigenValue(int state);
 
 	/** Get amplitude for given eigenvector \f$n\f$ and physical index
 	 * \f$x\f$: \f$\Psi_{n}(x)\f$.
@@ -159,6 +165,9 @@ class Diagonalizer : public Solver, public Communicator{
 	/** Maximum number of iterations in the self-consistency loop. */
 	int maxIterations;
 
+	/** Enables GPU acceleration for the solver. */
+	bool useGPUAcceleration;
+
 	/** SelfConsistencyCallback to call each time a diagonalization has
 	 *  been completed. */
 	SelfConsistencyCallback *selfConsistencyCallback;
@@ -172,12 +181,29 @@ class Diagonalizer : public Solver, public Communicator{
 	/** Diagonalizes the Hamiltonian. */
 	void solve();
 
+	/** Diagonalizes the input matrix using a GPU device.
+	 *  The output for the eigen vectors is
+	 *  written into the input matrix. 
+	 * 
+	 *  @param matrix Matrix nxn to diagonalize
+	 *  @param eigenValues vector of size n
+	 *
+	 *  @return Function overwrites matrix with eigenvectors
+	 * 			and returns eigenvalues in eigenValues     */
+	void solveGPU(CArray<std::complex<double>>& matrix, CArray<double>& eigenValues);
+
 	/** Setup the basis transformation. */
 	void setupBasisTransformation();
 
+	/** Setup the basis transformation on GPU. */
+	void setupBasisTransformationGPU();
+
 	/** Transform the Hamiltonian to an orthonormal basis. */
 	void transformToOrthonormalBasis();
 
+	/** Transform the Hamiltonian to an orthonormal basis on GPU. */	
+	void transformToOrthonormalBasisGPU();
+
 	/** Transform the eigen vectors to the original basis. */
 	void transformToOriginalBasis();
 };
@@ -216,10 +242,14 @@ inline const std::complex<double> Diagonalizer::getAmplitude(
 	return eigenVectors[model.getBasisSize()*state + model.getBasisIndex(index)];
 }
 
-inline const double Diagonalizer::getEigenValue(int state){
+inline double Diagonalizer::getEigenValue(int state){
 	return eigenValues[state];
 }
 
+inline void Diagonalizer::setUseGPUAcceleration(bool useGPUAcceleration){
+	this->useGPUAcceleration = useGPUAcceleration;
+}
+
 };	//End of namespace Solver
 };	//End of namespace TBTK
 

Lib/include/Utilities/TBTK/Exporter.h

@@ -183,17 +183,27 @@ class Exporter{
 	 *
 	 *  @param format The output format to use.*/
 	void setFormat(Format format);
+
+	/** Sets the number of significant digits to use in the save
+	 * file. The default is 6.
+	 *
+	 *  @param numberDigits The number of significant digits to use.*/
+	void setNumberSiginificantDigits(const unsigned& numberDigits);
 private:
 	/** The output format. */
 	Format format;
 
+	/** sets the number of significant digits */
+	unsigned numberSiginificantDigits;
+
 	/** Write double to output stream. */
 	template<typename DataType>
 	void write(std::ofstream &stream, const DataType &value) const;
 };
 
 inline Exporter::Exporter(){
 	format = Format::RowMajor;
+	numberSiginificantDigits = 6; // C++ standard
 }
 
 template<typename DataType>
@@ -241,6 +251,7 @@ void Exporter::save(
 	}
 
 	std::ofstream fout(filename);
+	fout << std::setprecision(numberSiginificantDigits);
 	if(!fout){
 		TBTKExit(
 			"Exporter::save()",
@@ -307,6 +318,10 @@ inline void Exporter::write<SpinMatrix>(
 	}
 }
 
+inline void Exporter::setNumberSiginificantDigits(const unsigned& numberDigits){
+	numberSiginificantDigits = numberDigits;
+}
+
 }; //End of namesapce TBTK
 
 #endif

Lib/src/Solver/Diagonalizer.cpp

@@ -30,6 +30,7 @@ namespace Solver{
 Diagonalizer::Diagonalizer() : Communicator(false){
 	maxIterations = 50;
 	selfConsistencyCallback = nullptr;
+	useGPUAcceleration = false;
 }
 
 void Diagonalizer::run(){
@@ -47,7 +48,14 @@ void Diagonalizer::run(){
 			Streams::out << "." << flush;
 		}
 
-		solve();
+		if(useGPUAcceleration){
+			eigenVectors = hamiltonian;
+			solveGPU(eigenVectors, eigenValues);
+			transformToOriginalBasis();
+		}
+		else{
+			solve();
+		}
 
 		if(selfConsistencyCallback){
 			if(selfConsistencyCallback->selfConsistencyCallback(*this))
@@ -71,7 +79,13 @@ void Diagonalizer::init(){
 	if(getGlobalVerbose() && getVerbose())
 		Streams::out << "\tBasis size: " << basisSize << "\n";
 
-	hamiltonian = CArray<complex<double>>((basisSize*(basisSize+1))/2);
+	if(useGPUAcceleration){
+		hamiltonian = CArray<complex<double>>(basisSize*basisSize);
+	}
+	else{
+		hamiltonian = CArray<complex<double>>((basisSize*(basisSize+1))/2);
+	}
+
 	eigenValues = CArray<double>(basisSize);
 	eigenVectors = CArray<complex<double>>(basisSize*basisSize);
 
@@ -97,12 +111,24 @@ void Diagonalizer::update(){
 		int to = model.getHoppingAmplitudeSet().getBasisIndex(
 			(*iterator).getToIndex()
 		);
-		if(from >= to)
-			hamiltonian[to + (from*(from+1))/2] += (*iterator).getAmplitude();
+		if(from >= to){
+			if(useGPUAcceleration){
+				hamiltonian[to + from*basisSize] += (*iterator).getAmplitude();
+			}
+			else{
+				hamiltonian[to + (from*(from+1))/2] += (*iterator).getAmplitude();
+			}
+		}
 	}
 
-	setupBasisTransformation();
-	transformToOrthonormalBasis();
+	if(useGPUAcceleration){
+		setupBasisTransformationGPU();
+		transformToOrthonormalBasisGPU();
+	}
+	else{
+		setupBasisTransformation();
+		transformToOrthonormalBasis();
+	}
 }
 
 //Lapack function for matrix diagonalization of triangular matrix.
@@ -118,7 +144,7 @@ extern "C" void zhpev_(char *jobz,		//'E' = Eigenvalues only, 'V' = Eigenvalues
 			int *info);		//0 = successful, <0 = -info value was illegal, >0 = info number of off-diagonal elements failed to converge.
 
 //Lapack function for matrix diagonalization of banded triangular matrix
-extern "C" void zhbeb_(
+extern "C" void zhbeb_( //TODO this function is not used?
 	char *jobz,		//'E' = Eigenvalues only, 'V' = Eigenvalues and eigenvectors.
 	char *uplo,		//'U' = Stored as upper triangular, 'L' = Stored as lower triangular.
 	int *n,			//n*n = Matrix size