elpa2.F90

!   This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaftrn,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!
!    More information can be found here:
!    http://elpa.rzg.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
!
! ELPA1 -- Faster replacements for ScaLAPACK symmetric eigenvalue routines
!
! Copyright of the original code rests with the authors inside the ELPA
! consortium. The copyright of any additional modifications shall rest
! with their original authors, but shall adhere to the licensing terms
! distributed along with the original code in the file "COPYING".



! ELPA2 -- 2-stage solver for ELPA
!
! Copyright of the original code rests with the authors inside the ELPA
! consortium. The copyright of any additional modifications shall rest
! with their original authors, but shall adhere to the licensing terms
! distributed along with the original code in the file "COPYING".


#include "config-f90.h"

module ELPA2

! Version 1.1.2, 2011-02-21

  use elpa_utilities
  USE ELPA1
  use elpa2_utilities
  use elpa_pdgeqrf

  use iso_c_binding

#ifdef WITH_GPU_VERSION
!  use cuda_routines
!  use cuda_c_kernel
!  use iso_c_binding
#endif
  implicit none

  PRIVATE ! By default, all routines contained are private

  ! The following routines are public:

  public :: solve_evp_real_2stage
  public :: solve_evp_complex_2stage

  public :: bandred_real
  public :: tridiag_band_real
  public :: trans_ev_tridi_to_band_real
  public :: trans_ev_band_to_full_real

  public :: bandred_complex
  public :: tridiag_band_complex
  public :: trans_ev_tridi_to_band_complex
  public :: trans_ev_band_to_full_complex

  public :: band_band_real
  public :: divide_band

  integer, public :: which_qr_decomposition = 1     ! defines, which QR-decomposition algorithm will be used
                                                    ! 0 for unblocked
                                                    ! 1 for blocked (maxrank: nblk)
!-------------------------------------------------------------------------------

  ! The following array contains the Householder vectors of the
  ! transformation band -> tridiagonal.
  ! It is allocated and set in tridiag_band_real and used in
  ! trans_ev_tridi_to_band_real.
  ! It must be deallocated by the user after trans_ev_tridi_to_band_real!

  real*8, allocatable :: hh_trans_real(:,:)
  complex*16, allocatable :: hh_trans_complex(:,:)

!-------------------------------------------------------------------------------

  include 'mpif.h'


!******
contains

  function solve_evp_real_2stage(na, nev, a, lda, ev, q, ldq, nblk,        &
                               matrixCols,                               &
                                 mpi_comm_rows, mpi_comm_cols,           &
                                 mpi_comm_all, THIS_REAL_ELPA_KERNEL_API,&
                                 useQR) result(success)

  !-------------------------------------------------------------------------------
  !  solve_evp_real_2stage: Solves the real eigenvalue problem with a 2 stage approach
  !
  !  Parameters
  !
  !  na          Order of matrix a
  !
  !  nev         Number of eigenvalues needed
  !
  !  a(lda,matrixCols)    Distributed matrix for which eigenvalues are to be computed.
  !              Distribution is like in Scalapack.
  !              The full matrix must be set (not only one half like in scalapack).
  !              Destroyed on exit (upper and lower half).
  !
  !  lda         Leading dimension of a
  !  matrixCols  local columns of matrix a and q
  !
  !  ev(na)      On output: eigenvalues of a, every processor gets the complete set
  !
  !  q(ldq,matrixCols)    On output: Eigenvectors of a
  !              Distribution is like in Scalapack.
  !              Must be always dimensioned to the full size (corresponding to (na,na))
  !              even if only a part of the eigenvalues is needed.
  !
  !  ldq         Leading dimension of q
  !
  !  nblk        blocksize of cyclic distribution, must be the same in both directions!
  !
  !  mpi_comm_rows
  !  mpi_comm_cols
  !              MPI-Communicators for rows/columns
  !  mpi_comm_all
  !              MPI-Communicator for the total processor set
  !
  !-------------------------------------------------------------------------------
#ifdef HAVE_DETAILED_TIMINGS
    use timings
#endif
    use cuda_functions
    use mod_check_for_gpu

    implicit none
    logical, intent(in), optional :: useQR
    logical                       :: useQRActual, useQREnvironment
    integer, intent(in), optional :: THIS_REAL_ELPA_KERNEL_API
    integer                       :: THIS_REAL_ELPA_KERNEL

    integer, intent(in)           :: na, nev, lda, ldq, matrixCols, mpi_comm_rows, &
                                     mpi_comm_cols, mpi_comm_all
    integer, intent(in)           :: nblk
    real*8, intent(inout)         :: a(lda,matrixCols), ev(na), q(ldq,matrixCols)

    integer                       :: my_pe, n_pes, my_prow, my_pcol, np_rows, np_cols, mpierr
    integer                       :: nbw, num_blocks
    real*8, allocatable           :: tmat(:,:,:), e(:)
    real*8                        :: ttt0, ttt1, ttts
    integer                       :: i
    logical                       :: success
    logical, save                 :: firstCall = .true.
    logical                       :: wantDebug
    integer                       :: istat
    character(200)                :: errorMessage
    logical                       :: useGPU
    integer                       :: numberOfGPUDevices


