Added code to use MAGMA GPU library& makefile support for PGI compiler. Not working yet!

cmake/Architectures.txt

@@ -14,7 +14,12 @@ elseif ($ENV{FC} MATCHES "gfortran.*")
 elseif (${sitename} MATCHES "iff.*")
    set (configfile "cmake/cmake.iff.config")
 elseif (${sitename} MATCHES "jrl.*")
-   set (configfile "cmake/cmake.jureca.config")
+   message("Identified jureca machine. Set GPU=yes if you want to use GPU version")
+   if ($ENV{GPU} MATCHES "yes")
+      set (configfile "cmake/cmake.jurecaGPU.config")
+   else()
+	set (configfile "cmake/cmake.jureca.config")
+   endif()
 elseif (${sitename} MATCHES "jj2.*")
    set (configfile "cmake/cmake.juropa.config")
 elseif (${sitename} MATCHES "j3.*")

cmake/cmake.jurecaGPU.config

+#Cmakefile for jureca with PGI compiler and GPU support
+
+set(CMAKE_Fortran_COMPILER pgfortran)
+enable_language(C Fortran)
+set(Fleur_uses_serial TRUE)
+
+
+find_package(LibXml2 REQUIRED)
+set(CMAKE_C_FLAGS "-I${LIBXML2_INCLUDE_DIR}")
+
+if ((NOT DEFINED ENV{NO_HDF}) AND (DEFINED ENV{HDF5_ROOT})) 
+  set (HDF5_INCL "-I$ENV{HDF5_ROOT}/include")
+  set(HDF5_LIBS "-L$ENV{HDF5_ROOT}/lib -lhdf5_fortran -lhdf5 -lz")
+  set(Fleur_uses_HDF5 TRUE)
+  add_definitions(-DCPP_HDF)
+else ()
+  message("No hdf support found")
+  message("You might want to set HDF5_ROOT")
+endif ()
+
+#At the moment we use a openmp+GPU version only
+set(Fleur_uses_MPI FALSE)
+
+set(LAPACK_LIBS "-Meh_frame -lxml2 -L/homec/jiff13/jiff1302/jureca/magma -lmagma -lstdc++ -lgcc_s -lblas -Mcudalib=cublas -Mcudalib=cusparse  -L/usr/local/software/jureca/Stages/2016a/software/imkl/11.3.2.181-pmvapich2c-2016a-GDR/lib/intel64/ -lcusparse -liomp5 -lirc -lpthread")
+set(N_STREAMS 3)
+set(CMAKE_Fortran_FLAGS "-I$ENV{HOME}/jureca/magma/magma-2.0.2/include/ ${HDF5_INCL} -Mr8intrinsics -Mr8 -Mcuda:kepler+ -ta:tesla:cuda7.5 -DUSE_STREAMS -DNUM_STREAMS=${N_STREAMS} -Minfo=accel -acc " )
+set(CMAKE_Fortran_FLAGS_RELEASE " -O3 " )    
+set(CMAKE_Fortran_FLAGS_DEBUG   " -O0 -g -C -Mchkstk -Mchkptr -traceback" )
+
+if (DEFINED ENV{SCALAPACK_ROOT})
+   set(Fleur_uses_SCALAPACK TRUE)
+   set(LAPACK_LIBS "${LAPACK_LIBS} -L${SCALAPACK_ROOT}/lib  -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64")
+endif() 
+
+if (DEFINED ENV{ELPA_ROOT})
+    set(Fleur_uses_ELPA TRUE)
+    add_definitions(-DCPP_ELPA_NEW)
+    set(LAPACK_LIBS "${LAPACK_LIBS} -L$ENV{ELPA_ROOT}/lib -lelpa_openmp -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64 -lstdc++")
+    set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -I$ENV{ELPA_ROOT}/include/elpa_openmp-2014.06.001/modules")
+endif()
+
+if (DEFINED ENV{ELEMENTAL_ROOT})
+   set(Fleur_uses_ELEMENTAL TRUE)
+   set(LAPACK_LIBS "${LAPACK_LIBS} fleur_elemental.o -L$ENV{ELEMENTAL_ROOT}/lib -lelemental -lpmrrr -llapack-addons -lstdc++") 
+endif()
+
+
+
+add_definitions(-DCPP_IO -DCPP_DOUBLE -DCPP_AIX -DCPP_MAGMA)

