elpa2_compute_real_template.X90

#if 0
!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), fomerly known as
!      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 Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!    This particular source code file contains additions, changes and
!    enhancements authored by Intel Corporation which is not part of
!    the ELPA consortium.
!
!    More information can be found here:
!    http://elpa.mpcdf.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".
#endif

    subroutine M_bandred_real_PRECISSION(na, a, a_dev, lda, nblk, nbw, matrixCols, numBlocks, mpi_comm_rows, mpi_comm_cols, &
                            tmat, tmat_dev, 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
      use elpa1_compute
#ifdef HAVE_DETAILED_TIMINGS
      use timings
#endif
#ifdef WITH_OPENMP
      use omp_lib
#endif
      use precision
      implicit none

      integer(kind=ik)                      :: na, lda, nblk, nbw, matrixCols, numBlocks, mpi_comm_rows, mpi_comm_cols
#ifdef USE_ASSUMED_SIZE
      real(kind=REAL_DATATYPE)              :: a(lda,*), tmat(nbw,nbw,*)
#else
      real(kind=REAL_DATATYPE)              :: a(lda,matrixCols), tmat(nbw,nbw,numBlocks)
#endif
      real(kind=REAL_DATATYPE)              :: eps
      logical, intent(in)                   :: useGPU

      integer(kind=ik)                      :: my_prow, my_pcol, np_rows, np_cols, mpierr
      integer(kind=ik)                      :: l_cols, l_rows, vmrCols
      integer(kind=ik)                      :: i, j, lcs, lce, lrs, lre, lc, lr, cur_pcol, n_cols, nrow
      integer(kind=ik)                      :: istep, ncol, lch, lcx, nlc, mynlc
      integer(kind=ik)                      :: tile_size, l_rows_tile, l_cols_tile

      real(kind=REAL_DATATYPE)              :: vnorm2, xf, aux1(nbw), aux2(nbw), vrl, tau, vav(nbw,nbw)

      real(kind=REAL_DATATYPE), allocatable :: tmpCUDA(:),  vmrCUDA(:),  umcCUDA(:)
      real(kind=REAL_DATATYPE), allocatable :: tmpCPU(:,:), vmrCPU(:,:), umcCPU(:,:)
      real(kind=REAL_DATATYPE), allocatable :: vr(:)
      ! needed for blocked QR decomposition
      integer(kind=ik)                      :: PQRPARAM(11), work_size
      real(kind=REAL_DATATYPE)              :: dwork_size(1)
      real(kind=REAL_DATATYPE), allocatable :: work_blocked(:), tauvector(:), blockheuristic(:)

      ! a_dev is passed from bandred_real to trans_ev_band
      integer(kind=C_intptr_T)              :: a_dev, vmr_dev, umc_dev, tmat_dev, vav_dev
#ifdef WITH_MPI
      integer(kind=ik), external            :: numroc
#endif
      integer(kind=ik)                      :: ierr
      integer(kind=ik)                      :: cur_l_rows, cur_l_cols, vmr_size, umc_size
      integer(kind=c_size_t)                :: lc_start, lc_end
      integer(kind=ik)                      :: lr_end
      integer(kind=ik)                      :: na_cols !, na_rows

      logical, intent(in)                   :: wantDebug
      logical, intent(out)                  :: success
      logical                               :: successCUDA
      integer(kind=ik)                      :: istat
      character(200)                        :: errorMessage

      logical, intent(in)                   :: useQR

      integer(kind=ik)                      :: mystart, myend, m_way, n_way, work_per_thread, m_id, n_id, n_threads, &
                                               ii, pp, transformChunkSize
#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("bandred_real" + M_PRECISSION_SUFFIX)
#endif
#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("mpi_communication")
#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)
#ifdef HAVE_DETAILED_TIMINGS
      call timer%stop("mpi_communication")
#endif
      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

! na_rows in used nowhere; only na_cols
      if (useGPU) then
#ifdef WITH_MPI
!        na_rows = numroc(na, nblk, my_prow, 0, np_rows)
        na_cols = numroc(na, nblk, my_pcol, 0, np_cols)
#else
!        na_rows = na
        na_cols = na
#endif
      endif ! useGPU

      ! 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(1:11),    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

          vmrCols = na

#ifdef USE_ASSUMED_SIZE_QR
          call M_qr_pdgeqrf_2dcomm_PRECISSION(a, lda, matrixCols, vmrCPU, max(l_rows,1), vmrCols, tauvector(1), na, tmat(1,1,1), &
                                 nbw, nbw, dwork_size, 1, -1, na, nbw, nblk, nblk, na, na, 1, 0, PQRPARAM(1:11), &
                                 mpi_comm_rows, mpi_comm_cols, blockheuristic)

#else
          call M_qr_pdgeqrf_2dcomm_PRECISSION(a(1:lda,1:matrixCols), matrixCols, lda, vmrCPU(1:max(l_rows,1),1:vmrCols), max(l_rows,1), &
                                 vmrCols, tauvector(1:na), na, tmat(1:nbw,1:nbw,1), nbw, &
                                 nbw, dwork_size(1:1), 1, -1, na, nbw, nblk, nblk, na, na, 1, 0, PQRPARAM(1:11), &
                                 mpi_comm_rows, mpi_comm_cols, blockheuristic)
#endif

          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 = M_CONST_0_0
          deallocate(vmrCPU, stat=istat, errmsg=errorMessage)
          if (istat .ne. 0) then
            print *,"bandred_real: error when deallocating vmrCPU "//errorMessage
            stop
          endif

        endif ! which_qr_decomposition

      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*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"bandred_real: error in cudaMalloc"
          stop
        endif
        successCUDA = cuda_malloc(tmat_dev, nbw*nbw*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"bandred_real: error in cudaMalloc"
          stop
        endif
        successCUDA = cuda_malloc(vav_dev, nbw*nbw*M_size_of_PRECISSION_real_datatape)
        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)*M_size_of_PRECISSION_real_datatape, 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*M_size_of_PRECISSION_real_datatape)
            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*M_size_of_PRECISSION_real_datatape)
            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) = M_CONST_0_0
        else
          vmrCPU(1:l_rows,1:n_cols) = M_CONST_0_0
        endif
        vr(:) = M_CONST_0_0
        tmat(:,:,istep) = M_CONST_0_0

        if (useGPU) then
          umcCUDA(1 : umc_size) = M_CONST_0_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*M_size_of_PRECISSION_real_datatape,         &
                                       (a_dev + ((lc_start-1) * lda*M_size_of_PRECISSION_real_datatape)),    &
                                       lda*M_size_of_PRECISSION_real_datatape, lr_end*M_size_of_PRECISSION_real_datatape, &
                                       (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
            vmrCols = 2*n_cols
#ifdef USE_ASSUMED_SIZE_QR
            call M_qr_pdgeqrf_2dcomm_PRECISSION(a, lda, matrixCols, vmrCPU, max(l_rows,1), vmrCols, tauvector(1), &
                                   na, tmat(1,1,istep), nbw, nbw, work_blocked, work_size,        &
                                     work_size, na, n_cols, nblk, nblk,        &
                                     istep*nbw+n_cols-nbw, istep*nbw+n_cols, 1,&
                                     0, PQRPARAM(1:11), mpi_comm_rows, mpi_comm_cols,&
                                     blockheuristic)

#else
            call M_qr_pdgeqrf_2dcomm_PRECISSION(a(1:lda,1:matrixCols), lda, matrixCols, vmrCPU(1:max(l_rows,1),1:vmrCols) ,   &
                                    max(l_rows,1), vmrCols, tauvector(1:na), na, &
                                     tmat(1:nbw,1:nbw,istep), nbw, nbw, work_blocked(1:work_size), work_size, &
                                     work_size, na, n_cols, nblk, nblk,        &
                                     istep*nbw+n_cols-nbw, istep*nbw+n_cols, 1,&
                                     0, PQRPARAM(1:11), mpi_comm_rows, mpi_comm_cols,&
                                     blockheuristic)
#endif
          endif

       else !useQR

         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) = M_CONST_0_0
             endif