#ifdef HAVE_DETAILED_TIMINGS
    call timer%start("solve_evp_real_2stage")
#endif
    call mpi_comm_rank(mpi_comm_all,my_pe,mpierr)
    call mpi_comm_size(mpi_comm_all,n_pes,mpierr)

    call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
    call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)
    call mpi_comm_rank(mpi_comm_cols,my_pcol,mpierr)
    call mpi_comm_size(mpi_comm_cols,np_cols,mpierr)


    wantDebug = .false.
    if (firstCall) then
      ! are debug messages desired?
      wantDebug = debug_messages_via_environment_variable()
      firstCall = .false.
    endif

    success = .true.

    useQRActual = .false.
    useGPU      = .false.

    ! set usage of qr decomposition via API call
    if (present(useQR)) then
      if (useQR) useQRActual = .true.
        if (.not.(useQR)) useQRACtual = .false.
    endif

    ! overwrite this with environment variable settings
    if (qr_decomposition_via_environment_variable(useQREnvironment)) then
      useQRActual = useQREnvironment
    endif

    if (useQRActual) then
      if (mod(na,nblk) .ne. 0) then
        if (wantDebug) then
          write(error_unit,*) "solve_evp_real_2stage: QR-decomposition: blocksize does not fit with matrixsize"
        endif
        print *, "Do not use QR-decomposition for this matrix and blocksize."
        success = .false.
        return
      endif
    endif


    if (present(THIS_REAL_ELPA_KERNEL_API)) then
      ! user defined kernel via the optional argument in the API call
      THIS_REAL_ELPA_KERNEL = THIS_REAL_ELPA_KERNEL_API
    else

      ! if kernel is not choosen via api
      ! check whether set by environment variable
      THIS_REAL_ELPA_KERNEL = get_actual_real_kernel()
    endif

    ! check whether choosen kernel is allowed
    if (check_allowed_real_kernels(THIS_REAL_ELPA_KERNEL)) then

      if (my_pe == 0) then
        write(error_unit,*) " "
        write(error_unit,*) "The choosen kernel ",REAL_ELPA_KERNEL_NAMES(THIS_REAL_ELPA_KERNEL)
        write(error_unit,*) "is not in the list of the allowed kernels!"
        write(error_unit,*) " "
        write(error_unit,*) "Allowed kernels are:"
        do i=1,size(REAL_ELPA_KERNEL_NAMES(:))
          if (AVAILABLE_REAL_ELPA_KERNELS(i) .ne. 0) then
            write(error_unit,*) REAL_ELPA_KERNEL_NAMES(i)
          endif
        enddo

        write(error_unit,*) " "
        write(error_unit,*) "The defaul kernel REAL_ELPA_KERNEL_GENERIC will be used !"
      endif
      THIS_REAL_ELPA_KERNEL = REAL_ELPA_KERNEL_GENERIC

    endif

    if (THIS_REAL_ELPA_KERNEL .eq. REAL_ELPA_KERNEL_GPU) then
      if (check_for_gpu(my_pe,numberOfGPUDevices)) then
        useGPU = .true.
      endif
      if (nblk .ne. 128) then
        print *,"At the moment GPU version needs blocksize 128"
        stop
      endif

      ! set the neccessary parameters
      cudaMemcpyHostToDevice   = cuda_memcpyHostToDevice()
      cudaMemcpyDeviceToHost   = cuda_memcpyDeviceToHost()
      cudaMemcpyDeviceToDevice = cuda_memcpyDeviceToDevice()
      cudaHostRegisterPortable = cuda_hostRegisterPortable()
      cudaHostRegisterMapped   = cuda_hostRegisterMapped()
    endif

    ! Choose bandwidth, must be a multiple of nblk, set to a value >= 32

    if (useGPU) then
      nbw = nblk
    else
      nbw = (31/nblk+1)*nblk
    endif

    num_blocks = (na-1)/nbw + 1

    allocate(tmat(nbw,nbw,num_blocks), stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_real_2stage: error when allocating tmat "//errorMessage
      stop
    endif

    ! Reduction full -> band

    ttt0 = MPI_Wtime()
    ttts = ttt0
    call bandred_real(na, a, lda, nblk, nbw, matrixCols, num_blocks, mpi_comm_rows, mpi_comm_cols, &
                      tmat, wantDebug, useGPU, success, useQRActual)
    if (.not.(success)) return
    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time bandred_real               :',ttt1-ttt0

     ! Reduction band -> tridiagonal

     allocate(e(na), stat=istat, errmsg=errorMessage)
     if (istat .ne. 0) then
       print *,"solve_evp_real_2stage: error when allocating e "//errorMessage
       stop
     endif

     ttt0 = MPI_Wtime()
     call tridiag_band_real(na, nbw, nblk, a, lda, ev, e, matrixCols, mpi_comm_rows, &
                            mpi_comm_cols, mpi_comm_all)
     ttt1 = MPI_Wtime()
     if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time tridiag_band_real          :',ttt1-ttt0

     call mpi_bcast(ev,na,MPI_REAL8,0,mpi_comm_all,mpierr)
     call mpi_bcast(e,na,MPI_REAL8,0,mpi_comm_all,mpierr)

     ttt1 = MPI_Wtime()
     time_evp_fwd = ttt1-ttts

     ! Solve tridiagonal system

     ttt0 = MPI_Wtime()
     call solve_tridi(na, nev, ev, e, q, ldq, nblk, matrixCols, mpi_comm_rows,  &
                      mpi_comm_cols, wantDebug, success)
     if (.not.(success)) return

     ttt1 = MPI_Wtime()
     if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
     write(error_unit,*) 'Time solve_tridi                :',ttt1-ttt0
     time_evp_solve = ttt1-ttt0
     ttts = ttt1

     deallocate(e, stat=istat, errmsg=errorMessage)
     if (istat .ne. 0) then
       print *,"solve_evp_real_2stage: error when deallocating e "//errorMessage
       stop
     endif
     ! Backtransform stage 1

     ttt0 = MPI_Wtime()

     call trans_ev_tridi_to_band_real(na, nev, nblk, nbw, q, ldq, matrixCols, mpi_comm_rows, &
                                      mpi_comm_cols, wantDebug, useGPU, success, THIS_REAL_ELPA_KERNEL)
     if (.not.(success)) return
     ttt1 = MPI_Wtime()
     if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time trans_ev_tridi_to_band_real:',ttt1-ttt0

     ! We can now deallocate the stored householder vectors
     deallocate(hh_trans_real, stat=istat, errmsg=errorMessage)
     if (istat .ne. 0) then
       print *,"solve_evp_real_2stage: error when deallocating hh_trans_real "//errorMessage
       stop
     endif


     ! Backtransform stage 2
     print *,"useGPU== ",useGPU
     ttt0 = MPI_Wtime()
     call trans_ev_band_to_full_real(na, nev, nblk, nbw, a, lda, tmat, q, ldq, matrixCols, num_blocks, mpi_comm_rows, &
                                     mpi_comm_cols, useGPU, useQRActual)
     ttt1 = MPI_Wtime()
     if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time trans_ev_band_to_full_real :',ttt1-ttt0
     time_evp_back = ttt1-ttts

     deallocate(tmat, stat=istat, errmsg=errorMessage)
     if (istat .ne. 0) then
       print *,"solve_evp_real_2stage: error when deallocating tmat"//errorMessage
       stop
     endif

#ifdef HAVE_DETAILED_TIMINGS
     call timer%stop("solve_evp_real_2stage")
#endif
1    format(a,f10.3)

   end function solve_evp_real_2stage

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

  function solve_evp_complex_2stage(na, nev, a, lda, ev, q, ldq, nblk, &
                                  matrixCols, mpi_comm_rows, mpi_comm_cols,      &
                                    mpi_comm_all, THIS_COMPLEX_ELPA_KERNEL_API) result(success)

  !-------------------------------------------------------------------------------
  !  solve_evp_complex_2stage: Solves the complex eigenvalue problem with a 2 stage approach
  !
  !  Parameters
  !
  !  na          Order of matrix a
  !
  !  nev         Number of eigenvalues needed
  !
  !  a(lda,matrixCols)    Distributed matrix for which eigenvalues are to be computed.
  !              Distribution is like in Scalapack.
  !              The full matrix must be set (not only one half like in scalapack).
  !              Destroyed on exit (upper and lower half).
  !
  !  lda         Leading dimension of a
  !  matrixCols  local columns of matrix a and q
  !
  !  ev(na)      On output: eigenvalues of a, every processor gets the complete set
  !
  !  q(ldq,matrixCols)    On output: Eigenvectors of a
  !              Distribution is like in Scalapack.
  !              Must be always dimensioned to the full size (corresponding to (na,na))
  !              even if only a part of the eigenvalues is needed.
  !
  !  ldq         Leading dimension of q
  !
  !  nblk        blocksize of cyclic distribution, must be the same in both directions!
  !
  !  mpi_comm_rows
  !  mpi_comm_cols
  !              MPI-Communicators for rows/columns
  !  mpi_comm_all
  !              MPI-Communicator for the total processor set
  !
  !-------------------------------------------------------------------------------
#ifdef HAVE_DETAILED_TIMINGS
    use timings
#endif
    use cuda_functions
    use mod_check_for_gpu
    implicit none
    integer, intent(in), optional :: THIS_COMPLEX_ELPA_KERNEL_API
    integer                       :: THIS_COMPLEX_ELPA_KERNEL
    integer, intent(in)           :: na, nev, lda, ldq, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols, mpi_comm_all
    complex*16, intent(inout)     :: a(lda,matrixCols), q(ldq,matrixCols)
    real*8, intent(inout)         :: ev(na)

    integer                       :: my_prow, my_pcol, np_rows, np_cols, mpierr, my_pe, n_pes
    integer                       :: l_cols, l_rows, l_cols_nev, nbw, num_blocks
    complex*16, allocatable       :: tmat(:,:,:)
    real*8, allocatable           :: q_real(:,:), e(:)
    real*8                        :: ttt0, ttt1, ttts
    integer                       :: i

    logical                       :: success, wantDebug
    logical, save                 :: firstCall = .true.
    integer                       :: istat
    character(200)                :: errorMessage
    logical                       :: useGPU
    integer                       :: numberOfGPUDevices


#ifdef HAVE_DETAILED_TIMINGS
    call timer%start("solve_evp_complex_2stage")
#endif
    call mpi_comm_rank(mpi_comm_all,my_pe,mpierr)
    call mpi_comm_size(mpi_comm_all,n_pes,mpierr)

    call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
    call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)
    call mpi_comm_rank(mpi_comm_cols,my_pcol,mpierr)
    call mpi_comm_size(mpi_comm_cols,np_cols,mpierr)

    useGPU = .false.
    wantDebug = .false.
    if (firstCall) then
      ! are debug messages desired?
      wantDebug = debug_messages_via_environment_variable()
      firstCall = .false.
    endif


    success = .true.

    if (present(THIS_COMPLEX_ELPA_KERNEL_API)) then
      ! user defined kernel via the optional argument in the API call
      THIS_COMPLEX_ELPA_KERNEL = THIS_COMPLEX_ELPA_KERNEL_API
    else
      ! if kernel is not choosen via api
      ! check whether set by environment variable
      THIS_COMPLEX_ELPA_KERNEL = get_actual_complex_kernel()
    endif

    ! check whether choosen kernel is allowed
    if (check_allowed_complex_kernels(THIS_COMPLEX_ELPA_KERNEL)) then

      if (my_pe == 0) then
        write(error_unit,*) " "
        write(error_unit,*) "The choosen kernel ",COMPLEX_ELPA_KERNEL_NAMES(THIS_COMPLEX_ELPA_KERNEL)
        write(error_unit,*) "is not in the list of the allowed kernels!"
        write(error_unit,*) " "
        write(error_unit,*) "Allowed kernels are:"
        do i=1,size(COMPLEX_ELPA_KERNEL_NAMES(:))
          if (AVAILABLE_COMPLEX_ELPA_KERNELS(i) .ne. 0) then
            write(error_unit,*) COMPLEX_ELPA_KERNEL_NAMES(i)
          endif
        enddo

        write(error_unit,*) " "
        write(error_unit,*) "The defaul kernel COMPLEX_ELPA_KERNEL_GENERIC will be used !"
      endif
      THIS_COMPLEX_ELPA_KERNEL = COMPLEX_ELPA_KERNEL_GENERIC
    endif

    if (THIS_COMPLEX_ELPA_KERNEL .eq. COMPLEX_ELPA_KERNEL_GPU) then
      if (check_for_gpu(my_pe, numberOfGPUDevices)) then
        useGPU=.true.
      endif
      if (nblk .ne. 128) then
        print *,"At the moment GPU version needs blocksize 128"
        stop
      endif

      ! set the neccessary parameters
      cudaMemcpyHostToDevice   = cuda_memcpyHostToDevice()
      cudaMemcpyDeviceToHost   = cuda_memcpyDeviceToHost()
      cudaMemcpyDeviceToDevice = cuda_memcpyDeviceToDevice()
      cudaHostRegisterPortable = cuda_hostRegisterPortable()
      cudaHostRegisterMapped   = cuda_hostRegisterMapped()
    endif

    ! Choose bandwidth, must be a multiple of nblk, set to a value >= 32

    nbw = (31/nblk+1)*nblk

    num_blocks = (na-1)/nbw + 1

    allocate(tmat(nbw,nbw,num_blocks), stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when allocating tmat"//errorMessage
      stop
    endif
    ! Reduction full -> band

    ttt0 = MPI_Wtime()
    ttts = ttt0
    call bandred_complex(na, a, lda, nblk, nbw, matrixCols, num_blocks, mpi_comm_rows, mpi_comm_cols, &
                         tmat, wantDebug, useGPU, success)
    if (.not.(success)) then
#ifdef HAVE_DETAILED_TIMINGS
      call timer%stop()
#endif
      return
    endif
    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time bandred_complex               :',ttt1-ttt0

    ! Reduction band -> tridiagonal

    allocate(e(na), stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when allocating e"//errorMessage
      stop
    endif


    ttt0 = MPI_Wtime()
    call tridiag_band_complex(na, nbw, nblk, a, lda, ev, e, matrixCols, mpi_comm_rows, mpi_comm_cols, mpi_comm_all)
    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time tridiag_band_complex          :',ttt1-ttt0

    call mpi_bcast(ev,na,MPI_REAL8,0,mpi_comm_all,mpierr)
    call mpi_bcast(e,na,MPI_REAL8,0,mpi_comm_all,mpierr)

    ttt1 = MPI_Wtime()
    time_evp_fwd = ttt1-ttts

    l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a and q
    l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local columns of q
    l_cols_nev = local_index(nev, my_pcol, np_cols, nblk, -1) ! Local columns corresponding to nev

    allocate(q_real(l_rows,l_cols), stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when allocating q_real"//errorMessage
      stop
    endif

    ! Solve tridiagonal system

    ttt0 = MPI_Wtime()
    call solve_tridi(na, nev, ev, e, q_real, ubound(q_real,dim=1), nblk, matrixCols, &
                     mpi_comm_rows, mpi_comm_cols, wantDebug, success)
    if (.not.(success)) return

    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times)  &
       write(error_unit,*) 'Time solve_tridi                   :',ttt1-ttt0
    time_evp_solve = ttt1-ttt0
    ttts = ttt1

    q(1:l_rows,1:l_cols_nev) = q_real(1:l_rows,1:l_cols_nev)

    deallocate(e, q_real, stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when deallocating e, q_real"//errorMessage
      stop
    endif


    ! Backtransform stage 1

    ttt0 = MPI_Wtime()
    call trans_ev_tridi_to_band_complex(na, nev, nblk, nbw, q, ldq,  &
                                        matrixCols, mpi_comm_rows, mpi_comm_cols,&
                                        wantDebug, useGPU, success,THIS_COMPLEX_ELPA_KERNEL)
    if (.not.(success)) return
    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time trans_ev_tridi_to_band_complex:',ttt1-ttt0

    ! We can now deallocate the stored householder vectors
    deallocate(hh_trans_complex, stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when deallocating hh_trans_complex"//errorMessage
      stop
    endif



    ! Backtransform stage 2

    ttt0 = MPI_Wtime()
    call trans_ev_band_to_full_complex(na, nev, nblk, nbw, a, lda, tmat, q, ldq, matrixCols, num_blocks, mpi_comm_rows, &
                                       mpi_comm_cols, useGPU)
    ttt1 = MPI_Wtime()
    if (my_prow==0 .and. my_pcol==0 .and. elpa_print_times) &
       write(error_unit,*) 'Time trans_ev_band_to_full_complex :',ttt1-ttt0
    time_evp_back = ttt1-ttts

    deallocate(tmat, stat=istat, errmsg=errorMessage)
    if (istat .ne. 0) then
      print *,"solve_evp_complex_2stage: error when deallocating tmat "//errorMessage
      stop
    endif

#ifdef HAVE_DETAILED_TIMINGS
    call timer%stop("solve_evp_complex_2stage")
#endif

1   format(a,f10.3)

  end function solve_evp_complex_2stage

!-------------------------------------------------------------------------------


  subroutine bandred_real(na, a, lda, nblk, nbw, matrixCols, numBlocks, mpi_comm_rows, mpi_comm_cols, &
                        tmat, wantDebug, useGPU, success, useQR)

  !-------------------------------------------------------------------------------
  !  bandred_real: Reduces a distributed symmetric matrix to band form
  !
  !  Parameters
  !
  !  na          Order of matrix
  !
  !  a(lda,matrixCols)    Distributed matrix which should be reduced.
  !              Distribution is like in Scalapack.
  !              Opposed to Scalapack, a(:,:) must be set completely (upper and lower half)
  !              a(:,:) is overwritten on exit with the band and the Householder vectors
  !              in the upper half.
  !
  !  lda         Leading dimension of a
  !  matrixCols  local columns of matrix a
  !
  !  nblk        blocksize of cyclic distribution, must be the same in both directions!
  !
  !  nbw         semi bandwith of output matrix
  !
  !  mpi_comm_rows
  !  mpi_comm_cols
  !              MPI-Communicators for rows/columns
  !
  !  tmat(nbw,nbw,numBlocks)    where numBlocks = (na-1)/nbw + 1
  !              Factors for the Householder vectors (returned), needed for back transformation
  !
  !-------------------------------------------------------------------------------


    use cuda_functions
    use iso_c_binding

#ifdef HAVE_DETAILED_TIMINGS
   use timings
#endif

   implicit none

   integer                  :: na, lda, nblk, nbw, matrixCols, numBlocks, mpi_comm_rows, mpi_comm_cols
   real*8                   :: a(lda,matrixCols), tmat(nbw,nbw,numBlocks)
   logical, intent(in)      :: useGPU

   integer                  :: my_prow, my_pcol, np_rows, np_cols, mpierr
   integer                  :: l_cols, l_rows
   integer                  :: i, j, lcs, lce, lre, lc, lr, cur_pcol, n_cols, nrow
   integer                  :: istep, ncol, lch, lcx, nlc
   integer                  :: tile_size, l_rows_tile, l_cols_tile

   real*8                   :: eps

   real*8                   :: vnorm2, xf, aux1(nbw), aux2(nbw), vrl, tau, vav(nbw,nbw)


   real*8, allocatable      :: tmpCUDA(:),  vmrCUDA(:),  umcCUDA(:)
   real*8, allocatable      :: tmpCPU(:,:), vmrCPU(:,:), umcCPU(:,:)
   real*8, allocatable      :: vr(:)


   ! needed for blocked QR decomposition
   integer                  :: PQRPARAM(11), work_size
   real*8                   :: dwork_size(1)
   real*8, allocatable      :: work_blocked(:), tauvector(:), blockheuristic(:)

   integer(kind=C_intptr_T) :: a_dev, vmr_dev, umc_dev, tmat_dev, vav_dev
   integer, external        :: numroc
   integer                  :: ierr
   integer                  :: cur_l_rows, cur_l_cols, vmr_size, umc_size
   integer(kind=c_size_t)   :: lc_start, lc_end
   integer                  :: lr_end
   integer                  :: na_rows, na_cols

   logical, intent(in)      :: wantDebug
   logical, intent(out)     :: success

   logical, intent(in)      :: useQR
   logical                  :: successCUDA
   integer                  :: istat
   character(200)           :: errorMessage

#ifdef HAVE_DETAILED_TIMINGS
   call timer%start("bandred_real")
#endif
   call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
   call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)
   call mpi_comm_rank(mpi_comm_cols,my_pcol,mpierr)
   call mpi_comm_size(mpi_comm_cols,np_cols,mpierr)
   success = .true.


   ! Semibandwith nbw must be a multiple of blocksize nblk
   if (mod(nbw,nblk)/=0) then
     if (my_prow==0 .and. my_pcol==0) then
       if (wantDebug) then
         write(error_unit,*) 'ELPA2_bandred_real: ERROR: nbw=',nbw,', nblk=',nblk
         write(error_unit,*) 'ELPA2_bandred_real: ELPA2 works only for nbw==n*nblk'
       endif
       success = .false.
       return
     endif
   endif

   if (useGPU) then
     na_rows = numroc(na, nblk, my_prow, 0, np_rows)
     na_cols = numroc(na, nblk, my_pcol, 0, np_cols)
   endif

   ! Matrix is split into tiles; work is done only for tiles on the diagonal or above

   tile_size = nblk*least_common_multiple(np_rows,np_cols) ! minimum global tile size
   tile_size = ((128*max(np_rows,np_cols)-1)/tile_size+1)*tile_size ! make local tiles at least 128 wide

   l_rows_tile = tile_size/np_rows ! local rows of a tile
   l_cols_tile = tile_size/np_cols ! local cols of a tile

   if (useQR) then
     if (useGPU) then
       print *,"qr decomposition at the moment not supported with GPU"
       stop
     endif

     if (which_qr_decomposition == 1) then
       call qr_pqrparam_init(pqrparam,    nblk,'M',0,   nblk,'M',0,   nblk,'M',1,'s')
       allocate(tauvector(na), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating tauvector "//errorMessage
         stop
       endif

       allocate(blockheuristic(nblk), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating blockheuristic "//errorMessage
         stop
       endif

       l_rows = local_index(na, my_prow, np_rows, nblk, -1)
       allocate(vmrCPU(max(l_rows,1),na), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating vmrCPU "//errorMessage
         stop
       endif

       call qr_pdgeqrf_2dcomm(a, lda, vmrCPU, max(l_rows,1), tauvector, tmat(1,1,1), nbw, dwork_size(1), -1, na, &
                             nbw, nblk, nblk, na, na, 1, 0, PQRPARAM, mpi_comm_rows, mpi_comm_cols, blockheuristic)
       work_size = dwork_size(1)
       allocate(work_blocked(work_size), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating work_blocked "//errorMessage
         stop
       endif

       work_blocked = 0.0d0
       deallocate(vmrCPU, stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when deallocating vmrCPU "//errorMessage
         stop
       endif

     endif

   endif ! useQr

   if (useGPU) then
     ! Here we convert the regular host array into a pinned host array
     successCUDA = cuda_malloc(a_dev, lda*na_cols*size_of_real_datatype)
     if (.not.(successCUDA)) then
       print *,"bandred_real: error in cudaMalloc"
       stop
     endif

     successCUDA = cuda_malloc(tmat_dev, nbw*nbw*size_of_real_datatype)
     if (.not.(successCUDA)) then
       print *,"bandred_real: error in cudaMalloc"
       stop
     endif

     successCUDA = cuda_malloc(vav_dev, nbw*nbw*size_of_real_datatype)
     if (.not.(successCUDA)) then
       print *,"bandred_real: error in cudaMalloc"
       stop
     endif

     cur_l_rows = 0
     cur_l_cols = 0

     successCUDA = cuda_memcpy(a_dev, loc(a(1,1)), (lda)*(na_cols)*size_of_real_datatype,cudaMemcpyHostToDevice)
     if (.not.(successCUDA)) then
       print *,"bandred_real: error in cudaMemcpy"
       stop
     endif
   endif ! useGPU

   do istep = (na-1)/nbw, 1, -1

     n_cols = MIN(na,(istep+1)*nbw) - istep*nbw ! Number of columns in current step

     ! Number of local columns/rows of remaining matrix
     l_cols = local_index(istep*nbw, my_pcol, np_cols, nblk, -1)
     l_rows = local_index(istep*nbw, my_prow, np_rows, nblk, -1)

     if (useGPU) then
       cur_l_rows = max(l_rows, 1)
       cur_l_cols = max(l_cols, 1)

       vmr_size = cur_l_rows * 2 * n_cols
       umc_size = cur_l_cols * 2 * n_cols

       ! Allocate vmr and umc only if the inew size exceeds their current capacity
       ! Added for FORTRAN CALLS
       if ((.not. allocated(vr)) .or. (l_rows + 1 .gt. ubound(vr, dim=1))) then
         if (allocated(vr)) then
           deallocate(vr, stat=istat, errmsg=errorMessage)
           if (istat .ne. 0) then
             print *,"bandred_real: error when deallocating vr "//errorMessage
             stop
           endif
         endif
         allocate(vr(l_rows + 1), stat=istat, errmsg=errorMessage)
         if (istat .ne. 0) then
           print *,"bandred_real: error when allocating vr "//errorMessage
           stop
         endif

       endif

       if ((.not. allocated(vmrCUDA)) .or. (vmr_size .gt. ubound(vmrCUDA, dim=1))) then
         if (allocated(vmrCUDA)) then
           deallocate(vmrCUDA, stat=istat, errmsg=errorMessage)
           if (istat .ne. 0) then
             print *,"bandred_real: error when allocating vmrCUDA "//errorMessage
             stop
           endif

           successCUDA = cuda_free(vmr_dev)
           if (.not.(successCUDA)) then
             print *,"bandred_real: error in cuda_free"
             stop
           endif
         endif

         allocate(vmrCUDA(vmr_size), stat=istat, errmsg=errorMessage)
         if (istat .ne. 0) then
           print *,"bandred_real: error when allocating vmrCUDA "//errorMessage
           stop
         endif

         successCUDA = cuda_malloc(vmr_dev, vmr_size*size_of_real_datatype)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMalloc"
           stop
         endif

       endif

       if ((.not. allocated(umcCUDA)) .or. (umc_size .gt. ubound(umcCUDA, dim=1))) then
         if (allocated(umcCUDA)) then
           deallocate(umcCUDA, stat=istat, errmsg=errorMessage)
           if (istat .ne. 0) then
             print *,"bandred_real: error when deallocating umcCUDA "//errorMessage
             stop
           endif

           successCUDA = cuda_free(umc_dev)
           if (.not.(successCUDA)) then
              print *,"bandred_real: error in cudaFree"
              stop
           endif

         endif

         allocate(umcCUDA(umc_size), stat=istat, errmsg=errorMessage)
         if (istat .ne. 0) then
           print *,"bandred_real: error when deallocating umcCUDA "//errorMessage
           stop
         endif

         successCUDA = cuda_malloc(umc_dev, umc_size*size_of_real_datatype)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMalloc"
           stop
         endif

       endif
     else ! GPU not used
       ! Allocate vmr and umc to their exact sizes so that they can be used in bcasts and reduces

       allocate(vmrCPU(max(l_rows,1),2*n_cols), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating vmrCPU "//errorMessage
         stop
       endif

       allocate(umcCPU(max(l_cols,1),2*n_cols), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating umcCPU "//errorMessage
         stop
       endif

       allocate(vr(l_rows+1), stat=istat, errmsg=errorMessage)
       if (istat .ne. 0) then
         print *,"bandred_real: error when allocating vr "//errorMessage
         stop
       endif
     endif ! use GPU

     if (useGPU) then
       vmrCUDA(1 : cur_l_rows * n_cols) = 0.
     else
       vmrCPU(1:l_rows,1:n_cols) = 0.
     endif

     vr(:) = 0
     tmat(:,:,istep) = 0

     if (useGPU) then
       umcCUDA(1 : umc_size) = 0.

       lc_start = local_index(istep*nbw+1, my_pcol, np_cols, nblk, -1)
       lc_end   = local_index(istep*nbw+n_cols, my_pcol, np_cols, nblk, -1)
       lr_end   = local_index((istep-1)*nbw + n_cols, my_prow, np_rows, nblk, -1)

       if(lc_start .le. 0) lc_start = 1

       ! Here we assume that the processor grid and the block grid are aligned
       cur_pcol = pcol(istep*nbw+1, nblk, np_cols)

       if(my_pcol == cur_pcol) then

         successCUDA = cuda_memcpy2d(loc(a(1, lc_start)), lda*size_of_real_datatype,         &
                                    (a_dev + ((lc_start-1) * lda*size_of_real_datatype)),    &
                                    lda*size_of_real_datatype, lr_end*size_of_real_datatype, &
                                    (lc_end - lc_start+1), cudaMemcpyDeviceToHost)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy2d"
           stop
         endif

       endif
     endif ! useGPU

     ! Reduce current block to lower triangular form

     if (useQR) then
       if (which_qr_decomposition == 1) then
         call qr_pdgeqrf_2dcomm(a, lda, vmrCPU, max(l_rows,1), tauvector(1), &
                                  tmat(1,1,istep), nbw, work_blocked,       &
                                  work_size, na, n_cols, nblk, nblk,        &
                                  istep*nbw+n_cols-nbw, istep*nbw+n_cols, 1,&
                                  0, PQRPARAM, mpi_comm_rows, mpi_comm_cols,&
                                  blockheuristic)
       endif
     else

       do lc = n_cols, 1, -1

         ncol = istep*nbw + lc ! absolute column number of householder vector
         nrow = ncol - nbw ! Absolute number of pivot row

         lr  = local_index(nrow, my_prow, np_rows, nblk, -1) ! current row length
         lch = local_index(ncol, my_pcol, np_cols, nblk, -1) ! HV local column number

         tau = 0

         if (nrow == 1) exit ! Nothing to do

         cur_pcol = pcol(ncol, nblk, np_cols) ! Processor column owning current block

         if (my_pcol==cur_pcol) then

           ! Get vector to be transformed; distribute last element and norm of
           ! remaining elements to all procs in current column

           vr(1:lr) = a(1:lr,lch) ! vector to be transformed

           if (my_prow==prow(nrow, nblk, np_rows)) then
             aux1(1) = dot_product(vr(1:lr-1),vr(1:lr-1))
             aux1(2) = vr(lr)
           else
             aux1(1) = dot_product(vr(1:lr),vr(1:lr))
             aux1(2) = 0.
           endif

           call mpi_allreduce(aux1,aux2,2,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)

           vnorm2 = aux2(1)
           vrl    = aux2(2)

           ! Householder transformation

           call hh_transform_real(vrl, vnorm2, xf, tau)

           ! Scale vr and store Householder vector for back transformation

           vr(1:lr) = vr(1:lr) * xf
           if (my_prow==prow(nrow, nblk, np_rows)) then
             a(1:lr-1,lch) = vr(1:lr-1)
             a(lr,lch) = vrl
             vr(lr) = 1.
           else
             a(1:lr,lch) = vr(1:lr)
           endif

         endif

         ! Broadcast Householder vector and tau along columns

         vr(lr+1) = tau
         call MPI_Bcast(vr,lr+1,MPI_REAL8,cur_pcol,mpi_comm_cols,mpierr)

         if (useGPU) then
           vmrCUDA(cur_l_rows * (lc - 1) + 1 : cur_l_rows * (lc - 1) + lr) = vr(1:lr)
         else
           vmrCPU(1:lr,lc) = vr(1:lr)
         endif

         tau = vr(lr+1)
         tmat(lc,lc,istep) = tau ! Store tau in diagonal of tmat

         ! Transform remaining columns in current block with Householder vector
         ! Local dot product

         aux1 = 0

         nlc = 0 ! number of local columns
         do j=1,lc-1
           lcx = local_index(istep*nbw+j, my_pcol, np_cols, nblk, 0)
           if (lcx>0) then
             nlc = nlc+1
             if (lr>0) aux1(nlc) = dot_product(vr(1:lr),a(1:lr,lcx))
           endif
         enddo