diagonalization/CMakeLists.txt

@@ -6,7 +6,8 @@ diagonalization/franza.F90
 diagonalization/geneigbrobl.F90
 diagonalization/locrectify.F90
 diagonalization/mod_loccoeff.f90
-diagonalization/zsymsecloc.F90)
+diagonalization/zsymsecloc.F90
+diagonalization/magma.F90 )
 
 if (${Fleur_uses_MPI})
    set(fleur_F90 ${fleur_F90}

diagonalization/eigen_diag.F90

@@ -25,6 +25,11 @@ MODULE m_eigen_diag
   INTEGER,PARAMETER:: diag_scalapack=3
 #else
   INTEGER,PARAMETER:: diag_scalapack=-3
+#endif
+#ifdef CPP_MAGMA
+  INTEGER,PARAMETER:: diag_magma=6
+#else
+  INTEGER,PARAMETER:: diag_magma=-6
 #endif
   INTEGER,PARAMETER:: diag_lapack=4
   INTEGER,PARAMETER:: diag_lapack2=5
@@ -36,6 +41,9 @@ CONTAINS
     USE m_alinemuff
     USE m_types
     USE m_franza
+#ifdef CPP_MAGMA
+    USE m_magma
+#endif
    
     IMPLICIT NONE
 #ifdef CPP_MPI    
@@ -141,6 +149,10 @@ CONTAINS
 #ifdef CPP_SCALAPACK
        CASE (diag_scalapack)
           CALL chani(lapw%nmat,dimension%nbasfcn/n_size,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm, a,b,z,eig,ne_found)
+#endif
+#ifdef CPP_MAGMA
+       CASE (diag_magma)
+          CALL magma_diag(lapw%nmat,a,b,z,eig,ne_found)
 #endif
        CASE (diag_lapack2)
           if (noco%l_ss) call juDFT_error("zsymsecloc not tested with noco%l_ss")
@@ -187,7 +199,11 @@ CONTAINS
        diag_solver=diag_scalapack
 #endif
     ELSE
+#ifdef CPP_MAGMA
+       diag_solver=diag_magma
+#else
        diag_solver=diag_lapack
+#endif
     ENDIF
 
     !check if a special solver was requested
@@ -196,7 +212,7 @@ CONTAINS
     IF (juDFT_was_argument("-elemental")) diag_solver=diag_elemental
     IF (juDFT_was_argument("-lapack")) diag_solver=diag_lapack
     IF (juDFT_was_argument("-lapack2")) diag_solver=diag_lapack2
-
+    IF (juDFT_was_argument("-magma")) diag_solver=diag_magma
     !Check if solver is possible
     if (diag_solver<0) call priv_solver_error(diag_solver,parallel)
     if (any((/diag_lapack,diag_lapack2/)==diag_solver).and.parallel) call priv_solver_error(diag_solver,parallel)
@@ -236,6 +252,8 @@ CONTAINS
        IF (parallel) THEN
           CALL juDFT_error("The LAPACK2 solver can not be used in parallel")
        ENDIF  
+    CASE (diag_magma)
+       CALL juDFT_error("You have not compiled with MAGMA support")
     CASE DEFAULT
        CALL judft_error("You have selected an unkown eigensolver")
     END SELECT

diagonalization/magma.F90

+MODULE m_magma
+  use m_juDFT
+  INTEGER,PARAMETER :: NGPU_CONST=1
+  !**********************************************************
+  !     Solve the generalized eigenvalue problem
+  !     using the MAGMA library for multiple GPUs
+  !**********************************************************
+CONTAINS
+  SUBROUTINE magma_diag(nsize,a,b,z,eig,ne)
+    use magma
+#include"cpp_double.h"
+    IMPLICIT NONE
+
+    ! ... Arguments ...
+
+    INTEGER, INTENT (IN) :: nsize
+  
+    REAL,    INTENT(OUT) :: eig(:)
+    INTEGER, INTENT(INOUT) :: ne
+#ifdef CPP_INVERSION
+    REAL, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:)
+    REAL, ALLOCATABLE, INTENT (INOUT) :: z(:,:)
+#else
+    COMPLEX, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:)
+    COMPLEX, ALLOCATABLE, INTENT (INOUT) :: z(:,:)
+#endif
+
+
+    ! ... Local Variables ..
+    INTEGER iind,ind1,ind2,info,lwork,liwork,lrwork,err,i,mout(1)
+    REAL eigTemp(nsize)
+    LOGICAL:: initialized=.false.
+
+    REAL,    ALLOCATABLE :: rwork(:)
+    INTEGER, ALLOCATABLE :: iwork(:)
+#ifdef CPP_INVERSION
+    REAL, ALLOCATABLE :: largea(:,:),largeb(:,:)
+#else
+    COMPLEX, ALLOCATABLE :: largea(:,:),largeb(:,:)
+    COMPLEX,ALLOCATABLE :: work(:)
+#endif
+
+    print *,"MAGMA start"
+    IF (.NOT.initialized) THEN
+       initialized=.true.
+       call magmaf_init()
+       print *,"MAGMA init"
+    ENDIF
+
+    !**********************************
+    !expand from packed to full storage
+    !**********************************
+    !hamiltonian
+    ALLOCATE ( largea(nsize,nsize), stat=err )
+    IF (err/=0)  CALL juDFT_error("error allocating largea",calledby="geneigprobl")
+    iind = 0
+    DO ind1 = 1, nsize
+       DO ind2 = 1, ind1
+          iind = iind+1
+          largea(ind2,ind1) = a(iind)
+       ENDDO
+    ENDDO
+    !save some storage by deallocation of unused array
+    !DEALLOCATE (a)
+    !metric
+    ALLOCATE ( largeb(nsize,nsize), stat=err )
+    IF (err/=0)  CALL juDFT_error("error allocating largeb",calledby ="geneigprobl")
+    iind=0
+    DO ind1 = 1, nsize
+       DO ind2 = 1, ind1
+          iind = iind+1
+          largeb(ind2,ind1) = b(iind)
+       ENDDO
+    ENDDO
+    !save some storage by deallocation of unused array
+    !DEALLOCATE (b)
+
+#ifdef CPP_INVERSION
+    call juDFT_error("REAL diagonalization not implemented in magma.F90")
+#else
+    !Query the workspace size 
+    allocate(work(1),rwork(1),iwork(1))
+    call magmaf_zhegvdx_2stage_m(NGPU_CONST,1,MagmaVec,MagmaRangeI,MagmaLower,nsize,largea,nsize,largeb,nsize,&
+         0.0,0.0,1,ne,mout,eigTemp,work,-1,rwork,-1,iwork,-1,err)
+    lwork=work(1)
+    lrwork=rwork(1)
+    liwork=iwork(1)
+    print*,"MAGMA:",lwork,lrwork,liwork
+    deallocate(work,rwork,iwork)
+    allocate(work(lwork),rwork(lrwork),iwork(liwork))
+    if (err/=0) call juDFT_error("Failed to allocate workspaces",calledby="magma.F90")
+    !Now the diagonalization
+    call magmaf_zhegvdx_2stage_m(NGPU_CONST,1,MagmaVec,MagmaRangeI,MagmaLower,nsize,largea,nsize,largeb,nsize,&
+         0.0,0.0,1,ne,mout,eigTemp,work,lwork,rwork,lrwork,iwork,liwork,err)
+    print*,"MAGMA info:",err
+    if (err/=0) call juDFT_error("Magma failed to diagonalize Hamiltonian")
+    print *,"MAGMA mout:",mout
+#endif
+
+    DO i = 1, ne
+       eig(i) = eigTemp(i)
+    END DO
+
+  END SUBROUTINE magma_diag
+END MODULE m_magma
+