#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
             call timer%start("mpi_communication")
#endif
             call mpi_allreduce(aux1, aux2, 2, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
             call timer%stop("mpi_communication")
#endif

#else /* WITH_MPI */
              aux2 = aux1 ! this should be optimized
#endif 

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

             ! Householder transformation
             call M_hh_transform_real_PRECISSION(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) = M_CONST_1_0
             else
               a(1:lr,lch) = vr(1:lr)
             endif

           endif

           ! Broadcast Householder vector and tau along columns

           vr(lr+1) = tau
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
           call timer%start("mpi_communication")
#endif
           call MPI_Bcast(vr, lr+1, M_MPI_REAL_PRECISSION, cur_pcol, mpi_comm_cols, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
           call timer%stop("mpi_communication")
#endif

#endif /* WITH_MPI */
           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

#ifdef WITH_OPENMP
           !Open up one omp region to avoid paying openmp overhead.
           !This does not help performance due to the addition of two openmp barriers around the MPI call,
           !But in the future this may be beneficial if these barriers are replaced with a faster implementation

           !$omp parallel private(mynlc, j, lcx, ii, pp ) shared(aux1)
           mynlc = 0 ! number of local columns

           !This loop does not have independent iterations,
           !'mynlc' is incremented each iteration, and it is difficult to remove this dependency
           !Thus each thread executes every iteration of the loop, except it only does the work if it 'owns' that iteration
           !That is, a thread only executes the work associated with an iteration if its thread id is congruent to
           !the iteration number modulo the number of threads
           do j=1,lc-1
             lcx = local_index(istep*nbw+j, my_pcol, np_cols, nblk, 0)
             if (lcx>0 ) then
               mynlc = mynlc+1
               if ( mod((j-1), omp_get_num_threads()) .eq. omp_get_thread_num() ) then
                   if (lr>0) aux1(mynlc) = dot_product(vr(1:lr),a(1:lr,lcx))
               endif
             endif
           enddo

           ! Get global dot products

           !$omp barrier
           !$omp single
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
           call timer%start("mpi_communication")
#endif
           if (mynlc>0) call mpi_allreduce(aux1, aux2, mynlc, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
           call timer%stop("mpi_communication")
#endif
#else /* WITH_MPI */
           if (mynlc>0) aux2 = aux1
#endif /* WITH_MPI */
           !$omp end single
           !$omp barrier

           ! Transform
           transformChunkSize=32
           mynlc = 0
           do j=1,lc-1
             lcx = local_index(istep*nbw+j, my_pcol, np_cols, nblk, 0)
             if (lcx>0) then
               mynlc = mynlc+1
               !This loop could be parallelized with an openmp pragma with static scheduling and chunk size 32
               !However, for some reason this is slower than doing it manually, so it is parallelized as below.
               do ii=omp_get_thread_num()*transformChunkSize,lr,omp_get_num_threads()*transformChunkSize
                  do pp = 1,transformChunkSize
                      if (pp + ii > lr) exit
                          a(ii+pp,lcx) = a(ii+pp,lcx) - tau*aux2(mynlc)*vr(ii+pp)
                  enddo
               enddo
             endif
           enddo
           !$omp end parallel

#else /* WITH_OPENMP */

           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

           ! Get global dot products
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
           call timer%start("mpi_communication")
#endif
           if (nlc>0) call mpi_allreduce(aux1, aux2, nlc, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
           call timer%stop("mpi_communication")
#endif
#else /* WITH_MPI */
           if (nlc>0) aux2=aux1
#endif /* WITH_MPI */
           ! Transform

           nlc = 0
           do j=1,lc-1
             lcx = local_index(istep*nbw+j, my_pcol, np_cols, nblk, 0)
             if (lcx>0) then
               nlc = nlc+1
               a(1:lr,lcx) = a(1:lr,lcx) - tau*aux2(nlc)*vr(1:lr)
             endif
           enddo
#endif /* WITH_OPENMP */
         enddo ! lc

         if (useGPU) then
           ! store column tiles back to GPU
           cur_pcol = pcol(istep*nbw+1, nblk, np_cols)
           if (my_pcol == cur_pcol) then
             successCUDA = cuda_memcpy2d((a_dev+((lc_start-1)*lda*M_size_of_PRECISSION_real_datatape)),          &
                                          lda*M_size_of_PRECISSION_real_datatape, loc(a(1, lc_start)),           &
                                          lda*M_size_of_PRECISSION_real_datatape,  lr_end*M_size_of_PRECISSION_real_datatape, &
                                          (lc_end - lc_start+1),cudaMemcpyHostToDevice)
             if (.not.(successCUDA)) then
               print *,"bandred_real: error in cudaMemcpy2d"
               stop
             endif

           endif
         endif

         ! Calculate scalar products of stored Householder vectors.
         ! This can be done in different ways, we use dsyrk

         vav = 0

         if (useGPU) then
           if (l_rows>0) &
             call M_PRECISSION_SYRK('U', 'T', n_cols, l_rows, M_CONST_1_0, vmrCUDA, cur_l_rows, M_CONST_0_0, vav, ubound(vav,dim=1))
         else
           if (l_rows>0) &
             call M_PRECISSION_SYRK('U', 'T', n_cols, l_rows, M_CONST_1_0, vmrCPU, ubound(vmrCPU,dim=1), M_CONST_0_0, vav, ubound(vav,dim=1))
         endif

         call M_symm_matrix_allreduce_PRECISSION(n_cols,vav, nbw, nbw,mpi_comm_rows)

         ! Calculate triangular matrix T for block Householder Transformation
         do lc=n_cols,1,-1
           tau = tmat(lc,lc,istep)
           if (lc<n_cols) then
             call M_PRECISSION_TRMV('U', 'T', 'N', n_cols-lc, tmat(lc+1,lc+1,istep), ubound(tmat,dim=1), vav(lc+1,lc), 1)
             tmat(lc,lc+1:n_cols,istep) = -tau * vav(lc+1:n_cols,lc)
           endif
         enddo
       endif

       ! Transpose vmr -> vmc (stored in umc, second half)

       if (useGPU) then
         call M_elpa_transpose_vectors_real_PRECISSION  (vmrCUDA, cur_l_rows, mpi_comm_rows, &
                                            umcCUDA(cur_l_cols * n_cols + 1), cur_l_cols, mpi_comm_cols, &
                                            1, istep*nbw, n_cols, nblk)
       else
         call M_elpa_transpose_vectors_real_PRECISSION  (vmrCPU, ubound(vmrCPU,dim=1), mpi_comm_rows, &
                                            umcCPU(1,n_cols+1), ubound(umcCPU,dim=1), mpi_comm_cols, &
                                            1, istep*nbw, n_cols, nblk)
       endif

       ! Calculate umc = A**T * vmr
       ! Note that the distributed A has to be transposed
       ! Opposed to direct tridiagonalization there is no need to use the cache locality
       ! of the tiles, so we can use strips of the matrix

       ! here the GPU version and CPU version diverged substantially, due to the newest
       ! optimizations due to Intel. The GPU version has to be re-written
       if (useGPU) then
         umcCUDA(1 : l_cols * n_cols) = M_CONST_0_0
         vmrCUDA(cur_l_rows * n_cols + 1 : cur_l_rows * n_cols * 2) = M_CONST_0_0

         if (l_cols>0 .and. l_rows>0) then
           successCUDA = cuda_memcpy(vmr_dev, loc(vmrCUDA(1)), vmr_size*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
           if (.not.(successCUDA)) then
             print *,"bandred_real: error in cudaMemcpy"
             stop
           endif
           successCUDA = cuda_memcpy(umc_dev, loc(umcCUDA(1)), umc_size*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
           if (.not.(successCUDA)) then
             print *,"bandred_real: error in cudaMemcpy"
             stop
           endif

           do i=0,(istep*nbw-1)/tile_size

             lcs = i*l_cols_tile+1
             lce = min(l_cols,(i+1)*l_cols_tile)
             if (lce<lcs) cycle

             lre = min(l_rows,(i+1)*l_rows_tile)
             call M_cublas_PRECISSION_gemm('T', 'N', lce-lcs+1, n_cols, lre, &
                               M_CONST_1_0, (a_dev + ((lcs-1)*lda*M_size_of_PRECISSION_real_datatape)), lda, vmr_dev,cur_l_rows, &
                               M_CONST_1_0, (umc_dev+ (lcs-1)*M_size_of_PRECISSION_real_datatape), cur_l_cols)
             if(i==0) cycle
             lre = min(l_rows,i*l_rows_tile)
             call M_cublas_PRECISSION_gemm('N', 'N', lre,n_cols, lce-lcs+1,&
                               M_CONST_1_0, (a_dev+ ((lcs-1)*lda*M_size_of_PRECISSION_real_datatape)), lda,                  &
                               (umc_dev+(cur_l_cols * n_cols+lcs-1)*M_size_of_PRECISSION_real_datatape), cur_l_cols, &
                               M_CONST_1_0, (vmr_dev+(cur_l_rows * n_cols)*M_size_of_PRECISSION_real_datatape), cur_l_rows)
           enddo
           successCUDA = cuda_memcpy(loc(vmrCUDA(1)), vmr_dev,vmr_size*M_size_of_PRECISSION_real_datatape,cudaMemcpyDeviceToHost)
           if (.not.(successCUDA)) then
             print *,"bandred_real: error in cudaMemcpy"
             stop
           endif
           successCUDA = cuda_memcpy(loc(umcCUDA(1)), umc_dev, umc_size*M_size_of_PRECISSION_real_datatape,cudaMemcpyDeviceToHost)
           if (.not.(successCUDA)) then
             print *,"bandred_real: error in cudaMemcpy"
             stop
           endif

         endif ! l_cols>0 .and. l_rows>0

       else ! do not useGPU version

         !Code for Algorithm 4

         n_way = 1
#ifdef WITH_OPENMP
         n_way = omp_get_max_threads()
#endif
         !umc(1:l_cols,1:n_cols) = 0.d0
         !vmr(1:l_rows,n_cols+1:2*n_cols) = 0
#ifdef WITH_OPENMP
         !$omp parallel private( i,lcs,lce,lrs,lre)
#endif
         if (n_way > 1) then
           !$omp do
           do i=1,min(l_cols_tile, l_cols)
             umcCPU(i,1:n_cols) = M_CONST_0_0
           enddo

           !$omp do
           do i=1,l_rows
             vmrCPU(i,n_cols+1:2*n_cols) = M_CONST_0_0
           enddo
           if (l_cols>0 .and. l_rows>0) then

             !SYMM variant 4
             !Partitioned Matrix Expression:
             ! Ct = Atl Bt + Atr Bb
             ! Cb = Atr' Bt + Abl Bb
             !
             !Loop invariant:
             ! Ct = Atl Bt + Atr Bb
             !
             !Update:
             ! C1 = A10'B0 + A11B1 + A21 B2
             !
             !This algorithm chosen because in this algoirhtm, the loop around the dgemm calls
             !is easily parallelized, and regardless of choise of algorithm,
             !the startup cost for parallelizing the dgemms inside the loop is too great

             !$omp do schedule(static,1)
             do i=0,(istep*nbw-1)/tile_size
               lcs = i*l_cols_tile+1                   ! local column start
               lce = min(l_cols, (i+1)*l_cols_tile)    ! local column end

               lrs = i*l_rows_tile+1                   ! local row start
               lre = min(l_rows, (i+1)*l_rows_tile)    ! local row end

               !C1 += [A11 A12] [B1
               !                 B2]
               if ( lre > lrs .and. l_cols > lcs ) then
                 call M_PRECISSION_GEMM('N', 'N', lre-lrs+1, n_cols, l_cols-lcs+1,          &
                            M_CONST_1_0, a(lrs,lcs), ubound(a,dim=1),                 &
                                  umcCPU(lcs,n_cols+1), ubound(umcCPU,dim=1),  &
                            M_CONST_0_0, vmrCPU(lrs,n_cols+1), ubound(vmrCPU,dim=1))
               endif

               ! C1 += A10' B0
               if ( lce > lcs .and. i > 0 ) then
                 call M_PRECISSION_GEMM('T', 'N', lce-lcs+1, n_cols, lrs-1,           &
                            M_CONST_1_0, a(1,lcs),   ubound(a,dim=1),           &
                                  vmrCPU(1,1),   ubound(vmrCPU,dim=1),   &
                            M_CONST_0_0, umcCPU(lcs,1), ubound(umcCPU,dim=1))
               endif
             enddo
           endif ! l_cols>0 .and. l_rows>0
         else ! n_way > 1
           umcCPU(1:l_cols,1:n_cols) = M_CONST_0_0
           vmrCPU(1:l_rows,n_cols+1:2*n_cols) = M_CONST_0_0
           if (l_cols>0 .and. l_rows>0) then
             do i=0,(istep*nbw-1)/tile_size

               lcs = i*l_cols_tile+1
               lce = min(l_cols,(i+1)*l_cols_tile)
               if (lce<lcs) cycle

               lre = min(l_rows,(i+1)*l_rows_tile)
               call M_PRECISSION_GEMM('T', 'N', lce-lcs+1, n_cols, lre, M_CONST_1_0, a(1,lcs), ubound(a,dim=1), &
                            vmrCPU, ubound(vmrCPU,dim=1), M_CONST_1_0, umcCPU(lcs,1), ubound(umcCPU,dim=1))

               if (i==0) cycle
                 lre = min(l_rows,i*l_rows_tile)
                 call M_PRECISSION_GEMM('N', 'N', lre, n_cols, lce-lcs+1, M_CONST_1_0, a(1,lcs), lda, &
                            umcCPU(lcs,n_cols+1), ubound(umcCPU,dim=1), M_CONST_1_0, vmrCPU(1,n_cols+1), ubound(vmrCPU,dim=1))
             enddo
           endif
         endif ! n_way > 1
#ifdef WITH_OPENMP
        !$omp end parallel
#endif
       endif ! do not useGPU version

       ! Sum up all ur(:) parts along rows and add them to the uc(:) parts
       ! on the processors containing the diagonal
       ! This is only necessary if ur has been calculated, i.e. if the
       ! global tile size is smaller than the global remaining matrix

       if (useGPU) then
         ! here the GPU version and CPU version divereged due to the same reasons as above

         if (tile_size < istep*nbw) then
           call M_elpa_reduce_add_vectors_real_PRECISSION  (vmrCUDA(cur_l_rows * n_cols + 1),cur_l_rows,mpi_comm_rows, &
                                               umcCUDA, cur_l_cols, mpi_comm_cols, &
                                               istep*nbw, n_cols, nblk)
         endif

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

#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
           call timer%start("mpi_communication")
#endif

           call mpi_allreduce(umcCUDA, tmpCUDA, l_cols*n_cols, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, ierr)
           umcCUDA(1 : l_cols * n_cols) = tmpCUDA(1 : l_cols * n_cols)
#ifdef HAVE_DETAILED_TIMINGS
           call timer%stop("mpi_communication")
#endif
#else /* WITH_MPI */

!           tmpCUDA(1 : l_cols * n_cols) = umcCUDA(1 : l_cols * n_cols)

#endif /* WITH_MPI */

           if (allocated(tmpCUDA)) then
             deallocate(tmpCUDA, stat=istat, errmsg=errorMessage)
             if (istat .ne. 0) then
               print *,"bandred_real: error when deallocating tmpCUDA "//errorMessage
               stop
             endif
           endif
         endif ! l_cols

         ! U = U * Tmat**T
         successCUDA = cuda_memcpy(umc_dev, loc(umcCUDA(1)), umc_size*M_size_of_PRECISSION_real_datatape, cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif
         successCUDA = cuda_memcpy(tmat_dev,loc(tmat(1,1,istep)),nbw*nbw*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif
         call M_cublas_PRECISSION_trmm('Right', 'Upper', 'Trans', 'Nonunit', l_cols, n_cols, &
                           M_CONST_1_0, tmat_dev, nbw, umc_dev, cur_l_cols)
         ! VAV = Tmat * V**T * A * V * Tmat**T = (U*Tmat**T)**T * V * Tmat**T
         successCUDA = cuda_memcpy(vav_dev,loc(vav(1,1)), nbw*nbw*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif
         call M_cublas_PRECISSION_gemm('T', 'N', n_cols, n_cols, l_cols, &
                           M_CONST_1_0, umc_dev, cur_l_cols, (umc_dev+(cur_l_cols * n_cols )*M_size_of_PRECISSION_real_datatape),cur_l_cols, &
                           M_CONST_0_0, vav_dev, nbw)

         call M_cublas_PRECISSION_trmm('Right', 'Upper', 'Trans', 'Nonunit', n_cols, n_cols, &
                           M_CONST_1_0, tmat_dev, nbw, vav_dev, nbw)

         successCUDA = cuda_memcpy(loc(vav(1,1)), vav_dev, nbw*nbw*M_size_of_PRECISSION_real_datatape, cudaMemcpyDeviceToHost)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif

         call M_symm_matrix_allreduce_PRECISSION(n_cols,vav, nbw,nbw,mpi_comm_cols)

         successCUDA = cuda_memcpy(vav_dev, loc(vav(1,1)), nbw*nbw*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif

         ! U = U - 0.5 * V * VAV
         call M_cublas_PRECISSION_gemm('N', 'N', l_cols, n_cols, n_cols,&
                           -M_CONST_0_5, (umc_dev+(cur_l_cols * n_cols )*M_size_of_PRECISSION_real_datatape),cur_l_cols, vav_dev,nbw,&
                           M_CONST_1_0, umc_dev, cur_l_cols)

         successCUDA = cuda_memcpy(loc(umcCUDA(1)), umc_dev, umc_size*M_size_of_PRECISSION_real_datatape, cudaMemcpyDeviceToHost)

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

         ! Transpose umc -> umr (stored in vmr, second half)
         call M_elpa_transpose_vectors_real_PRECISSION  (umcCUDA, cur_l_cols, mpi_comm_cols, &
                                            vmrCUDA(cur_l_rows * n_cols + 1), cur_l_rows, mpi_comm_rows, &
                                            1, istep*nbw, n_cols, nblk)

         successCUDA = cuda_memcpy(vmr_dev, loc(vmrCUDA(1)), vmr_size*M_size_of_PRECISSION_real_datatape, cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif

         successCUDA = cuda_memcpy(umc_dev, loc(umcCUDA(1)), umc_size*M_size_of_PRECISSION_real_datatape, cudaMemcpyHostToDevice)
         if (.not.(successCUDA)) then
           print *,"bandred_real: error in cudaMemcpy"
           stop
         endif

         ! A = A - V*U**T - U*V**T
         do i=0,(istep*nbw-1)/tile_size
           lcs = i*l_cols_tile+1
           lce = min(l_cols,(i+1)*l_cols_tile)
           lre = min(l_rows,(i+1)*l_rows_tile)
           if (lce<lcs .or. lre<1) cycle
           call M_cublas_PRECISSION_gemm('N', 'T', lre, lce-lcs+1, 2*n_cols, -M_CONST_1_0, &
                             vmr_dev, cur_l_rows, (umc_dev +(lcs-1)*M_size_of_PRECISSION_real_datatape), cur_l_cols, &
                             M_CONST_1_0, (a_dev+(lcs-1)*lda*M_size_of_PRECISSION_real_datatape), lda)
         enddo

       else ! do not useGPU

         ! Or if we used the Algorithm 4
         if (tile_size < istep*nbw .or. n_way > 1) then
           call M_elpa_reduce_add_vectors_real_PRECISSION  (vmrCPU(1,n_cols+1),ubound(vmrCPU,dim=1),mpi_comm_rows, &
                                             umcCPU, ubound(umcCPU,dim=1), mpi_comm_cols, &
                                             istep*nbw, n_cols, nblk)
         endif

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

#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
           call timer%start("mpi_communication")
#endif
           call mpi_allreduce(umcCPU, tmpCPU, l_cols*n_cols, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
           umcCPU(1:l_cols,1:n_cols) = tmpCPU(1:l_cols,1:n_cols)
#ifdef HAVE_DETAILED_TIMINGS
           call timer%stop("mpi_communication")
#endif
#else /* WITH_MPI */
!           tmpCPU(1:l_cols,1:n_cols) = umcCPU(1:l_cols,1:n_cols)
#endif /* WITH_MPI */

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

         ! U = U * Tmat**T
         call M_PRECISSION_TRMM('Right', 'Upper', 'Trans', 'Nonunit', l_cols,n_cols, M_CONST_1_0, tmat(1,1,istep), ubound(tmat,dim=1), &
                    umcCPU, ubound(umcCPU,dim=1))

         ! VAV = Tmat * V**T * A * V * Tmat**T = (U*Tmat**T)**T * V * Tmat**T

         call M_PRECISSION_GEMM('T', 'N', n_cols, n_cols, l_cols, M_CONST_1_0, umcCPU, ubound(umcCPU,dim=1), umcCPU(1,n_cols+1), &
                    ubound(umcCPU,dim=1), M_CONST_0_0, vav, ubound(vav,dim=1))

         call M_PRECISSION_TRMM('Right', 'Upper', 'Trans', 'Nonunit', n_cols, n_cols, M_CONST_1_0, tmat(1,1,istep),    &
                    ubound(tmat,dim=1), vav, ubound(vav,dim=1))

         call M_symm_matrix_allreduce_PRECISSION(n_cols,vav, nbw, nbw ,mpi_comm_cols)

         ! U = U - 0.5 * V * VAV
         call M_PRECISSION_GEMM('N', 'N', l_cols, n_cols, n_cols, -M_CONST_0_5, umcCPU(1,n_cols+1), ubound(umcCPU,dim=1), vav, &
                     ubound(vav,dim=1), M_CONST_1_0, umcCPU, ubound(umcCPU,dim=1))

         ! Transpose umc -> umr (stored in vmr, second half)
         call M_elpa_transpose_vectors_real_PRECISSION(umcCPU, ubound(umcCPU,dim=1), mpi_comm_cols, &
                                         vmrCPU(1,n_cols+1), ubound(vmrCPU,dim=1), mpi_comm_rows, &
                                         1, istep*nbw, n_cols, nblk)

         ! A = A - V*U**T - U*V**T

#ifdef WITH_OPENMP
         !$omp parallel private( ii, i, lcs, lce, lre, n_way, m_way, m_id, n_id, work_per_thread, mystart, myend  )
         n_threads = omp_get_num_threads()
         if (mod(n_threads, 2) == 0) then
           n_way = 2
         else
           n_way = 1
         endif

         m_way = n_threads / n_way

         m_id = mod(omp_get_thread_num(),  m_way)
         n_id = omp_get_thread_num() / m_way

         do ii=n_id*tile_size,(istep*nbw-1),tile_size*n_way
           i = ii / tile_size
           lcs = i*l_cols_tile+1
           lce = min(l_cols,(i+1)*l_cols_tile)
           lre = min(l_rows,(i+1)*l_rows_tile)
           if (lce<lcs .or. lre<1) cycle

           !Figure out this thread's range
           work_per_thread = lre / m_way
           if (work_per_thread * m_way < lre) work_per_thread = work_per_thread + 1
           mystart = m_id * work_per_thread + 1
           myend   = mystart + work_per_thread - 1
           if ( myend > lre ) myend = lre
           if ( myend-mystart+1 < 1) cycle
           call M_PRECISSION_GEMM('N', 'T', myend-mystart+1, lce-lcs+1, 2*n_cols, -M_CONST_1_0, &
                      vmrCPU(mystart, 1), ubound(vmrCPU,1), umcCPU(lcs,1), ubound(umcCPU,1), &
                       M_CONST_1_0, a(mystart,lcs), ubound(a,1))
         enddo
         !$omp end parallel

#else /* WITH_OPENMP */

         do i=0,(istep*nbw-1)/tile_size
           lcs = i*l_cols_tile+1
           lce = min(l_cols,(i+1)*l_cols_tile)
           lre = min(l_rows,(i+1)*l_rows_tile)
           if (lce<lcs .or. lre<1) cycle
           call M_PRECISSION_GEMM('N', 'T', lre,lce-lcs+1, 2*n_cols, -M_CONST_1_0, &
                       vmrCPU, ubound(vmrCPU,dim=1), umcCPU(lcs,1), ubound(umcCPU,dim=1), &
                       M_CONST_1_0, a(1,lcs), lda)
         enddo
#endif /* WITH_OPENMP */

       endif ! useGPU

       if (.not.(useGPU)) 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

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

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

       endif !useGPU

     enddo ! istep

     if (useGPU) then
       ! this is not needed since a_dev is passed along from one subroutine to the other
       successCUDA = cuda_memcpy ( loc (a), a_dev, lda*na_cols*M_size_of_PRECISSION_real_datatape,cudaMemcpyDeviceToHost)
       if (.not.(successCUDA)) then
         print *,"bandred_real: error in cudaMemcpy"
         stop
       endif

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

!#ifdef WITH_MPI
! it should be possible to keep tmat dev on the device and not copy it arround
! this is not necessary tmat_dev is passed (unchanged) from one routine to the other
       successCUDA = cuda_free(tmat_dev)
       if (.not.(successCUDA)) then
         print *,"bandred_real: error in cudaFree"
         stop
       endif
!#endif

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

     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

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

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

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

       successCUDA = cuda_free(umc_dev)
       if (.not.(successCUDA)) then
         print *,"bandred_real: error in cudaFree"
         stop
       endif
       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
       endif
       if (allocated(vmrCUDA)) then
         deallocate(vmrCUDA, stat=istat, errmsg=errorMessage)
         if (istat .ne. 0) then
           print *,"bandred_real: error when deallocating vmrCUDA "//errorMessage
           stop
         endif
       endif

     endif ! useGPU

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

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

#ifdef HAVE_DETAILED_TIMINGS
     call timer%stop("bandred_real" + M_PRECISSION_SUFFIX)
#endif

   end subroutine M_bandred_real_PRECISSION ! slower for gpu on 10000 10000 ???

    subroutine M_symm_matrix_allreduce_PRECISSION(n,a,lda,ldb,comm)
    !-------------------------------------------------------------------------------
    !  symm_matrix_allreduce: Does an mpi_allreduce for a symmetric matrix A.
    !  On entry, only the upper half of A needs to be set
    !  On exit, the complete matrix is set
    !-------------------------------------------------------------------------------
#ifdef HAVE_DETAILED_TIMINGS
      use timings
#endif
      use precision
      implicit none
      integer(kind=ik)             :: n, lda, ldb, comm
#ifdef USE_ASSUMED_SIZE
      real(kind=REAL_DATATYPE)     :: a(lda,*)
#else
      real(kind=REAL_DATATYPE)     :: a(lda,ldb)
#endif
      integer(kind=ik)             :: i, nc, mpierr
      real(kind=REAL_DATATYPE)     :: h1(n*n), h2(n*n)

#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("symm_matrix_allreduce" + M_PRECISSION_SUFFIX)
#endif

      nc = 0
      do i=1,n
        h1(nc+1:nc+i) = a(1:i,i)
        nc = nc+i
      enddo

#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("mpi_communication")
#endif
      call mpi_allreduce(h1, h2, nc, M_MPI_REAL_PRECISSION, MPI_SUM, comm, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
      call timer%stop("mpi_communication")
#endif
      nc = 0
      do i=1,n
        a(1:i,i) = h2(nc+1:nc+i)
        a(i,1:i-1) = a(1:i-1,i)
        nc = nc+i
      enddo

#else /* WITH_MPI */
!      h2=h1

      nc = 0
      do i=1,n
        a(1:i,i) = h1(nc+1:nc+i)
        a(i,1:i-1) = a(1:i-1,i)
        nc = nc+i
      enddo

#endif /* WITH_MPI */
!      nc = 0
!      do i=1,n
!        a(1:i,i) = h2(nc+1:nc+i)
!        a(i,1:i-1) = a(1:i-1,i)
!        nc = nc+i
!      enddo

#ifdef HAVE_DETAILED_TIMINGS
      call timer%stop("symm_matrix_allreduce" + M_PRECISSION_SUFFIX)
#endif

    end subroutine M_symm_matrix_allreduce_PRECISSION

    subroutine M_trans_ev_band_to_full_real_PRECISSION(na, nqc, nblk, nbw, a, a_dev, lda, tmat, tmat_dev, q, q_dev, ldq, matrixCols, &
                                                       numBlocks, mpi_comm_rows, mpi_comm_cols, useGPU, useQR)
    !-------------------------------------------------------------------------------
    !  trans_ev_band_to_full_real:
    !  Transforms the eigenvectors of a band matrix back to the eigenvectors of the original matrix
    !
    !  Parameters
    !
    !  na          Order of matrix a, number of rows of matrix q
    !
    !  nqc         Number of columns of matrix q
    !
    !  nblk        blocksize of cyclic distribution, must be the same in both directions!
    !
    !  nbw         semi bandwith
    !
    !  a(lda,matrixCols)    Matrix containing the Householder vectors (i.e. matrix a after bandred_real)
    !              Distribution is like in Scalapack.
    !
    !  lda         Leading dimension of a
    !  matrixCols  local columns of matrix a and q
    !
    !  tmat(nbw,nbw,numBlocks) Factors returned by bandred_real
    !
    !  q           On input: Eigenvectors of band matrix
    !              On output: Transformed eigenvectors
    !              Distribution is like in Scalapack.
    !
    !  ldq         Leading dimension of q
    !
    !  mpi_comm_rows
    !  mpi_comm_cols
    !              MPI-Communicators for rows/columns
    !
    !-------------------------------------------------------------------------------
#ifdef HAVE_DETAILED_TIMINGS
      use timings
#endif
      use precision
      use cuda_functions
      use iso_c_binding

      implicit none

      integer(kind=ik)                       :: na, nqc, lda, ldq, nblk, nbw, matrixCols, numBlocks, mpi_comm_rows, mpi_comm_cols
#ifdef USE_ASSUMED_SIZE
      real(kind=REAL_DATATYPE)               :: a(lda,*), q(ldq,*), tmat(nbw,nbw,*)
#else
      real(kind=REAL_DATATYPE)               :: a(lda,matrixCols), q(ldq,matrixCols), tmat(nbw, nbw, numBlocks)
#endif

      integer(kind=C_intptr_T)               :: a_dev ! passed from bandred_real at the moment not used since copied in bandred_real

      integer(kind=ik)                       :: my_prow, my_pcol, np_rows, np_cols, mpierr
      integer(kind=ik)                       :: max_blocks_row, max_blocks_col, max_local_rows, &
                                                max_local_cols
      integer(kind=ik)                       :: l_cols, l_rows, l_colh, n_cols
      integer(kind=ik)                       :: istep, lc, ncol, nrow, nb, ns

      real(kind=REAL_DATATYPE), allocatable  :: tmp1(:), tmp2(:), hvb(:), hvm(:,:)

      ! hvm_dev is fist used and set in this routine
      ! q is changed in trans_ev_tridi on the host, copied to device and passed here. this can be adapted
      ! tmp_dev is first used in this routine
      ! tmat_dev is passed along from bandred_real
      integer(kind=C_intptr_T)               :: hvm_dev, q_dev, tmp_dev, tmat_dev

      integer(kind=ik)                       :: i

      real(kind=REAL_DATATYPE), allocatable  :: tmat_complete(:,:), t_tmp(:,:), t_tmp2(:,:)
      integer(kind=ik)                       :: cwy_blocking, t_blocking, t_cols, t_rows
      logical, intent(in)                    :: useQR, useGPU
      integer(kind=ik)                       :: istat
      character(200)                         :: errorMessage
      logical                                :: successCUDA

#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("trans_ev_band_to_full_real" + M_PRECISSION_SUFFIX)
#endif
#ifdef HAVE_DETAILED_TIMINGS
      call timer%start("mpi_communication")
#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)
#ifdef HAVE_DETAILED_TIMINGS
      call timer%stop("mpi_communication")
#endif
      max_blocks_row = ((na -1)/nblk)/np_rows + 1  ! Rows of A
      max_blocks_col = ((nqc-1)/nblk)/np_cols + 1  ! Columns of q!

      max_local_rows = max_blocks_row*nblk
      max_local_cols = max_blocks_col*nblk

      if (useGPU) then

        ! here the GPU and CPU version diverged: the CPU version now always uses the useQR path which
        ! is not implemented in the GPU version
        allocate(tmp1(max_local_cols*nbw), stat=istat, errmsg=errorMessage)
        if (istat .ne. 0) then
          print *,"trans_ev_band_to_full_real: error when allocating tmp1 "//errorMessage
          stop
        endif

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

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

        allocate(hvm(max_local_rows,nbw), stat=istat, errmsg=errorMessage)
        if (istat .ne. 0) then
          print *,"trans_ev_band_to_full_real: error when allocating hvm "//errorMessage
          stop
        endif
        successCUDA = cuda_malloc(hvm_dev, (max_local_rows)*nbw*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
          stop
        endif
        successCUDA = cuda_malloc(tmp_dev, (max_local_cols)*nbw*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
          stop
        endif

!#ifdef WITH_MPI
! it should be possible to keep tmat dev on the device and not copy it around
! already existent on GPU
        successCUDA = cuda_malloc(tmat_dev, nbw*nbw*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
          stop
        endif
!#endif

! q_dev already living on device
!        successCUDA = cuda_malloc(q_dev, ldq*matrixCols*M_size_of_PRECISSION_real_datatape)
!        if (.not.(successCUDA)) then
!          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
!          stop
!        endif

  !      q_temp(:,:) = 0.0
  !      q_temp(1:ldq,1:na_cols) = q(1:ldq,1:na_cols)

!        ! copy q_dev to device, maybe this can be avoided if q_dev can be kept on device in trans_ev_tridi_to_band
!        successCUDA = cuda_memcpy(q_dev, loc(q), (ldq)*(matrixCols)*M_size_of_PRECISSION_real_datatape, cudaMemcpyHostToDevice)
!        if (.not.(successCUDA)) then
!          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
!          stop
!        endif

        ! if MPI is NOT used the following steps could be done on the GPU and memory transfers could be avoided
        successCUDA = cuda_memset(hvm_dev, 0, (max_local_rows)*(nbw)*M_size_of_PRECISSION_real_datatape)
        if (.not.(successCUDA)) then
          print *,"trans_ev_band_to_full_real: error in cudaMalloc"
          stop
        endif
        hvm = M_CONST_0_0   ! Must be set to 0 !!!
        hvb = M_CONST_0_0   ! Safety only
        
        l_cols = local_index(nqc, my_pcol, np_cols, nblk, -1) ! Local columns of q

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

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

          ! Broadcast all Householder vectors for current step compressed in hvb

          nb = 0
          ns = 0

          do lc = 1, n_cols
            ncol = istep*nbw + lc ! absolute column number of householder vector
            nrow = ncol - nbw ! absolute number of pivot row

            l_rows = local_index(nrow-1, my_prow, np_rows, nblk, -1) ! row length for bcast
            l_colh = local_index(ncol  , my_pcol, np_cols, nblk, -1) ! HV local column number

            if (my_pcol==pcol(ncol, nblk, np_cols)) hvb(nb+1:nb+l_rows) = a(1:l_rows,l_colh)

            nb = nb+l_rows

            if (lc==n_cols .or. mod(ncol,nblk)==0) then
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
              call timer%start("mpi_communication")
#endif
              call MPI_Bcast(hvb(ns+1), nb-ns, M_MPI_REAL_PRECISSION, pcol(ncol, nblk, np_cols), mpi_comm_cols, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
              call timer%stop("mpi_communication")
#endif

#endif /* WITH_MPI */
              ns = nb
            endif
          enddo

          ! Expand compressed Householder vectors into matrix hvm

          nb = 0
          do lc = 1, n_cols
            nrow = (istep-1)*nbw+lc ! absolute number of pivot row
            l_rows = local_index(nrow-1, my_prow, np_rows, nblk, -1) ! row length for bcast

            hvm(1:l_rows,lc) = hvb(nb+1:nb+l_rows)
            if (my_prow==prow(nrow, nblk, np_rows)) hvm(l_rows+1,lc) = M_CONST_1_0

            nb = nb+l_rows
          enddo
          successCUDA = cuda_memcpy(hvm_dev, loc(hvm), ((max_local_rows)*nbw*M_size_of_PRECISSION_real_datatape),cudaMemcpyHostToDevice)

          if (.not.(successCUDA)) then
            print *,"trans_ev_band_to_full_real: error in cudaMemcpy"
            stop

          endif

          l_rows = local_index(MIN(na,(istep+1)*nbw), my_prow, np_rows, nblk, -1)

          ! Q = Q - V * T**T * V**T * Q

          if (l_rows>0) then
            call M_cublas_PRECISSION_gemm('T', 'N', n_cols, l_cols, l_rows, M_CONST_1_0, hvm_dev, max_local_rows, &
                              q_dev, ldq , M_CONST_0_0, tmp_dev, n_cols)

#ifdef WITH_MPI
            ! copy data from device to host for a later MPI_ALLREDUCE

            ! copy to host maybe this can be avoided this is needed if MPI is used (allreduce)
            successCUDA = cuda_memcpy(loc(tmp1), tmp_dev, l_cols*n_cols*M_size_of_PRECISSION_real_datatape, cudaMemcpyDeviceToHost)
            if (.not.(successCUDA)) then
              print *,"trans_ev_band_to_full_real: error in cudaMemcpy"
              stop
            endif
#endif /* WITH_MPI */

          else ! l_rows>0

#ifdef WITH_MPI
            tmp1(1:l_cols*n_cols) = 0
#else

            ! if MPI is not used (we do not need to transfer because of MPI_ALLREDUCE) we can set to zero on the device
#ifdef WITH_GPU_VERSION

            successCUDA = cuda_memset(tmp_dev, 0, l_cols*n_cols*M_size_of_PRECISSION_real_datatape)
            if (.not.(successCUDA)) then
              print *,"trans_ev_band_to_full_real: error in cudaMemset"
              stop
             endif
#endif

#endif /* WITH_MPI */

          endif ! l_rows>0

          !#ifdef WITH_GPU_VERSION
          !       istat = cuda_memcpy(loc(tmp1), tmp_dev, max_local_cols*nbw*size_of_real_datatype,cudaMemcpyDeviceToHost)
          !       if (istat .ne. 0) then
          !         print *,"error in cudaMemcpy"
          !         stop
          !       endif
          !#endif
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
          call timer%start("mpi_communication")
#endif
          call mpi_allreduce(tmp1, tmp2, n_cols*l_cols, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
          call timer%stop("mpi_communication")
#endif
#else /* WITH_MPI */
!          tmp2(1:n_cols*l_cols) = tmp1(1:n_cols*l_cols)
#endif /* WITH_MPI */
          !#ifdef WITH_GPU_VERSION
          !       istat = cuda_memcpy(tmp_dev, loc(tmp2), max_local_cols*nbw*size_of_real_datatype,cudaMemcpyHostToDevice)
          !       if (istat .ne. 0) then
          !         print *,"error in cudaMemcpy"
          !         stop
          !       endif
          !#endif

          if (l_rows>0) then
#ifdef WITH_MPI
            ! after the mpi_allreduce we have to copy back to the device
            ! copy back to device
            successCUDA = cuda_memcpy(tmp_dev, loc(tmp2), n_cols*l_cols*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
            if (.not.(successCUDA)) then
              print *,"trans_ev_band_to_full_real: error in cudaMemcpy"
              stop
            endif
#endif /* WITH_MPI */

!#ifdef WITH_MPI
! it should be possible to keep tmat on the device and not copy it aroud
!            ! copy to device, maybe this can be avoided tmat is input from bandred_real

            successCUDA = cuda_memcpy(tmat_dev, loc(tmat(1,1,istep)), nbw*nbw*M_size_of_PRECISSION_real_datatape,cudaMemcpyHostToDevice)
            if (.not.(successCUDA)) then
              print *,"trans_ev_band_to_full_real: error in cudaMemcpy"
              stop
            endif
!#endif /* WITH_MPI */

            call M_cublas_PRECISSION_trmm('L', 'U', 'T', 'N', n_cols, l_cols, M_CONST_1_0, tmat_dev, nbw, tmp_dev, n_cols)
            call M_cublas_PRECISSION_gemm('N', 'N', l_rows, l_cols, n_cols, -M_CONST_1_0, hvm_dev, max_local_rows, &
                              tmp_dev, n_cols, M_CONST_1_0, q_dev, ldq)

            ! copy to host maybe this can be avoided
            ! this is not necessary hvm is not used anymore
            successCUDA = cuda_memcpy(loc(hvm), hvm_dev, ((max_local_rows)*nbw*M_size_of_PRECISSION_real_datatape),cudaMemcpyDeviceToHost)
            if (.not.(successCUDA)) then
              print *,"trans_ev_band_to_full_real: error in cudaMemcpy"
              stop
            endif

          endif ! l_rows > 0
          !#ifdef WITH_GPU_VERSION
          !       istat = cuda_memcpy(loc(hvm), hvm_dev, ((max_local_rows)*nbw*size_of_real_datatype),cudaMemcpyDeviceToHost)
          !       if (istat .ne. 0) then
          !         print *,"error in cudaMemcpy"
          !         stop
          !       endif
          !
          !#endif
        enddo ! istep

      else ! do not useGPU

        ! t_blocking was formerly 2; 3 is a better choice
        t_blocking = 3 ! number of matrices T (tmat) which are aggregated into a new (larger) T matrix (tmat_complete) and applied at once

        ! we only use the t_blocking if we could call it fully, this is might be better but needs to benchmarked.
!       if ( na >= ((t_blocking+1)*nbw) ) then
        cwy_blocking = t_blocking * nbw

        allocate(tmp1(max_local_cols*cwy_blocking))
        allocate(tmp2(max_local_cols*cwy_blocking))
        allocate(hvb(max_local_rows*cwy_blocking))
        allocate(hvm(max_local_rows,cwy_blocking))
        allocate(tmat_complete(cwy_blocking,cwy_blocking))
        allocate(t_tmp(cwy_blocking,nbw))
        allocate(t_tmp2(cwy_blocking,nbw))
!        else
!          allocate(tmp1(max_local_cols*nbw))
!          allocate(tmp2(max_local_cols*nbw))
!          allocate(hvb(max_local_rows*nbw))
!          allocate(hvm(max_local_rows,nbw))
!        endif
        hvm = M_CONST_0_0   ! Must be set to 0 !!!
        hvb = M_CONST_0_0   ! Safety only
        l_cols = local_index(nqc, my_pcol, np_cols, nblk, -1) ! Local columns of q

!       if ( na >= ((t_blocking+1)*nbw) ) then

        do istep=1,((na-1)/nbw-1)/t_blocking + 1
          ! This the call when using  na >= ((t_blocking+1)*nbw)
          !      n_cols = MIN(na,istep*cwy_blocking+nbw) - (istep-1)*cwy_blocking - nbw ! Number of columns in current step
          ! As an alternative we add some special case handling if na < cwy_blocking
          IF (na < cwy_blocking) THEN
            n_cols = MAX(0, na-nbw)
            IF ( n_cols .eq. 0 ) THEN
              EXIT
            END IF
          ELSE
            n_cols = MIN(na,istep*cwy_blocking+nbw) - (istep-1)*cwy_blocking - nbw ! Number of columns in current step
          END IF
          ! Broadcast all Householder vectors for current step compressed in hvb

          nb = 0
          ns = 0

          do lc = 1, n_cols
            ncol = (istep-1)*cwy_blocking + nbw + lc ! absolute column number of householder vector
            nrow = ncol - nbw ! absolute number of pivot row

            l_rows = local_index(nrow-1, my_prow, np_rows, nblk, -1) ! row length for bcast
            l_colh = local_index(ncol  , my_pcol, np_cols, nblk, -1) ! HV local column number

            if (my_pcol==pcol(ncol, nblk, np_cols)) hvb(nb+1:nb+l_rows) = a(1:l_rows,l_colh)

            nb = nb+l_rows

            if (lc==n_cols .or. mod(ncol,nblk)==0) then
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
              call timer%start("mpi_communication")
#endif
              call MPI_Bcast(hvb(ns+1), nb-ns, M_MPI_REAL_PRECISSION, pcol(ncol, nblk, np_cols), mpi_comm_cols, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
              call timer%stop("mpi_communication")
#endif

#endif /* WITH_MPI */
              ns = nb
            endif
          enddo

          ! Expand compressed Householder vectors into matrix hvm

          nb = 0
          do lc = 1, n_cols
            nrow = (istep-1)*cwy_blocking + lc ! absolute number of pivot row
            l_rows = local_index(nrow-1, my_prow, np_rows, nblk, -1) ! row length for bcast

            hvm(1:l_rows,lc) = hvb(nb+1:nb+l_rows)
            if (my_prow==prow(nrow, nblk, np_rows)) hvm(l_rows+1,lc) = M_CONST_1_0

            nb = nb+l_rows
          enddo

          l_rows = local_index(MIN(na,(istep+1)*cwy_blocking), my_prow, np_rows, nblk, -1)

          ! compute tmat2 out of tmat(:,:,)
          tmat_complete = 0
          do i = 1, t_blocking
            t_cols = MIN(nbw, n_cols - (i-1)*nbw)
            if (t_cols <= 0) exit
            t_rows = (i - 1) * nbw
            tmat_complete(t_rows+1:t_rows+t_cols,t_rows+1:t_rows+t_cols) = tmat(1:t_cols,1:t_cols,(istep-1)*t_blocking + i)
            if (i > 1) then
              call M_PRECISSION_GEMM('T', 'N', t_rows, t_cols, l_rows, M_CONST_1_0, hvm(1,1), max_local_rows, hvm(1,(i-1)*nbw+1), &
                        max_local_rows, M_CONST_0_0, t_tmp, cwy_blocking)
#ifdef WITH_MPI
#ifdef HAVE_DETAILED_TIMINGS
              call timer%start("mpi_communication")
#endif

              call mpi_allreduce(t_tmp, t_tmp2, cwy_blocking*nbw, M_MPI_REAL_PRECISSION, MPI_SUM, mpi_comm_rows, mpierr)
#ifdef HAVE_DETAILED_TIMINGS
              call timer%stop("mpi_communication")
#endif

              call M_PRECISSION_TRMM('L', 'U', 'N', 'N', t_rows, t_cols, M_CONST_1_0, tmat_complete, cwy_blocking, t_tmp2, cwy_blocking)
              call M_PRECISSION_TRMM('R', 'U', 'N', 'N', t_rows, t_cols, -M_CONST_1_0, tmat_complete(t_rows+1,t_rows+1), cwy_blocking, &
                         t_tmp2, cwy_blocking)
              tmat_complete(1:t_rows,t_rows+1:t_rows+t_cols) = t_tmp2(1:t_rows,1:t_cols)

#else
!              t_tmp2(1:cwy_blocking,1:nbw) = t_tmp(1:cwy_blocking,1:nbw)

              call M_PRECISSION_TRMM('L', 'U', 'N', 'N', t_rows, t_cols, M_CONST_1_0, tmat_complete, cwy_blocking, t_tmp, cwy_blocking)
              call M_PRECISSION_TRMM('R',