! 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".

module ELPA1

! Version 1.1.2, 2011-02-21

  implicit none

  PRIVATE ! By default, all routines contained are private

  ! The following routines are public:

!  public :: get_elpa_row_col_comms     ! Sets MPI row/col communicators

!  public :: solve_evp_real             ! Driver routine for real eigenvalue problem
!  public :: solve_evp_complex          ! Driver routine for complex eigenvalue problem

  public :: tridiag_real               ! Transform real symmetric matrix to tridiagonal form
  public :: trans_ev_real              ! Transform eigenvectors of a tridiagonal matrix back
  public :: mult_at_b_real             ! Multiply real matrices A**T * B

  public :: tridiag_complex            ! Transform complex hermitian matrix to tridiagonal form
  public :: trans_ev_complex           ! Transform eigenvectors of a tridiagonal matrix back
  public :: mult_ah_b_complex          ! Multiply complex matrices A**H * B

  public :: solve_tridi                ! Solve tridiagonal eigensystem with divide and conquer method

  public :: cholesky_real              ! Cholesky factorization of a real matrix
  public :: invert_trm_real            ! Invert real triangular matrix

  public :: cholesky_complex           ! Cholesky factorization of a complex matrix
  public :: invert_trm_complex         ! Invert complex triangular matrix

  public :: local_index                ! Get local index of a block cyclic distributed matrix
  public :: least_common_multiple      ! Get least common multiple

  public :: hh_transform_real
  public :: hh_transform_complex

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

  ! Timing results, set by every call to solve_evp_xxx

  real*8, public :: time_evp_fwd    ! forward transformations (to tridiagonal form)
  real*8, public :: time_evp_solve  ! time for solving the tridiagonal system
  real*8, public :: time_evp_back   ! time for back transformations of eigenvectors

  ! Set elpa_print_times to .true. for explicit timing outputs

  logical, public :: elpa_print_times = .false.

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

  include 'mpif.h'

contains

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

subroutine tridiag_real(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols, d, e, tau)

!-------------------------------------------------------------------------------
!  tridiag_real: Reduces a distributed symmetric matrix to tridiagonal form
!                (like Scalapack Routine PDSYTRD)
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be reduced.
!              Distribution is like in Scalapack.
!              Opposed to PDSYTRD, a(:,:) must be set completely (upper and lower half)
!              a(:,:) is overwritten on exit with the Householder vectors
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!  d(na)       Diagonal elements (returned), identical on all processors
!
!  e(na)       Off-Diagonal elements (returned), identical on all processors
!
!  tau(na)     Factors for the Householder vectors (returned), needed for back transformation
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   real*8 a(lda,*), d(na), e(na), tau(na)

   integer, parameter :: max_stored_rows = 32

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer totalblocks, max_blocks_row, max_blocks_col, max_local_rows, max_local_cols
   integer l_cols, l_rows, nstor
   integer istep, i, j, lcs, lce, lrs, lre
   integer tile_size, l_rows_tile, l_cols_tile

   real*8 vav, vnorm2, x, aux(2*max_stored_rows), aux1(2), aux2(2), vrl, xf

   real*8, allocatable:: tmp(:), vr(:), vc(:), ur(:), uc(:), vur(:,:), uvc(:,:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   ! 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


   totalblocks = (na-1)/nblk + 1
   max_blocks_row = (totalblocks-1)/np_rows + 1
   max_blocks_col = (totalblocks-1)/np_cols + 1

   max_local_rows = max_blocks_row*nblk
   max_local_cols = max_blocks_col*nblk

   allocate(tmp(MAX(max_local_rows,max_local_cols)))
   allocate(vr(max_local_rows+1))
   allocate(ur(max_local_rows))
   allocate(vc(max_local_cols))
   allocate(uc(max_local_cols))

   tmp = 0
   vr = 0
   ur = 0
   vc = 0
   uc = 0

   allocate(vur(max_local_rows,2*max_stored_rows))
   allocate(uvc(max_local_cols,2*max_stored_rows))

   d(:) = 0
   e(:) = 0
   tau(:) = 0

   nstor = 0

   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a
   if(my_prow==prow(na) .and. my_pcol==pcol(na)) d(na) = a(l_rows,l_cols)

   do istep=na,3,-1

      ! Calculate number of local rows and columns of the still remaining matrix
      ! on the local processor

      l_rows = local_index(istep-1, my_prow, np_rows, nblk, -1)
      l_cols = local_index(istep-1, my_pcol, np_cols, nblk, -1)

      ! Calculate vector for Householder transformation on all procs
      ! owning column istep

      if(my_pcol==pcol(istep)) then

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

         vr(1:l_rows) = a(1:l_rows,l_cols+1)
         if(nstor>0 .and. l_rows>0) then
            call DGEMV('N',l_rows,2*nstor,1.d0,vur,ubound(vur,1), &
                       uvc(l_cols+1,1),ubound(uvc,1),1.d0,vr,1)
         endif

         if(my_prow==prow(istep-1)) then
            aux1(1) = dot_product(vr(1:l_rows-1),vr(1:l_rows-1))
            aux1(2) = vr(l_rows)
         else
            aux1(1) = dot_product(vr(1:l_rows),vr(1:l_rows))
            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(istep))

         ! Scale vr and store Householder vector for back transformation

         vr(1:l_rows) = vr(1:l_rows) * xf
         if(my_prow==prow(istep-1)) then
            vr(l_rows) = 1.
            e(istep-1) = vrl
         endif
         a(1:l_rows,l_cols+1) = vr(1:l_rows) ! store Householder vector for back transformation

      endif

      ! Broadcast the Householder vector (and tau) along columns

      if(my_pcol==pcol(istep)) vr(l_rows+1) = tau(istep)
      call MPI_Bcast(vr,l_rows+1,MPI_REAL8,pcol(istep),mpi_comm_cols,mpierr)
      tau(istep) =  vr(l_rows+1)

      ! Transpose Householder vector vr -> vc

      call elpa_transpose_vectors  (vr, ubound(vr,1), mpi_comm_rows, &
                                    vc, ubound(vc,1), mpi_comm_cols, &
                                    1, istep-1, 1, nblk)


      ! Calculate u = (A + VU**T + UV**T)*v

      ! For cache efficiency, we use only the upper half of the matrix tiles for this,
      ! thus the result is partly in uc(:) and partly in ur(:)

      uc(1:l_cols) = 0
      ur(1:l_rows) = 0
      if(l_rows>0 .and. l_cols>0) then

         do i=0,(istep-2)/tile_size
            lcs = i*l_cols_tile+1
            lce = min(l_cols,(i+1)*l_cols_tile)
            if(lce<lcs) cycle
            do j=0,i
               lrs = j*l_rows_tile+1
               lre = min(l_rows,(j+1)*l_rows_tile)
               if(lre<lrs) cycle
               call DGEMV('T',lre-lrs+1,lce-lcs+1,1.d0,a(lrs,lcs),lda,vr(lrs),1,1.d0,uc(lcs),1)
               if(i/=j) call DGEMV('N',lre-lrs+1,lce-lcs+1,1.d0,a(lrs,lcs),lda,vc(lcs),1,1.d0,ur(lrs),1)
            enddo
         enddo

         if(nstor>0) then
            call DGEMV('T',l_rows,2*nstor,1.d0,vur,ubound(vur,1),vr,1,0.d0,aux,1)
            call DGEMV('N',l_cols,2*nstor,1.d0,uvc,ubound(uvc,1),aux,1,1.d0,uc,1)
         endif

      endif

      ! 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(tile_size < istep-1) then
         call elpa_reduce_add_vectors  (ur, ubound(ur,1), mpi_comm_rows, &
                                        uc, ubound(uc,1), mpi_comm_cols, &
                                        istep-1, 1, nblk)
      endif

      ! Sum up all the uc(:) parts, transpose uc -> ur

      if(l_cols>0) then
         tmp(1:l_cols) = uc(1:l_cols)
         call mpi_allreduce(tmp,uc,l_cols,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
      endif

      call elpa_transpose_vectors  (uc, ubound(uc,1), mpi_comm_cols, &
                                    ur, ubound(ur,1), mpi_comm_rows, &
                                    1, istep-1, 1, nblk)

      ! calculate u**T * v (same as v**T * (A + VU**T + UV**T) * v )

      x = 0
      if(l_cols>0) x = dot_product(vc(1:l_cols),uc(1:l_cols))
      call mpi_allreduce(x,vav,1,MPI_REAL8,MPI_SUM,mpi_comm_cols,mpierr)

      ! store u and v in the matrices U and V
      ! these matrices are stored combined in one here

      do j=1,l_rows
         vur(j,2*nstor+1) = tau(istep)*vr(j)
         vur(j,2*nstor+2) = 0.5*tau(istep)*vav*vr(j) - ur(j)
      enddo
      do j=1,l_cols
         uvc(j,2*nstor+1) = 0.5*tau(istep)*vav*vc(j) - uc(j)
         uvc(j,2*nstor+2) = tau(istep)*vc(j)
      enddo

      nstor = nstor+1

      ! If the limit of max_stored_rows is reached, calculate A + VU**T + UV**T

      if(nstor==max_stored_rows .or. istep==3) then

         do i=0,(istep-2)/tile_size
            lcs = i*l_cols_tile+1
            lce = min(l_cols,(i+1)*l_cols_tile)
            lrs = 1
            lre = min(l_rows,(i+1)*l_rows_tile)
            if(lce<lcs .or. lre<lrs) cycle
            call dgemm('N','T',lre-lrs+1,lce-lcs+1,2*nstor,1.d0, &
                       vur(lrs,1),ubound(vur,1),uvc(lcs,1),ubound(uvc,1), &
                       1.d0,a(lrs,lcs),lda)
         enddo

         nstor = 0

      endif

      if(my_prow==prow(istep-1) .and. my_pcol==pcol(istep-1)) then
         if(nstor>0) a(l_rows,l_cols) = a(l_rows,l_cols) &
                        + dot_product(vur(l_rows,1:2*nstor),uvc(l_cols,1:2*nstor))
         d(istep-1) = a(l_rows,l_cols)
      endif

   enddo

   ! Store e(1) and d(1)

   if(my_prow==prow(1) .and. my_pcol==pcol(2)) e(1) = a(1,l_cols) ! use last l_cols value of loop above
   if(my_prow==prow(1) .and. my_pcol==pcol(1)) d(1) = a(1,1)

   deallocate(tmp, vr, ur, vc, uc, vur, uvc)

   ! distribute the arrays d and e to all processors

   allocate(tmp(na))
   tmp = d
   call mpi_allreduce(tmp,d,na,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
   tmp = d
   call mpi_allreduce(tmp,d,na,MPI_REAL8,MPI_SUM,mpi_comm_cols,mpierr)
   tmp = e
   call mpi_allreduce(tmp,e,na,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
   tmp = e
   call mpi_allreduce(tmp,e,na,MPI_REAL8,MPI_SUM,mpi_comm_cols,mpierr)
   deallocate(tmp)

end subroutine tridiag_real

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

subroutine trans_ev_real(na, nqc, a, lda, tau, q, ldq, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  trans_ev_real: Transforms the eigenvectors of a tridiagonal matrix back
!                 to the eigenvectors of the original matrix
!                 (like Scalapack Routine PDORMTR)
!
!  Parameters
!
!  na          Order of matrix a, number of rows of matrix q
!
!  nqc         Number of columns of matrix q
!
!  a(lda,*)    Matrix containing the Householder vectors (i.e. matrix a after tridiag_real)
!              Distribution is like in Scalapack.
!
!  lda         Leading dimension of a
!
!  tau(na)     Factors of the Householder vectors
!
!  q           On input: Eigenvectors of tridiagonal matrix
!              On output: Transformed eigenvectors
!              Distribution is like in Scalapack.
!
!  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
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, nqc, lda, ldq, nblk, mpi_comm_rows, mpi_comm_cols
   real*8 a(lda,*), q(ldq,*), tau(na)

   integer :: max_stored_rows

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer totalblocks, max_blocks_row, max_blocks_col, max_local_rows, max_local_cols
   integer l_cols, l_rows, l_colh, nstor
   integer istep, i, n, nc, ic, ics, ice, nb, cur_pcol

   real*8, allocatable:: tmp1(:), tmp2(:), hvb(:), hvm(:,:)
   real*8, allocatable:: tmat(:,:), h1(:), h2(:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)


   totalblocks = (na-1)/nblk + 1
   max_blocks_row = (totalblocks-1)/np_rows + 1
   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


   max_stored_rows = (63/nblk+1)*nblk

   allocate(tmat(max_stored_rows,max_stored_rows))
   allocate(h1(max_stored_rows*max_stored_rows))
   allocate(h2(max_stored_rows*max_stored_rows))
   allocate(tmp1(max_local_cols*max_stored_rows))
   allocate(tmp2(max_local_cols*max_stored_rows))
   allocate(hvb(max_local_rows*nblk))
   allocate(hvm(max_local_rows,max_stored_rows))

   hvm = 0   ! Must be set to 0 !!!
   hvb = 0   ! Safety only

   l_cols = local_index(nqc, my_pcol, np_cols, nblk, -1) ! Local columns of q

   nstor = 0

   do istep=1,na,nblk

      ics = MAX(istep,3)
      ice = MIN(istep+nblk-1,na)
      if(ice<ics) cycle

      cur_pcol = pcol(istep)

      nb = 0
      do ic=ics,ice

         l_colh = local_index(ic  , my_pcol, np_cols, nblk, -1) ! Column of Householder vector
         l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector


         if(my_pcol==cur_pcol) then
            hvb(nb+1:nb+l_rows) = a(1:l_rows,l_colh)
            if(my_prow==prow(ic-1)) then
               hvb(nb+l_rows) = 1.
            endif
         endif

         nb = nb+l_rows
      enddo

      if(nb>0) &
         call MPI_Bcast(hvb,nb,MPI_REAL8,cur_pcol,mpi_comm_cols,mpierr)

      nb = 0
      do ic=ics,ice
         l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector
         hvm(1:l_rows,nstor+1) = hvb(nb+1:nb+l_rows)
         nstor = nstor+1
         nb = nb+l_rows
      enddo

      ! Please note: for smaller matix sizes (na/np_rows<=256), a value of 32 for nstor is enough!
      if(nstor+nblk>max_stored_rows .or. istep+nblk>na .or. (na/np_rows<=256 .and. nstor>=32)) then

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

         tmat = 0
         if(l_rows>0) &
            call dsyrk('U','T',nstor,l_rows,1.d0,hvm,ubound(hvm,1),0.d0,tmat,max_stored_rows)

         nc = 0
         do n=1,nstor-1
            h1(nc+1:nc+n) = tmat(1:n,n+1)
            nc = nc+n
         enddo

         if(nc>0) call mpi_allreduce(h1,h2,nc,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)

         ! Calculate triangular matrix T

         nc = 0
         tmat(1,1) = tau(ice-nstor+1)
         do n=1,nstor-1
            call dtrmv('L','T','N',n,tmat,max_stored_rows,h2(nc+1),1)
            tmat(n+1,1:n) = -h2(nc+1:nc+n)*tau(ice-nstor+n+1)
            tmat(n+1,n+1) = tau(ice-nstor+n+1)
            nc = nc+n
         enddo

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

         if(l_rows>0) then
            call dgemm('T','N',nstor,l_cols,l_rows,1.d0,hvm,ubound(hvm,1), &
                       q,ldq,0.d0,tmp1,nstor)
         else
            tmp1(1:l_cols*nstor) = 0
         endif
         call mpi_allreduce(tmp1,tmp2,nstor*l_cols,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
         if(l_rows>0) then
            call dtrmm('L','L','N','N',nstor,l_cols,1.0d0,tmat,max_stored_rows,tmp2,nstor)
            call dgemm('N','N',l_rows,l_cols,nstor,-1.d0,hvm,ubound(hvm,1), &
                       tmp2,nstor,1.d0,q,ldq)
         endif
         nstor = 0
      endif

   enddo

   deallocate(tmat, h1, h2, tmp1, tmp2, hvb, hvm)


end subroutine trans_ev_real

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

subroutine mult_at_b_real(uplo_a, uplo_c, na, ncb, a, lda, b, ldb, nblk, mpi_comm_rows, mpi_comm_cols, c, ldc)

!-------------------------------------------------------------------------------
!  mult_at_b_real:  Performs C := A**T * B
!
!      where:  A is a square matrix (na,na) which is optionally upper or lower triangular
!              B is a (na,ncb) matrix
!              C is a (na,ncb) matrix where optionally only the upper or lower
!              triangle may be computed
!
!  Parameters
!
!  uplo_a      'U' if A is upper triangular
!              'L' if A is lower triangular
!              anything else if A is a full matrix
!              Please note: This pertains to the original A (as set in the calling program)
!              whereas the transpose of A is used for calculations
!              If uplo_a is 'U' or 'L', the other triangle is not used at all,
!              i.e. it may contain arbitrary numbers
!
!  uplo_c      'U' if only the upper diagonal part of C is needed
!              'L' if only the upper diagonal part of C is needed
!              anything else if the full matrix C is needed
!              Please note: Even when uplo_c is 'U' or 'L', the other triangle may be
!              written to a certain extent, i.e. one shouldn't rely on the content there!
!
!  na          Number of rows/columns of A, number of rows of B and C
!
!  ncb         Number of columns  of B and C
!
!  a           Matrix A
!
!  lda         Leading dimension of a
!
!  b           Matrix B
!
!  ldb         Leading dimension of b
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!  c           Matrix C
!
!  ldc         Leading dimension of c
!
!-------------------------------------------------------------------------------

   implicit none

   character*1 uplo_a, uplo_c

   integer na, ncb, lda, ldb, nblk, mpi_comm_rows, mpi_comm_cols, ldc
   real*8 a(lda,*), b(ldb,*), c(ldc,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_rows_np
   integer np, n, nb, nblk_mult, lrs, lre, lcs, lce
   integer gcol_min, gcol, goff
   integer nstor, nr_done, noff, np_bc, n_aux_bc, nvals
   integer, allocatable :: lrs_save(:), lre_save(:)

   logical a_lower, a_upper, c_lower, c_upper

   real*8, allocatable:: aux_mat(:,:), aux_bc(:), tmp1(:,:), tmp2(:,:)


   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)

   l_rows = local_index(na,  my_prow, np_rows, nblk, -1) ! Local rows of a and b
   l_cols = local_index(ncb, my_pcol, np_cols, nblk, -1) ! Local cols of b

   ! Block factor for matrix multiplications, must be a multiple of nblk

   if(na/np_rows<=256) then
      nblk_mult = (31/nblk+1)*nblk
   else
      nblk_mult = (63/nblk+1)*nblk
   endif

   allocate(aux_mat(l_rows,nblk_mult))
   allocate(aux_bc(l_rows*nblk))
   allocate(lrs_save(nblk))
   allocate(lre_save(nblk))

   a_lower = .false.
   a_upper = .false.
   c_lower = .false.
   c_upper = .false.

   if(uplo_a=='u' .or. uplo_a=='U') a_upper = .true.
   if(uplo_a=='l' .or. uplo_a=='L') a_lower = .true.
   if(uplo_c=='u' .or. uplo_c=='U') c_upper = .true.
   if(uplo_c=='l' .or. uplo_c=='L') c_lower = .true.

   ! Build up the result matrix by processor rows

   do np = 0, np_rows-1

      ! In this turn, procs of row np assemble the result

      l_rows_np = local_index(na, np, np_rows, nblk, -1) ! local rows on receiving processors

      nr_done = 0 ! Number of rows done
      aux_mat = 0
      nstor = 0   ! Number of columns stored in aux_mat

      ! Loop over the blocks on row np

      do nb=0,(l_rows_np-1)/nblk

         goff  = nb*np_rows + np ! Global offset in blocks corresponding to nb

         ! Get the processor column which owns this block (A is transposed, so we need the column)
         ! and the offset in blocks within this column.
         ! The corresponding block column in A is then broadcast to all for multiplication with B

         np_bc = MOD(goff,np_cols)
         noff = goff/np_cols
         n_aux_bc = 0

         ! Gather up the complete block column of A on the owner

         do n = 1, min(l_rows_np-nb*nblk,nblk) ! Loop over columns to be broadcast

            gcol = goff*nblk + n ! global column corresponding to n
            if(nstor==0 .and. n==1) gcol_min = gcol

            lrs = 1       ! 1st local row number for broadcast
            lre = l_rows  ! last local row number for broadcast
            if(a_lower) lrs = local_index(gcol, my_prow, np_rows, nblk, +1)
            if(a_upper) lre = local_index(gcol, my_prow, np_rows, nblk, -1)

            if(lrs<=lre) then
               nvals = lre-lrs+1
               if(my_pcol == np_bc) aux_bc(n_aux_bc+1:n_aux_bc+nvals) = a(lrs:lre,noff*nblk+n)
               n_aux_bc = n_aux_bc + nvals
            endif

            lrs_save(n) = lrs
            lre_save(n) = lre

         enddo

         ! Broadcast block column

         call MPI_Bcast(aux_bc,n_aux_bc,MPI_REAL8,np_bc,mpi_comm_cols,mpierr)

         ! Insert what we got in aux_mat

         n_aux_bc = 0
         do n = 1, min(l_rows_np-nb*nblk,nblk)
            nstor = nstor+1
            lrs = lrs_save(n)
            lre = lre_save(n)
            if(lrs<=lre) then
               nvals = lre-lrs+1
               aux_mat(lrs:lre,nstor) = aux_bc(n_aux_bc+1:n_aux_bc+nvals)
               n_aux_bc = n_aux_bc + nvals
            endif
         enddo

         ! If we got nblk_mult columns in aux_mat or this is the last block
         ! do the matrix multiplication

         if(nstor==nblk_mult .or. nb*nblk+nblk >= l_rows_np) then

            lrs = 1       ! 1st local row number for multiply
            lre = l_rows  ! last local row number for multiply
            if(a_lower) lrs = local_index(gcol_min, my_prow, np_rows, nblk, +1)
            if(a_upper) lre = local_index(gcol, my_prow, np_rows, nblk, -1)

            lcs = 1       ! 1st local col number for multiply
            lce = l_cols  ! last local col number for multiply
            if(c_upper) lcs = local_index(gcol_min, my_pcol, np_cols, nblk, +1)
            if(c_lower) lce = MIN(local_index(gcol, my_pcol, np_cols, nblk, -1),l_cols)

            if(lcs<=lce) then
               allocate(tmp1(nstor,lcs:lce),tmp2(nstor,lcs:lce))
               if(lrs<=lre) then
                  call dgemm('T','N',nstor,lce-lcs+1,lre-lrs+1,1.d0,aux_mat(lrs,1),ubound(aux_mat,1), &
                             b(lrs,lcs),ldb,0.d0,tmp1,nstor)
               else
                  tmp1 = 0
               endif

               ! Sum up the results and send to processor row np
               call mpi_reduce(tmp1,tmp2,nstor*(lce-lcs+1),MPI_REAL8,MPI_SUM,np,mpi_comm_rows,mpierr)

               ! Put the result into C
               if(my_prow==np) c(nr_done+1:nr_done+nstor,lcs:lce) = tmp2(1:nstor,lcs:lce)

               deallocate(tmp1,tmp2)
            endif

            nr_done = nr_done+nstor
            nstor=0
            aux_mat(:,:)=0
         endif
      enddo
   enddo

   deallocate(aux_mat, aux_bc, lrs_save, lre_save)

end subroutine mult_at_b_real

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

subroutine tridiag_complex(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols, d, e, tau)

!-------------------------------------------------------------------------------
!  tridiag_complex: Reduces a distributed hermitian matrix to tridiagonal form
!                   (like Scalapack Routine PZHETRD)
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be reduced.
!              Distribution is like in Scalapack.
!              Opposed to PZHETRD, a(:,:) must be set completely (upper and lower half)
!              a(:,:) is overwritten on exit with the Householder vectors
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!  d(na)       Diagonal elements (returned), identical on all processors
!
!  e(na)       Off-Diagonal elements (returned), identical on all processors
!
!  tau(na)     Factors for the Householder vectors (returned), needed for back transformation
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   complex*16 a(lda,*), tau(na)
   real*8 d(na), e(na)

   integer, parameter :: max_stored_rows = 32

   complex*16, parameter :: CZERO = (0.d0,0.d0), CONE = (1.d0,0.d0)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer totalblocks, max_blocks_row, max_blocks_col, max_local_rows, max_local_cols
   integer l_cols, l_rows, nstor
   integer istep, i, j, lcs, lce, lrs, lre
   integer tile_size, l_rows_tile, l_cols_tile

   real*8 vnorm2
   complex*16 vav, xc, aux(2*max_stored_rows),  aux1(2), aux2(2), vrl, xf

   complex*16, allocatable:: tmp(:), vr(:), vc(:), ur(:), uc(:), vur(:,:), uvc(:,:)
   real*8, allocatable:: tmpr(:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   ! 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


   totalblocks = (na-1)/nblk + 1
   max_blocks_row = (totalblocks-1)/np_rows + 1
   max_blocks_col = (totalblocks-1)/np_cols + 1

   max_local_rows = max_blocks_row*nblk
   max_local_cols = max_blocks_col*nblk

   allocate(tmp(MAX(max_local_rows,max_local_cols)))
   allocate(vr(max_local_rows+1))
   allocate(ur(max_local_rows))
   allocate(vc(max_local_cols))
   allocate(uc(max_local_cols))

   tmp = 0
   vr = 0
   ur = 0
   vc = 0
   uc = 0

   allocate(vur(max_local_rows,2*max_stored_rows))
   allocate(uvc(max_local_cols,2*max_stored_rows))

   d(:) = 0
   e(:) = 0
   tau(:) = 0

   nstor = 0

   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a
   if(my_prow==prow(na) .and. my_pcol==pcol(na)) d(na) = a(l_rows,l_cols)

   do istep=na,3,-1

      ! Calculate number of local rows and columns of the still remaining matrix
      ! on the local processor

      l_rows = local_index(istep-1, my_prow, np_rows, nblk, -1)
      l_cols = local_index(istep-1, my_pcol, np_cols, nblk, -1)

      ! Calculate vector for Householder transformation on all procs
      ! owning column istep

      if(my_pcol==pcol(istep)) then

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

         vr(1:l_rows) = a(1:l_rows,l_cols+1)
         if(nstor>0 .and. l_rows>0) then
            aux(1:2*nstor) = conjg(uvc(l_cols+1,1:2*nstor))
            call ZGEMV('N',l_rows,2*nstor,CONE,vur,ubound(vur,1), &
                       aux,1,CONE,vr,1)
         endif

         if(my_prow==prow(istep-1)) then
            aux1(1) = dot_product(vr(1:l_rows-1),vr(1:l_rows-1))
            aux1(2) = vr(l_rows)
         else
            aux1(1) = dot_product(vr(1:l_rows),vr(1:l_rows))
            aux1(2) = 0.
         endif

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

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

         ! Householder transformation

         call hh_transform_complex(vrl, vnorm2, xf, tau(istep))

         ! Scale vr and store Householder vector for back transformation

         vr(1:l_rows) = vr(1:l_rows) * xf
         if(my_prow==prow(istep-1)) then
            vr(l_rows) = 1.
            e(istep-1) = vrl
         endif
         a(1:l_rows,l_cols+1) = vr(1:l_rows) ! store Householder vector for back transformation

      endif

      ! Broadcast the Householder vector (and tau) along columns

      if(my_pcol==pcol(istep)) vr(l_rows+1) = tau(istep)
      call MPI_Bcast(vr,l_rows+1,MPI_DOUBLE_COMPLEX,pcol(istep),mpi_comm_cols,mpierr)
      tau(istep) =  vr(l_rows+1)

      ! Transpose Householder vector vr -> vc

      call elpa_transpose_vectors  (vr, 2*ubound(vr,1), mpi_comm_rows, &
                                    vc, 2*ubound(vc,1), mpi_comm_cols, &
                                    1, 2*(istep-1), 1, 2*nblk)

      ! Calculate u = (A + VU**T + UV**T)*v

      ! For cache efficiency, we use only the upper half of the matrix tiles for this,
      ! thus the result is partly in uc(:) and partly in ur(:)

      uc(1:l_cols) = 0
      ur(1:l_rows) = 0
      if(l_rows>0 .and. l_cols>0) then

         do i=0,(istep-2)/tile_size
            lcs = i*l_cols_tile+1
            lce = min(l_cols,(i+1)*l_cols_tile)
            if(lce<lcs) cycle
            do j=0,i
               lrs = j*l_rows_tile+1
               lre = min(l_rows,(j+1)*l_rows_tile)
               if(lre<lrs) cycle
               call ZGEMV('C',lre-lrs+1,lce-lcs+1,CONE,a(lrs,lcs),lda,vr(lrs),1,CONE,uc(lcs),1)
               if(i/=j) call ZGEMV('N',lre-lrs+1,lce-lcs+1,CONE,a(lrs,lcs),lda,vc(lcs),1,CONE,ur(lrs),1)
            enddo
         enddo

         if(nstor>0) then
            call ZGEMV('C',l_rows,2*nstor,CONE,vur,ubound(vur,1),vr,1,CZERO,aux,1)
            call ZGEMV('N',l_cols,2*nstor,CONE,uvc,ubound(uvc,1),aux,1,CONE,uc,1)
         endif

      endif

      ! 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(tile_size < istep-1) then
         call elpa_reduce_add_vectors  (ur, 2*ubound(ur,1), mpi_comm_rows, &
                                        uc, 2*ubound(uc,1), mpi_comm_cols, &
                                        2*(istep-1), 1, 2*nblk)
      endif

      ! Sum up all the uc(:) parts, transpose uc -> ur

      if(l_cols>0) then
         tmp(1:l_cols) = uc(1:l_cols)
         call mpi_allreduce(tmp,uc,l_cols,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_comm_rows,mpierr)
      endif

      call elpa_transpose_vectors  (uc, 2*ubound(uc,1), mpi_comm_cols, &
                                    ur, 2*ubound(ur,1), mpi_comm_rows, &
                                    1, 2*(istep-1), 1, 2*nblk)

      ! calculate u**T * v (same as v**T * (A + VU**T + UV**T) * v )

      xc = 0
      if(l_cols>0) xc = dot_product(vc(1:l_cols),uc(1:l_cols))
      call mpi_allreduce(xc,vav,1,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_comm_cols,mpierr)

      ! store u and v in the matrices U and V
      ! these matrices are stored combined in one here

      do j=1,l_rows
         vur(j,2*nstor+1) = conjg(tau(istep))*vr(j)
         vur(j,2*nstor+2) = 0.5*conjg(tau(istep))*vav*vr(j) - ur(j)
      enddo
      do j=1,l_cols
         uvc(j,2*nstor+1) = 0.5*conjg(tau(istep))*vav*vc(j) - uc(j)
         uvc(j,2*nstor+2) = conjg(tau(istep))*vc(j)
      enddo

      nstor = nstor+1

      ! If the limit of max_stored_rows is reached, calculate A + VU**T + UV**T

      if(nstor==max_stored_rows .or. istep==3) then

         do i=0,(istep-2)/tile_size
            lcs = i*l_cols_tile+1
            lce = min(l_cols,(i+1)*l_cols_tile)
            lrs = 1
            lre = min(l_rows,(i+1)*l_rows_tile)
            if(lce<lcs .or. lre<lrs) cycle
            call ZGEMM('N','C',lre-lrs+1,lce-lcs+1,2*nstor,CONE, &
                       vur(lrs,1),ubound(vur,1),uvc(lcs,1),ubound(uvc,1), &
                       CONE,a(lrs,lcs),lda)
         enddo

         nstor = 0

      endif

      if(my_prow==prow(istep-1) .and. my_pcol==pcol(istep-1)) then
         if(nstor>0) a(l_rows,l_cols) = a(l_rows,l_cols) &
                        + dot_product(vur(l_rows,1:2*nstor),uvc(l_cols,1:2*nstor))
         d(istep-1) = a(l_rows,l_cols)
      endif

   enddo

   ! Store e(1) and d(1)

   if(my_pcol==pcol(2)) then
      if(my_prow==prow(1)) then
         ! We use last l_cols value of loop above
         vrl = a(1,l_cols)
         call hh_transform_complex(vrl, 0.d0, xf, tau(2))
         e(1) = vrl
         a(1,l_cols) = 1. ! for consistency only
      endif
      call mpi_bcast(tau(2),1,MPI_DOUBLE_COMPLEX,prow(1),mpi_comm_rows,mpierr)
   endif
   call mpi_bcast(tau(2),1,MPI_DOUBLE_COMPLEX,pcol(2),mpi_comm_cols,mpierr)

   if(my_prow==prow(1) .and. my_pcol==pcol(1)) d(1) = a(1,1)

   deallocate(tmp, vr, ur, vc, uc, vur, uvc)

   ! distribute the arrays d and e to all processors

   allocate(tmpr(na))
   tmpr = d
   call mpi_allreduce(tmpr,d,na,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
   tmpr = d
   call mpi_allreduce(tmpr,d,na,MPI_REAL8,MPI_SUM,mpi_comm_cols,mpierr)
   tmpr = e
   call mpi_allreduce(tmpr,e,na,MPI_REAL8,MPI_SUM,mpi_comm_rows,mpierr)
   tmpr = e
   call mpi_allreduce(tmpr,e,na,MPI_REAL8,MPI_SUM,mpi_comm_cols,mpierr)
   deallocate(tmpr)

end subroutine tridiag_complex

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

subroutine trans_ev_complex(na, nqc, a, lda, tau, q, ldq, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  trans_ev_complex: Transforms the eigenvectors of a tridiagonal matrix back
!                    to the eigenvectors of the original matrix
!                    (like Scalapack Routine PZUNMTR)
!
!  Parameters
!
!  na          Order of matrix a, number of rows of matrix q
!
!  nqc         Number of columns of matrix q
!
!  a(lda,*)    Matrix containing the Householder vectors (i.e. matrix a after tridiag_complex)
!              Distribution is like in Scalapack.
!
!  lda         Leading dimension of a
!
!  tau(na)     Factors of the Householder vectors
!
!  q           On input: Eigenvectors of tridiagonal matrix
!              On output: Transformed eigenvectors
!              Distribution is like in Scalapack.
!
!  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
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, nqc, lda, ldq, nblk, mpi_comm_rows, mpi_comm_cols
   complex*16 a(lda,*), q(ldq,*), tau(na)

   integer :: max_stored_rows

   complex*16, parameter :: CZERO = (0.d0,0.d0), CONE = (1.d0,0.d0)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer totalblocks, max_blocks_row, max_blocks_col, max_local_rows, max_local_cols
   integer l_cols, l_rows, l_colh, nstor
   integer istep, i, n, nc, ic, ics, ice, nb, cur_pcol

   complex*16, allocatable:: tmp1(:), tmp2(:), hvb(:), hvm(:,:)
   complex*16, allocatable:: tmat(:,:), h1(:), h2(:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)


   totalblocks = (na-1)/nblk + 1
   max_blocks_row = (totalblocks-1)/np_rows + 1
   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


   max_stored_rows = (63/nblk+1)*nblk

   allocate(tmat(max_stored_rows,max_stored_rows))
   allocate(h1(max_stored_rows*max_stored_rows))
   allocate(h2(max_stored_rows*max_stored_rows))
   allocate(tmp1(max_local_cols*max_stored_rows))
   allocate(tmp2(max_local_cols*max_stored_rows))
   allocate(hvb(max_local_rows*nblk))
   allocate(hvm(max_local_rows,max_stored_rows))

   hvm = 0   ! Must be set to 0 !!!
   hvb = 0   ! Safety only

   l_cols = local_index(nqc, my_pcol, np_cols, nblk, -1) ! Local columns of q

   nstor = 0

   ! In the complex case tau(2) /= 0
   if(my_prow == prow(1)) then
      q(1,1:l_cols) = q(1,1:l_cols)*((1.d0,0.d0)-tau(2))
   endif

   do istep=1,na,nblk

      ics = MAX(istep,3)
      ice = MIN(istep+nblk-1,na)
      if(ice<ics) cycle

      cur_pcol = pcol(istep)

      nb = 0
      do ic=ics,ice

         l_colh = local_index(ic  , my_pcol, np_cols, nblk, -1) ! Column of Householder vector
         l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector


         if(my_pcol==cur_pcol) then
            hvb(nb+1:nb+l_rows) = a(1:l_rows,l_colh)
            if(my_prow==prow(ic-1)) then
               hvb(nb+l_rows) = 1.
            endif
         endif

         nb = nb+l_rows
      enddo

      if(nb>0) &
         call MPI_Bcast(hvb,nb,MPI_DOUBLE_COMPLEX,cur_pcol,mpi_comm_cols,mpierr)

      nb = 0
      do ic=ics,ice
         l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector
         hvm(1:l_rows,nstor+1) = hvb(nb+1:nb+l_rows)
         nstor = nstor+1
         nb = nb+l_rows
      enddo

      ! Please note: for smaller matix sizes (na/np_rows<=256), a value of 32 for nstor is enough!
      if(nstor+nblk>max_stored_rows .or. istep+nblk>na .or. (na/np_rows<=256 .and. nstor>=32)) then

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

         tmat = 0
         if(l_rows>0) &
            call zherk('U','C',nstor,l_rows,CONE,hvm,ubound(hvm,1),CZERO,tmat,max_stored_rows)

         nc = 0
         do n=1,nstor-1
            h1(nc+1:nc+n) = tmat(1:n,n+1)
            nc = nc+n
         enddo

         if(nc>0) call mpi_allreduce(h1,h2,nc,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_comm_rows,mpierr)

         ! Calculate triangular matrix T

         nc = 0
         tmat(1,1) = tau(ice-nstor+1)
         do n=1,nstor-1
            call ztrmv('L','C','N',n,tmat,max_stored_rows,h2(nc+1),1)
            tmat(n+1,1:n) = -conjg(h2(nc+1:nc+n))*tau(ice-nstor+n+1)
            tmat(n+1,n+1) = tau(ice-nstor+n+1)
            nc = nc+n
         enddo

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

         if(l_rows>0) then
            call zgemm('C','N',nstor,l_cols,l_rows,CONE,hvm,ubound(hvm,1), &
                       q,ldq,CZERO,tmp1,nstor)
         else
            tmp1(1:l_cols*nstor) = 0
         endif
         call mpi_allreduce(tmp1,tmp2,nstor*l_cols,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_comm_rows,mpierr)
         if(l_rows>0) then
            call ztrmm('L','L','N','N',nstor,l_cols,CONE,tmat,max_stored_rows,tmp2,nstor)
            call zgemm('N','N',l_rows,l_cols,nstor,-CONE,hvm,ubound(hvm,1), &
                       tmp2,nstor,CONE,q,ldq)
         endif
         nstor = 0
      endif

   enddo

   deallocate(tmat, h1, h2, tmp1, tmp2, hvb, hvm)


end subroutine trans_ev_complex

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

subroutine mult_ah_b_complex(uplo_a, uplo_c, na, ncb, a, lda, b, ldb, nblk, mpi_comm_rows, mpi_comm_cols, c, ldc)

!-------------------------------------------------------------------------------
!  mult_ah_b_complex:  Performs C := A**H * B
!
!      where:  A is a square matrix (na,na) which is optionally upper or lower triangular
!              B is a (na,ncb) matrix
!              C is a (na,ncb) matrix where optionally only the upper or lower
!              triangle may be computed
!
!  Parameters
!
!  uplo_a      'U' if A is upper triangular
!              'L' if A is lower triangular
!              anything else if A is a full matrix
!              Please note: This pertains to the original A (as set in the calling program)
!              whereas the transpose of A is used for calculations
!              If uplo_a is 'U' or 'L', the other triangle is not used at all,
!              i.e. it may contain arbitrary numbers
!
!  uplo_c      'U' if only the upper diagonal part of C is needed
!              'L' if only the upper diagonal part of C is needed
!              anything else if the full matrix C is needed
!              Please note: Even when uplo_c is 'U' or 'L', the other triangle may be
!              written to a certain extent, i.e. one shouldn't rely on the content there!
!
!  na          Number of rows/columns of A, number of rows of B and C
!
!  ncb         Number of columns  of B and C
!
!  a           Matrix A
!
!  lda         Leading dimension of a
!
!  b           Matrix B
!
!  ldb         Leading dimension of b
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!  c           Matrix C
!
!  ldc         Leading dimension of c
!
!-------------------------------------------------------------------------------

   implicit none

   character*1 uplo_a, uplo_c

   integer na, ncb, lda, ldb, nblk, mpi_comm_rows, mpi_comm_cols, ldc
   complex*16 a(lda,*), b(ldb,*), c(ldc,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_rows_np
   integer np, n, nb, nblk_mult, lrs, lre, lcs, lce
   integer gcol_min, gcol, goff
   integer nstor, nr_done, noff, np_bc, n_aux_bc, nvals
   integer, allocatable :: lrs_save(:), lre_save(:)

   logical a_lower, a_upper, c_lower, c_upper

   complex*16, allocatable:: aux_mat(:,:), aux_bc(:), tmp1(:,:), tmp2(:,:)


   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)

   l_rows = local_index(na,  my_prow, np_rows, nblk, -1) ! Local rows of a and b
   l_cols = local_index(ncb, my_pcol, np_cols, nblk, -1) ! Local cols of b

   ! Block factor for matrix multiplications, must be a multiple of nblk

   if(na/np_rows<=256) then
      nblk_mult = (31/nblk+1)*nblk
   else
      nblk_mult = (63/nblk+1)*nblk
   endif

   allocate(aux_mat(l_rows,nblk_mult))
   allocate(aux_bc(l_rows*nblk))
   allocate(lrs_save(nblk))
   allocate(lre_save(nblk))

   a_lower = .false.
   a_upper = .false.
   c_lower = .false.
   c_upper = .false.

   if(uplo_a=='u' .or. uplo_a=='U') a_upper = .true.
   if(uplo_a=='l' .or. uplo_a=='L') a_lower = .true.
   if(uplo_c=='u' .or. uplo_c=='U') c_upper = .true.
   if(uplo_c=='l' .or. uplo_c=='L') c_lower = .true.

   ! Build up the result matrix by processor rows

   do np = 0, np_rows-1

      ! In this turn, procs of row np assemble the result

      l_rows_np = local_index(na, np, np_rows, nblk, -1) ! local rows on receiving processors

      nr_done = 0 ! Number of rows done
      aux_mat = 0
      nstor = 0   ! Number of columns stored in aux_mat

      ! Loop over the blocks on row np

      do nb=0,(l_rows_np-1)/nblk

         goff  = nb*np_rows + np ! Global offset in blocks corresponding to nb

         ! Get the processor column which owns this block (A is transposed, so we need the column)
         ! and the offset in blocks within this column.
         ! The corresponding block column in A is then broadcast to all for multiplication with B

         np_bc = MOD(goff,np_cols)
         noff = goff/np_cols
         n_aux_bc = 0

         ! Gather up the complete block column of A on the owner

         do n = 1, min(l_rows_np-nb*nblk,nblk) ! Loop over columns to be broadcast

            gcol = goff*nblk + n ! global column corresponding to n
            if(nstor==0 .and. n==1) gcol_min = gcol

            lrs = 1       ! 1st local row number for broadcast
            lre = l_rows  ! last local row number for broadcast
            if(a_lower) lrs = local_index(gcol, my_prow, np_rows, nblk, +1)
            if(a_upper) lre = local_index(gcol, my_prow, np_rows, nblk, -1)

            if(lrs<=lre) then
               nvals = lre-lrs+1
               if(my_pcol == np_bc) aux_bc(n_aux_bc+1:n_aux_bc+nvals) = a(lrs:lre,noff*nblk+n)
               n_aux_bc = n_aux_bc + nvals
            endif

            lrs_save(n) = lrs
            lre_save(n) = lre

         enddo

         ! Broadcast block column

         call MPI_Bcast(aux_bc,n_aux_bc,MPI_DOUBLE_COMPLEX,np_bc,mpi_comm_cols,mpierr)

         ! Insert what we got in aux_mat

         n_aux_bc = 0
         do n = 1, min(l_rows_np-nb*nblk,nblk)
            nstor = nstor+1
            lrs = lrs_save(n)
            lre = lre_save(n)
            if(lrs<=lre) then
               nvals = lre-lrs+1
               aux_mat(lrs:lre,nstor) = aux_bc(n_aux_bc+1:n_aux_bc+nvals)
               n_aux_bc = n_aux_bc + nvals
            endif
         enddo

         ! If we got nblk_mult columns in aux_mat or this is the last block
         ! do the matrix multiplication

         if(nstor==nblk_mult .or. nb*nblk+nblk >= l_rows_np) then

            lrs = 1       ! 1st local row number for multiply
            lre = l_rows  ! last local row number for multiply
            if(a_lower) lrs = local_index(gcol_min, my_prow, np_rows, nblk, +1)
            if(a_upper) lre = local_index(gcol, my_prow, np_rows, nblk, -1)

            lcs = 1       ! 1st local col number for multiply
            lce = l_cols  ! last local col number for multiply
            if(c_upper) lcs = local_index(gcol_min, my_pcol, np_cols, nblk, +1)
            if(c_lower) lce = MIN(local_index(gcol, my_pcol, np_cols, nblk, -1),l_cols)

            if(lcs<=lce) then
               allocate(tmp1(nstor,lcs:lce),tmp2(nstor,lcs:lce))
               if(lrs<=lre) then
                  call zgemm('C','N',nstor,lce-lcs+1,lre-lrs+1,(1.d0,0.d0),aux_mat(lrs,1),ubound(aux_mat,1), &
                             b(lrs,lcs),ldb,(0.d0,0.d0),tmp1,nstor)
               else
                  tmp1 = 0
               endif

               ! Sum up the results and send to processor row np
               call mpi_reduce(tmp1,tmp2,nstor*(lce-lcs+1),MPI_DOUBLE_COMPLEX,MPI_SUM,np,mpi_comm_rows,mpierr)

               ! Put the result into C
               if(my_prow==np) c(nr_done+1:nr_done+nstor,lcs:lce) = tmp2(1:nstor,lcs:lce)

               deallocate(tmp1,tmp2)
            endif

            nr_done = nr_done+nstor
            nstor=0
            aux_mat(:,:)=0
         endif
      enddo
   enddo

   deallocate(aux_mat, aux_bc, lrs_save, lre_save)

end subroutine mult_ah_b_complex

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

subroutine solve_tridi( na, nev, d, e, q, ldq, nblk, mpi_comm_rows, mpi_comm_cols )

   implicit none

   integer  na, nev, ldq, nblk, mpi_comm_rows, mpi_comm_cols
   real*8 d(na), e(na), q(ldq,*)

   double precision sum1,sum2

   integer i, j, n, np, nc, nev1, l_cols, l_rows
   integer my_prow, my_pcol, np_rows, np_cols, mpierr

   integer, allocatable :: limits(:), l_col(:), p_col(:), l_col_bc(:), p_col_bc(:)


   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)

   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

   ! Set Q to 0

   q(1:l_rows, 1:l_cols) = 0.

   ! Get the limits of the subdivisons, each subdivison has as many cols
   ! as fit on the respective processor column

   allocate(limits(0:np_cols))

   limits(0) = 0
   do np=0,np_cols-1
      nc = local_index(na, np, np_cols, nblk, -1) ! number of columns on proc column np

      ! Check for the case that a column has have zero width.
      ! This is not supported!
      ! Scalapack supports it but delivers no results for these columns,
      ! which is rather annoying
      if(nc==0) then
         print *,'ERROR: Problem contains processor column with zero width'
         call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
      endif

      limits(np+1) = limits(np) + nc
   enddo

   ! Subdivide matrix by subtracting rank 1 modifications

   do i=1,np_cols-1
      n = limits(i)
      d(n) = d(n)-abs(e(n))
      d(n+1) = d(n+1)-abs(e(n))
   enddo

   ! Solve sub problems on processsor columns

   nc = limits(my_pcol) ! column after which my problem starts

   if(np_cols>1) then
      nev1 = l_cols ! all eigenvectors are needed
   else
      nev1 = MIN(nev,l_cols)
   endif
   call solve_tridi_col(l_cols, nev1, nc, d(nc+1), e(nc+1), q, ldq, nblk, mpi_comm_rows)

   ! If there is only 1 processor column, we are done

   if(np_cols==1) then
      deallocate(limits)
      return
   endif

   ! Set index arrays for Q columns

   ! Dense distribution scheme:

   allocate(l_col(na))
   allocate(p_col(na))

   n = 0
   do np=0,np_cols-1
      nc = local_index(na, np, np_cols, nblk, -1)
      do i=1,nc
         n = n+1
         l_col(n) = i
         p_col(n) = np
      enddo
   enddo

   ! Block cyclic distribution scheme, only nev columns are set:

   allocate(l_col_bc(na))
   allocate(p_col_bc(na))
   p_col_bc(:) = -1
   l_col_bc(:) = -1

   do i = 0, na-1, nblk*np_cols
      do j = 0, np_cols-1
         do n = 1, nblk
            if(i+j*nblk+n <= MIN(nev,na)) then
               p_col_bc(i+j*nblk+n) = j
               l_col_bc(i+j*nblk+n) = i/np_cols + n
            endif
         enddo
      enddo
   enddo

   ! Recursively merge sub problems

   call merge_recursive(0, np_cols)

   deallocate(limits,l_col,p_col,l_col_bc,p_col_bc)

contains
recursive subroutine merge_recursive(np_off, nprocs)

   implicit none

   ! noff is always a multiple of nblk_ev
   ! nlen-noff is always > nblk_ev

   integer np_off, nprocs
   integer np1, np2, noff, nlen, nmid, n, mpi_status(mpi_status_size)

   if(nprocs<=1) then
      ! Safety check only
      print *,"INTERNAL error merge_recursive: nprocs=",nprocs
      call mpi_abort(MPI_COMM_WORLD,1,mpierr)
   endif

   ! Split problem into 2 subproblems of size np1 / np2

   np1 = nprocs/2
   np2 = nprocs-np1

   if(np1 > 1) call merge_recursive(np_off, np1)
   if(np2 > 1) call merge_recursive(np_off+np1, np2)


   noff = limits(np_off)
   nmid = limits(np_off+np1) - noff
   nlen = limits(np_off+nprocs) - noff

   if(my_pcol==np_off) then
      do n=np_off+np1,np_off+nprocs-1
         call mpi_send(d(noff+1),nmid,MPI_REAL8,n,1,mpi_comm_cols,mpierr)
      enddo
   endif
   if(my_pcol>=np_off+np1 .and. my_pcol<np_off+nprocs) then
      call mpi_recv(d(noff+1),nmid,MPI_REAL8,np_off,1,mpi_comm_cols,mpi_status,mpierr)
   endif

   if(my_pcol==np_off+np1) then
      do n=np_off,np_off+np1-1
         call mpi_send(d(noff+nmid+1),nlen-nmid,MPI_REAL8,n,1,mpi_comm_cols,mpierr)
      enddo
   endif
   if(my_pcol>=np_off .and. my_pcol<np_off+np1) then
      call mpi_recv(d(noff+nmid+1),nlen-nmid,MPI_REAL8,np_off+np1,1,mpi_comm_cols,mpi_status,mpierr)
   endif

   if(nprocs == np_cols) then

      ! Last merge, result distribution must be block cyclic, noff==0,
      ! p_col_bc is set so that only nev eigenvalues are calculated

      call merge_systems(nlen, nmid, d(noff+1), e(noff+nmid), q, ldq, noff, &
                         nblk, mpi_comm_rows, mpi_comm_cols, l_col, p_col, &
                         l_col_bc, p_col_bc, np_off, nprocs )

   else

      ! Not last merge, leave dense column distribution

      call merge_systems(nlen, nmid, d(noff+1), e(noff+nmid), q, ldq, noff, &
                         nblk, mpi_comm_rows, mpi_comm_cols, l_col(noff+1), p_col(noff+1), &
                         l_col(noff+1), p_col(noff+1), np_off, nprocs )
   endif

end subroutine merge_recursive

end subroutine solve_tridi

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

subroutine solve_tridi_col( na, nev, nqoff, d, e, q, ldq, nblk, mpi_comm_rows )

   ! Solves the symmetric, tridiagonal eigenvalue problem on one processor column
   ! with the divide and conquer method.
   ! Works best if the number of processor rows is a power of 2!

   implicit none
   integer  na, nev, nqoff, ldq, nblk, mpi_comm_rows
   real*8 d(na), e(na), q(ldq,*)

   integer, parameter:: min_submatrix_size = 16 ! Minimum size of the submatrices to be used

   real*8, allocatable :: qmat1(:,:), qmat2(:,:)

   real*8 sum1,sum2

   integer i, j, n, np
   integer ndiv, noff, nmid, nlen, max_size
   integer my_prow, np_rows, mpierr

   integer, allocatable :: limits(:), l_col(:), p_col_i(:), p_col_o(:)

   call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
   call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)

   ! Calculate the number of subdivisions needed.

   n = na
   ndiv = 1
   do while(2*ndiv<=np_rows .and. n>2*min_submatrix_size)
      n = ((n+3)/4)*2 ! the bigger one of the two halves, we want EVEN boundaries
      ndiv = ndiv*2
   enddo

   ! If there is only 1 processor row and not all eigenvectors are needed
   ! and the matrix size is big enough, then use 2 subdivisions
   ! so that merge_systems is called once and only the needed
   ! eigenvectors are calculated for the final problem.

   if(np_rows==1 .and. nev<na .and. na>2*min_submatrix_size) ndiv = 2

   allocate(limits(0:ndiv))

   limits(0) = 0
   limits(ndiv) = na

   n = ndiv
   do while(n>1)
      n = n/2 ! n is always a power of 2
      do i=0,ndiv-1,2*n
         ! We want to have even boundaries (for cache line alignments)
         limits(i+n) = limits(i) + ((limits(i+2*n)-limits(i)+3)/4)*2
      enddo
   enddo

   ! Calculate the maximum size of a subproblem

   max_size = 0
   do i=1,ndiv
      max_size = MAX(max_size,limits(i)-limits(i-1))
   enddo

   ! Subdivide matrix by subtracting rank 1 modifications

   do i=1,ndiv-1
      n = limits(i)
      d(n) = d(n)-abs(e(n))
      d(n+1) = d(n+1)-abs(e(n))
   enddo

   if(np_rows==1)    then

      ! For 1 processor row there may be 1 or 2 subdivisions

      do n=0,ndiv-1
         noff = limits(n)        ! Start of subproblem
         nlen = limits(n+1)-noff ! Size of subproblem

         call solve_tridi_single(nlen,d(noff+1),e(noff+1),q(nqoff+noff+1,noff+1),ubound(q,1))

      enddo

   else

      ! Solve sub problems in parallel with solve_tridi_single
      ! There is at maximum 1 subproblem per processor

      allocate(qmat1(max_size,max_size))
      allocate(qmat2(max_size,max_size))

      qmat1 = 0 ! Make sure that all elements are defined

      if(my_prow < ndiv) then

         noff = limits(my_prow)        ! Start of subproblem
         nlen = limits(my_prow+1)-noff ! Size of subproblem

         call solve_tridi_single(nlen,d(noff+1),e(noff+1),qmat1,ubound(qmat1,1))
      endif

      ! Fill eigenvectors in qmat1 into global matrix q

      do np = 0, ndiv-1

         noff = limits(np)
         nlen = limits(np+1)-noff

         call MPI_Bcast(d(noff+1),nlen,MPI_REAL8,np,mpi_comm_rows,mpierr)
         qmat2 = qmat1
         call MPI_Bcast(qmat2,max_size*max_size,MPI_REAL8,np,mpi_comm_rows,mpierr)

         do i=1,nlen
            call distribute_global_column(qmat2(1,i), q(1,noff+i), nqoff+noff, nlen, my_prow, np_rows, nblk)
         enddo
      enddo

      deallocate(qmat1, qmat2)

   endif

   ! Allocate and set index arrays l_col and p_col

   allocate(l_col(na), p_col_i(na),  p_col_o(na))

   do i=1,na
      l_col(i) = i
      p_col_i(i) = 0
      p_col_o(i) = 0
   enddo

   ! Merge subproblems

   n = 1
   do while(n<ndiv) ! if ndiv==1, the problem was solved by single call to solve_tridi_single

      do i=0,ndiv-1,2*n

         noff = limits(i)
         nmid = limits(i+n) - noff
         nlen = limits(i+2*n) - noff

         if(nlen == na) then
           ! Last merge, set p_col_o=-1 for unneeded (output) eigenvectors
           p_col_o(nev+1:na) = -1
         endif

         call merge_systems(nlen, nmid, d(noff+1), e(noff+nmid), q, ldq, nqoff+noff, nblk, &
                            mpi_comm_rows, mpi_comm_self, l_col(noff+1), p_col_i(noff+1), &
                            l_col(noff+1), p_col_o(noff+1), 0, 1)

      enddo

      n = 2*n

   enddo

   deallocate(limits, l_col, p_col_i, p_col_o)

end subroutine solve_tridi_col

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

subroutine solve_tridi_single(nlen, d, e, q, ldq)

   ! Solves the symmetric, tridiagonal eigenvalue problem on a single processor.

   implicit none

   integer nlen, ldq
   real*8 d(nlen), e(nlen), q(ldq,nlen)

   double precision sum1,sum2

   real*8, allocatable :: work(:), qtmp(:), w(:)
   real*8 dtmp,vl,vu,abstol

   integer i, j, lwork, info, mpierr, m, ldz, il, iu
   integer, allocatable :: iwork(:), ifail(:)

   lwork = 1 + 4*nlen + nlen**2
   allocate(work(lwork))

   call dsteqr('I',nlen,d,e,q,ldq,work,info)

   ! If DSTEQR fails also, we don't know what to do further ...
   if(info /= 0) then
      print '(a,i8,a)','ERROR: Lapack routine DSTEQR failed, info= ',info,', Aborting!'
      call mpi_abort(mpi_comm_world,0,mpierr)
   endif

   deallocate(work)

end subroutine solve_tridi_single

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


subroutine merge_systems( na, nm, d, e, q, ldq, nqoff, nblk, mpi_comm_rows, mpi_comm_cols, &
                          l_col, p_col, l_col_out, p_col_out, npc_0, npc_n)

   implicit none

   integer  na, nm, ldq, nqoff, nblk, mpi_comm_rows, mpi_comm_cols, npc_0, npc_n
   integer  l_col(na), p_col(na), l_col_out(na), p_col_out(na)
   real*8 d(na), e, q(ldq,*)

   integer, parameter :: max_strip=128

   real*8 beta, sig, s, c, t, tau, rho, eps, tol, dlamch, dlapy2, qtrans(2,2), dmax, zmax, d1new, d2new, sum1
   real*8 z(na), d1(na), d2(na), z1(na), delta(na), dbase(na), ddiff(na), ev_scale(na), tmp(na)
   real*8 d1u(na), zu(na), d1l(na), zl(na)
   real*8, allocatable :: qtmp1(:,:), qtmp2(:,:), ev(:,:)

   integer i, j, na1, na2, l_rows, l_cols, l_rqs, l_rqe, l_rqm, ns, info
   integer l_rnm, nnzu, nnzl, ndef, ncnt, max_local_cols, l_cols_qreorg, np, l_idx, nqcols1, nqcols2
   integer my_proc, n_procs, my_prow, my_pcol, np_rows, np_cols, mpierr, mpi_status(mpi_status_size)
   integer np_next, np_prev, np_rem
   integer idx(na), idx1(na), idx2(na)
   integer coltyp(na), idxq1(na), idxq2(na)

   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)

   ! If my processor column isn't in the requested set, do nothing

   if(my_pcol<npc_0 .or. my_pcol>=npc_0+npc_n) return

   ! Determine number of "next" and "prev" column for ring sends

   if(my_pcol == npc_0+npc_n-1) then
      np_next = npc_0
   else
      np_next = my_pcol + 1
   endif

   if(my_pcol == npc_0) then
      np_prev = npc_0+npc_n-1
   else
      np_prev = my_pcol - 1
   endif

   call check_monotony(nm,d,'Input1')
   call check_monotony(na-nm,d(nm+1),'Input2')

   ! Get global number of processors and my processor number.
   ! Please note that my_proc does not need to match any real processor number,
   ! it is just used for load balancing some loops.

   n_procs = np_rows*npc_n
   my_proc = my_prow*npc_n + (my_pcol-npc_0) ! Row major

   ! Local limits of the rows of Q

   l_rqs = local_index(nqoff+1 , my_prow, np_rows, nblk, +1) ! First row of Q
   l_rqm = local_index(nqoff+nm, my_prow, np_rows, nblk, -1) ! Last row <= nm
   l_rqe = local_index(nqoff+na, my_prow, np_rows, nblk, -1) ! Last row of Q

   l_rnm  = l_rqm-l_rqs+1 ! Number of local rows <= nm
   l_rows = l_rqe-l_rqs+1 ! Total number of local rows

   ! My number of local columns

   l_cols = COUNT(p_col(1:na)==my_pcol)

   ! Get max number of local columns

   max_local_cols = 0
   do np = npc_0, npc_0+npc_n-1
      max_local_cols = MAX(max_local_cols,COUNT(p_col(1:na)==np))
   enddo

   ! Calculations start here

   beta = abs(e)
   sig  = sign(1.d0,e)

   ! Calculate rank-1 modifier z

   z(:) = 0

   if(MOD((nqoff+nm-1)/nblk,np_rows)==my_prow) then
      ! nm is local on my row
      do i = 1, na
         if(p_col(i)==my_pcol) z(i) = q(l_rqm,l_col(i))
      enddo
   endif

   if(MOD((nqoff+nm)/nblk,np_rows)==my_prow) then
      ! nm+1 is local on my row
      do i = 1, na
         if(p_col(i)==my_pcol) z(i) = z(i) + sig*q(l_rqm+1,l_col(i))
      enddo
   endif

   call global_gather(z, na)

   ! Normalize z so that norm(z) = 1.  Since z is the concatenation of
   ! two normalized vectors, norm2(z) = sqrt(2).

   z = z/sqrt(2.0d0)
   rho = 2.*beta

   ! Calculate index for merging both systems by ascending eigenvalues

   call DLAMRG( nm, na-nm, d, 1, 1, idx )

   ! Calculate the allowable deflation tolerance

   zmax = maxval(abs(z))
   dmax = maxval(abs(d))
   EPS = DLAMCH( 'Epsilon' )
   TOL = 8.*EPS*MAX(dmax,zmax)

   ! If the rank-1 modifier is small enough, no more needs to be done
   ! except to reorganize D and Q

   IF( RHO*zmax <= TOL ) THEN

      ! Rearrange eigenvalues

      tmp = d
      do i=1,na
         d(i) = tmp(idx(i))
      enddo

      ! Rearrange eigenvectors

      call resort_ev(idx)
      return
   ENDIF

   ! Merge and deflate system

   na1 = 0
   na2 = 0

   ! COLTYP:
   ! 1 : non-zero in the upper half only;
   ! 2 : dense;
   ! 3 : non-zero in the lower half only;
   ! 4 : deflated.

   coltyp(1:nm) = 1
   coltyp(nm+1:na) = 3

   do i=1,na

      if(rho*abs(z(idx(i))) <= tol) then

         ! Deflate due to small z component.

         na2 = na2+1
         d2(na2)   = d(idx(i))
         idx2(na2) = idx(i)
         coltyp(idx(i)) = 4

      else if(na1>0) then

         ! Check if eigenvalues are close enough to allow deflation.

         S = Z(idx(i))
         C = Z1(na1)

         ! Find sqrt(a**2+b**2) without overflow or
         ! destructive underflow.

         TAU = DLAPY2( C, S )
         T = D1(na1) - D(idx(i))
         C = C / TAU
         S = -S / TAU
         IF( ABS( T*C*S ) <= TOL ) THEN

            ! Deflation is possible.

            na2 = na2+1

            Z1(na1) = TAU

            d2new = D(idx(i))*C**2 + D1(na1)*S**2
            d1new = D(idx(i))*S**2 + D1(na1)*C**2

            ! D(idx(i)) >= D1(na1) and C**2 + S**2 == 1.0
            ! This means that after the above transformation it must be
            !    D1(na1) <= d1new <= D(idx(i))
            !    D1(na1) <= d2new <= D(idx(i))
            !
            ! D1(na1) may get bigger but it is still smaller than the next D(idx(i+1))
            ! so there is no problem with sorting here.
            ! d2new <= D(idx(i)) which means that it might be smaller than D2(na2-1)
            ! which makes a check (and possibly a resort) necessary.
            !
            ! The above relations may not hold exactly due to numeric differences
            ! so they have to be enforced in order not to get troubles with sorting.


            if(d1new<D1(na1)  ) d1new = D1(na1)
            if(d1new>D(idx(i))) d1new = D(idx(i))
  
          if(d2new<D1(na1)  ) d2new = D1(na1)
            if(d2new>D(idx(i))) d2new = D(idx(i))

            D1(na1) = d1new

            do j=na2-1,1,-1
               if(d2new<d2(j)) then
                  d2(j+1)   = d2(j)
                  idx2(j+1) = idx2(j)
               else
                  exit ! Loop
               endif
            enddo

            d2(j+1)   = d2new
            idx2(j+1) = idx(i)

            qtrans(1,1) = C; qtrans(1,2) =-S
            qtrans(2,1) = S; qtrans(2,2) = C

            call transform_columns(idx(i), idx1(na1))

            if(coltyp(idx(i))==1 .and. coltyp(idx1(na1))/=1) coltyp(idx1(na1)) = 2
            if(coltyp(idx(i))==3 .and. coltyp(idx1(na1))/=3) coltyp(idx1(na1)) = 2

            coltyp(idx(i)) = 4

         else
            na1 = na1+1
            d1(na1) = d(idx(i))
            z1(na1) = z(idx(i))
            idx1(na1) = idx(i)
         endif
      else
         na1 = na1+1
         d1(na1) = d(idx(i))
         z1(na1) = z(idx(i))
         idx1(na1) = idx(i)
      endif

   enddo
   call check_monotony(na1,d1,'Sorted1')
   call check_monotony(na2,d2,'Sorted2')

   if(na1==1 .or. na1==2) then
      ! if(my_proc==0) print *,'--- Remark solve_tridi: na1==',na1,' proc==',myid

      if(na1==1) then
         d(1) = d1(1) + rho*z1(1)**2 ! solve secular equation
      else ! na1==2
         call DLAED5(1, d1, z1, qtrans(1,1), rho, d(1))
         call DLAED5(2, d1, z1, qtrans(1,2), rho, d(2))
         call transform_columns(idx1(1), idx1(2))
      endif

      ! Add the deflated eigenvalues
      d(na1+1:na) = d2(1:na2)

      ! Calculate arrangement of all eigenvalues  in output

      call DLAMRG( na1, na-na1, d, 1, 1, idx )

      ! Rearrange eigenvalues

      tmp = d
      do i=1,na
         d(i) = tmp(idx(i))
      enddo

      ! Rearrange eigenvectors

      do i=1,na
         if(idx(i)<=na1) then
            idxq1(i) = idx1(idx(i))
         else
            idxq1(i) = idx2(idx(i)-na1)
         endif
      enddo

      call resort_ev(idxq1)

   else if(na1>2) then

      ! Solve secular equation

      z(1:na1) = 1
      dbase(1:na1) = 0
      ddiff(1:na1) = 0
      info = 0

      DO i = my_proc+1, na1, n_procs ! work distributed over all processors

         ! DLAED4 is unstable, thus we use bisection
         call solve_secular_equation(na1, i, d1, z1, delta, rho, s)

         ! Compute updated z

         do j=1,na1
            if(i/=j)  z(j) = z(j)*( delta(j) / (d1(j)-d1(i)) )
         enddo
         z(i) = z(i)*delta(i)

         ! store dbase/ddiff

         if(i<na1) then
            if(abs(delta(i+1)) < abs(delta(i))) then
               dbase(i) = d1(i+1)
               ddiff(i) = delta(i+1)
            else
               dbase(i) = d1(i)
               ddiff(i) = delta(i)
            endif
         else
            dbase(i) = d1(i)
            ddiff(i) = delta(i)
         endif

      enddo

      call global_product(z, na1)
      z(1:na1) = SIGN( SQRT( -z(1:na1) ), z1(1:na1) )

      call global_gather(dbase, na1)
      call global_gather(ddiff, na1)

      d(1:na1) = dbase(1:na1) - ddiff(1:na1)

      ! Calculate scale factors for eigenvectors

      ev_scale(:) = 0

      DO i = my_proc+1, na1, n_procs ! work distributed over all processors

         ! tmp(1:na1) = z(1:na1) / delta(1:na1,i)  ! original code
         ! tmp(1:na1) = z(1:na1) / (d1(1:na1)-d(i))! bad results

         ! All we want to calculate is tmp = (d1(1:na1)-dbase(i))+ddiff(i)
         ! in exactly this order, but we want to prevent compiler optimization

         tmp(1:na1) = d1(1:na1)-dbase(i)
         call v_add_s(tmp,na1,ddiff(i))
         tmp(1:na1) = z(1:na1) / (tmp(1:na1)+1.0e-26)
         ev_scale(i) = 1.0/(sqrt(dot_product(tmp(1:na1),tmp(1:na1)))+1.0e-26)

      enddo
      call global_gather(ev_scale, na1)

      ! Add the deflated eigenvalues
      d(na1+1:na) = d2(1:na2)

      ! Calculate arrangement of all eigenvalues  in output

      call DLAMRG( na1, na-na1, d, 1, 1, idx )

      ! Rearrange eigenvalues

      tmp = d
      do i=1,na
         d(i) = tmp(idx(i))
      enddo

      call check_monotony(na,d,'Output')

      ! Eigenvector calculations


      ! Calculate the number of columns in the new local matrix Q
      ! which are updated from non-deflated/deflated eigenvectors.
      ! idxq1/2 stores the global column numbers.

      nqcols1 = 0 ! number of non-deflated eigenvectors
      nqcols2 = 0 ! number of deflated eigenvectors
      DO i = 1, na
         if(p_col_out(i)==my_pcol) then
            if(idx(i)<=na1) then
               nqcols1 = nqcols1+1
               idxq1(nqcols1) = i
            else
               nqcols2 = nqcols2+1
               idxq2(nqcols2) = i
            endif
         endif
      enddo

      allocate(ev(max_local_cols,MIN(max_strip,MAX(1,nqcols1))))
      allocate(qtmp1(MAX(1,l_rows),max_local_cols))
      allocate(qtmp2(MAX(1,l_rows),MIN(max_strip,MAX(1,nqcols1))))

      ! Gather nonzero upper/lower components of old matrix Q
      ! which are needed for multiplication with new eigenvectors

      qtmp1 = 0 ! May contain empty (unset) parts
      qtmp2 = 0 ! Not really needed

      nnzu = 0
      nnzl = 0
      do i = 1, na1
         l_idx = l_col(idx1(i))
         if(p_col(idx1(i))==my_pcol) then
            if(coltyp(idx1(i))==1 .or. coltyp(idx1(i))==2) then
               nnzu = nnzu+1
               qtmp1(1:l_rnm,nnzu) = q(l_rqs:l_rqm,l_idx)
            endif
            if(coltyp(idx1(i))==3 .or. coltyp(idx1(i))==2) then
               nnzl = nnzl+1
               qtmp1(l_rnm+1:l_rows,nnzl) = q(l_rqm+1:l_rqe,l_idx)
            endif
         endif
      enddo

      ! Gather deflated eigenvalues behind nonzero components

      ndef = max(nnzu,nnzl)
      do i = 1, na2
         l_idx = l_col(idx2(i))
         if(p_col(idx2(i))==my_pcol) then
            ndef = ndef+1
            qtmp1(1:l_rows,ndef) = q(l_rqs:l_rqe,l_idx)
         endif
      enddo

      l_cols_qreorg = ndef ! Number of columns in reorganized matrix

      ! Set (output) Q to 0, it will sum up new Q

      DO i = 1, na
         if(p_col_out(i)==my_pcol) q(l_rqs:l_rqe,l_col_out(i)) = 0
      enddo

      np_rem = my_pcol

      do np = 1, npc_n

         ! Do a ring send of qtmp1

         if(np>1) then

            if(np_rem==npc_0) then
               np_rem = npc_0+npc_n-1
            else
               np_rem = np_rem-1
            endif

            call MPI_Sendrecv_replace(qtmp1, l_rows*max_local_cols, MPI_REAL8, &
                                      np_next, 1111, np_prev, 1111, &
                                      mpi_comm_cols, mpi_status, mpierr)
         endif

         ! Gather the parts in d1 and z which are fitting to qtmp1.
         ! This also delivers nnzu/nnzl for proc np_rem

         nnzu = 0
         nnzl = 0
         do i=1,na1
            if(p_col(idx1(i))==np_rem) then
               if(coltyp(idx1(i))==1 .or. coltyp(idx1(i))==2) then
                  nnzu = nnzu+1
                  d1u(nnzu) = d1(i)
                  zu (nnzu) = z (i)
               endif
               if(coltyp(idx1(i))==3 .or. coltyp(idx1(i))==2) then
                  nnzl = nnzl+1
                  d1l(nnzl) = d1(i)
                  zl (nnzl) = z (i)
               endif
            endif
         enddo

         ! Set the deflated eigenvectors in Q (comming from proc np_rem)

         ndef = MAX(nnzu,nnzl) ! Remote counter in input matrix
         do i = 1, na
            j = idx(i)
            if(j>na1) then
               if(p_col(idx2(j-na1))==np_rem) then
                  ndef = ndef+1
                  if(p_col_out(i)==my_pcol) &
                     q(l_rqs:l_rqe,l_col_out(i)) = qtmp1(1:l_rows,ndef)
               endif
            endif
         enddo

         do ns = 0, nqcols1-1, max_strip ! strimining loop

            ncnt = MIN(max_strip,nqcols1-ns) ! number of columns in this strip

            ! Get partial result from (output) Q

            do i = 1, ncnt
               qtmp2(1:l_rows,i) = q(l_rqs:l_rqe,l_col_out(idxq1(i+ns)))
            enddo

            ! Compute eigenvectors of the rank-1 modified matrix.
            ! Parts for multiplying with upper half of Q:

            do i = 1, ncnt
               j = idx(idxq1(i+ns))
               ! Calculate the j-th eigenvector of the deflated system
               ! See above why we are doing it this way!
               tmp(1:nnzu) = d1u(1:nnzu)-dbase(j)
               call v_add_s(tmp,nnzu,ddiff(j))
               ev(1:nnzu,i) = zu(1:nnzu) / (tmp(1:nnzu)+1.0e-26) * ev_scale(j)
            enddo

            ! Multiply old Q with eigenvectors (upper half)

            if(l_rnm>0 .and. ncnt>0 .and. nnzu>0) &
               call dgemm('N','N',l_rnm,ncnt,nnzu,1.d0,qtmp1,ubound(qtmp1,1),ev,ubound(ev,1), &
                          1.d0,qtmp2(1,1),ubound(qtmp2,1))

            ! Compute eigenvectors of the rank-1 modified matrix.
            ! Parts for multiplying with lower half of Q:

            do i = 1, ncnt
               j = idx(idxq1(i+ns))
               ! Calculate the j-th eigenvector of the deflated system
               ! See above why we are doing it this way!
               tmp(1:nnzl) = d1l(1:nnzl)-dbase(j)
               call v_add_s(tmp,nnzl,ddiff(j))
               ev(1:nnzl,i) = zl(1:nnzl) / (tmp(1:nnzl)+1.0e-26) * ev_scale(j)
            enddo

            ! Multiply old Q with eigenvectors (lower half)

            if(l_rows-l_rnm>0 .and. ncnt>0 .and. nnzl>0) &
               call dgemm('N','N',l_rows-l_rnm,ncnt,nnzl,1.d0,qtmp1(l_rnm+1,1),ubound(qtmp1,1),ev,ubound(ev,1), &
                          1.d0,qtmp2(l_rnm+1,1),ubound(qtmp2,1))

            ! Put partial result into (output) Q

            do i = 1, ncnt
               q(l_rqs:l_rqe,l_col_out(idxq1(i+ns))) = qtmp2(1:l_rows,i)
            enddo
         enddo
      enddo

      deallocate(ev, qtmp1, qtmp2)

   endif

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

contains
subroutine resort_ev(idx_ev)

   implicit none

   integer idx_ev(*)
   integer i, nc, pc1, pc2, lc1, lc2, l_cols_out

   real*8, allocatable :: qtmp(:,:)

   if(l_rows==0) return ! My processor column has no work to do

   ! Resorts eigenvectors so that q_new(:,i) = q_old(:,idx_ev(i))

   l_cols_out = COUNT(p_col_out(1:na)==my_pcol)
   allocate(qtmp(l_rows,l_cols_out))


   nc = 0

   do i=1,na

      pc1 = p_col(idx_ev(i))
      lc1 = l_col(idx_ev(i))
      pc2 = p_col_out(i)

      if(pc2<0) cycle ! This column is not needed in output

      if(pc2==my_pcol) nc = nc+1 ! Counter for output columns

      if(pc1==my_pcol) then
         if(pc2==my_pcol) then
            ! send and recieve column are local
            qtmp(1:l_rows,nc) = q(l_rqs:l_rqe,lc1)
         else
            call mpi_send(q(l_rqs,lc1),l_rows,MPI_REAL8,pc2,mod(i,4096),mpi_comm_cols,mpierr)
         endif
      else if(pc2==my_pcol) then
         call mpi_recv(qtmp(1,nc),l_rows,MPI_REAL8,pc1,mod(i,4096),mpi_comm_cols,mpi_status,mpierr)
      endif
   enddo

   ! Insert qtmp into (output) q

   nc = 0

   do i=1,na

      pc2 = p_col_out(i)
      lc2 = l_col_out(i)

      if(pc2==my_pcol) then
         nc = nc+1
         q(l_rqs:l_rqe,lc2) = qtmp(1:l_rows,nc)
      endif
   enddo

   deallocate(qtmp)

end subroutine resort_ev

subroutine transform_columns(col1, col2)

   implicit none

   integer col1, col2
   integer pc1, pc2, lc1, lc2

   if(l_rows==0) return ! My processor column has no work to do

   pc1 = p_col(col1)
   lc1 = l_col(col1)
   pc2 = p_col(col2)
   lc2 = l_col(col2)

   if(pc1==my_pcol) then
      if(pc2==my_pcol) then
         ! both columns are local
         tmp(1:l_rows)      = q(l_rqs:l_rqe,lc1)*qtrans(1,1) + q(l_rqs:l_rqe,lc2)*qtrans(2,1)
         q(l_rqs:l_rqe,lc2) = q(l_rqs:l_rqe,lc1)*qtrans(1,2) + q(l_rqs:l_rqe,lc2)*qtrans(2,2)
         q(l_rqs:l_rqe,lc1) = tmp(1:l_rows)
      else
         call mpi_sendrecv(q(l_rqs,lc1),l_rows,MPI_REAL8,pc2,1, &
                           tmp,l_rows,MPI_REAL8,pc2,1, &
                           mpi_comm_cols,mpi_status,mpierr)
         q(l_rqs:l_rqe,lc1) = q(l_rqs:l_rqe,lc1)*qtrans(1,1) + tmp(1:l_rows)*qtrans(2,1)
      endif
   else if(pc2==my_pcol) then
      call mpi_sendrecv(q(l_rqs,lc2),l_rows,MPI_REAL8,pc1,1, &
                        tmp,l_rows,MPI_REAL8,pc1,1, &
                        mpi_comm_cols,mpi_status,mpierr)
      q(l_rqs:l_rqe,lc2) = tmp(1:l_rows)*qtrans(1,2) + q(l_rqs:l_rqe,lc2)*qtrans(2,2)
   endif

end subroutine transform_columns

subroutine global_gather(z, n)

   ! This routine sums up z over all processors.
   ! It should only be used for gathering distributed results,
   ! i.e. z(i) should be nonzero on exactly 1 processor column,
   ! otherways the results may be numerically different on different columns

   implicit none

   integer n
   real*8 z(n)

   real*8 tmp(n)

   if(npc_n==1 .and. np_rows==1) return ! nothing to do

   ! Do an mpi_allreduce over processor rows

   call mpi_allreduce(z, tmp, n, MPI_REAL8, MPI_SUM, mpi_comm_rows, mpierr)

   ! If only 1 processor column, we are done
   if(npc_n==1) then
      z(:) = tmp(:)
      return
   endif

   ! If all processor columns are involved, we can use mpi_allreduce
   if(npc_n==np_cols) then
      call mpi_allreduce(tmp, z, n, MPI_REAL8, MPI_SUM, mpi_comm_cols, mpierr)
      return
   endif

   ! Do a ring send over processor columns
   z(:) = 0
   do np = 1, npc_n
      z(:) = z(:) + tmp(:)
      call MPI_Sendrecv_replace(z, n, MPI_REAL8, np_next, 1111, np_prev, 1111, &
                                mpi_comm_cols, mpi_status, mpierr)
   enddo

end subroutine global_gather

subroutine global_product(z, n)

   ! This routine calculates the global product of z.

   implicit none

   integer n
   real*8 z(n)

   real*8 tmp(n)

   if(npc_n==1 .and. np_rows==1) return ! nothing to do

   ! Do an mpi_allreduce over processor rows

   call mpi_allreduce(z, tmp, n, MPI_REAL8, MPI_PROD, mpi_comm_rows, mpierr)

   ! If only 1 processor column, we are done
   if(npc_n==1) then
      z(:) = tmp(:)
      return
   endif

   ! If all processor columns are involved, we can use mpi_allreduce
   if(npc_n==np_cols) then
      call mpi_allreduce(tmp, z, n, MPI_REAL8, MPI_PROD, mpi_comm_cols, mpierr)
      return
   endif

   ! We send all vectors to the first proc, do the product there
   ! and redistribute the result.

   if(my_pcol == npc_0) then
      z(1:n) = tmp(1:n)
      do np = npc_0+1, npc_0+npc_n-1
         call mpi_recv(tmp,n,MPI_REAL8,np,1111,mpi_comm_cols,mpi_status,mpierr)
         z(1:n) = z(1:n)*tmp(1:n)
      enddo
      do np = npc_0+1, npc_0+npc_n-1
         call mpi_send(z,n,MPI_REAL8,np,1111,mpi_comm_cols,mpierr)
      enddo
   else
      call mpi_send(tmp,n,MPI_REAL8,npc_0,1111,mpi_comm_cols,mpierr)
      call mpi_recv(z  ,n,MPI_REAL8,npc_0,1111,mpi_comm_cols,mpi_status,mpierr)
   endif

end subroutine global_product

subroutine check_monotony(n,d,text)

! This is a test routine for checking if the eigenvalues are monotonically increasing.
! It is for debug purposes only, an error should never be triggered!

   implicit none

   integer n
   real*8 d(n)
   character*(*) text

   integer i

   do i=1,n-1
      if(d(i+1)<d(i)) then
         print '(a,a,i8,2g25.17)','Monotony error on ',text,i,d(i),d(i+1)
         call mpi_abort(mpi_comm_world,0,mpierr)
      endif
   enddo

end subroutine check_monotony

end subroutine merge_systems

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

subroutine v_add_s(v,n,s)
   implicit none
   integer n
   real*8 v(n),s

   v(:) = v(:) + s
end subroutine v_add_s

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

subroutine distribute_global_column(g_col, l_col, noff, nlen, my_prow, np_rows, nblk)

   implicit none

   real*8 g_col(nlen), l_col(*)
   integer noff, nlen, my_prow, np_rows, nblk

   integer nbs, nbe, jb, g_off, l_off, js, je

   nbs = noff/(nblk*np_rows)
   nbe = (noff+nlen-1)/(nblk*np_rows)

   do jb = nbs, nbe

      g_off = jb*nblk*np_rows + nblk*my_prow
      l_off = jb*nblk

      js = MAX(noff+1-g_off,1)
      je = MIN(noff+nlen-g_off,nblk)

      if(je<js) cycle

      l_col(l_off+js:l_off+je) = g_col(g_off+js-noff:g_off+je-noff)

  enddo

end subroutine distribute_global_column

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

subroutine solve_secular_equation(n, i, d, z, delta, rho, dlam)

!-------------------------------------------------------------------------------
! This routine solves the secular equation of a symmetric rank 1 modified
! diagonal matrix:
!
!    1. + rho*SUM(z(:)**2/(d(:)-x)) = 0
!
! It does the same as the LAPACK routine DLAED4 but it uses a bisection technique
! which is more robust (it always yields a solution) but also slower
! than the algorithm used in DLAED4.
!
! The same restictions than in DLAED4 hold, namely:
!
!   rho > 0   and   d(i+1) > d(i)
!
! but this routine will not terminate with error if these are not satisfied
! (it will normally converge to a pole in this case).
!
! The output in DELTA(j) is always (D(j) - lambda_I), even for the cases
! N=1 and N=2 which is not compatible with DLAED4.
! Thus this routine shouldn't be used for these cases as a simple replacement
! of DLAED4.
!
! The arguments are the same as in DLAED4 (with the exception of the INFO argument):
!
!
!  N      (input) INTEGER
!         The length of all arrays.
!
!  I      (input) INTEGER
!         The index of the eigenvalue to be computed.  1 <= I <= N.
!
!  D      (input) DOUBLE PRECISION array, dimension (N)
!         The original eigenvalues.  It is assumed that they are in
!         order, D(I) < D(J)  for I < J.
!
!  Z      (input) DOUBLE PRECISION array, dimension (N)
!         The components of the updating vector.
!
!  DELTA  (output) DOUBLE PRECISION array, dimension (N)
!         DELTA contains (D(j) - lambda_I) in its  j-th component.
!         See remark above about DLAED4 compatibility!
!
!  RHO    (input) DOUBLE PRECISION
!         The scalar in the symmetric updating formula.
!
!  DLAM   (output) DOUBLE PRECISION
!         The computed lambda_I, the I-th updated eigenvalue.
!-------------------------------------------------------------------------------


   implicit none

   integer n, i
   real*8 d(n), z(n), delta(n), rho, dlam

   integer iter
   real*8 a, b, x, y, dshift

   ! In order to obtain sufficient numerical accuracy we have to shift the problem
   ! either by d(i) or d(i+1), whichever is closer to the solution

   ! Upper and lower bound of the shifted solution interval are a and b


   if(i==n) then

      ! Special case: Last eigenvalue
      ! We shift always by d(n), lower bound is d(n),
      ! upper bound is determined by a guess:

      dshift = d(n)
      delta(:) = d(:) - dshift

      a = 0. ! delta(n)
      b = rho*SUM(z(:)**2) + 1. ! rho*SUM(z(:)**2) is the lower bound for the guess

   else

      ! Other eigenvalues: lower bound is d(i), upper bound is d(i+1)
      ! We check the sign of the function in the midpoint of the interval
      ! in order to determine if eigenvalue is more close to d(i) or d(i+1)

      x = 0.5*(d(i)+d(i+1))
      y = 1. + rho*SUM(z(:)**2/(d(:)-x))

      if(y>0) then
         ! solution is next to d(i)
         dshift = d(i)
      else
         ! solution is next to d(i+1)
         dshift = d(i+1)
      endif

      delta(:) = d(:) - dshift
      a = delta(i)
      b = delta(i+1)

   endif

   ! Bisection:

   do iter=1,200

      ! Interval subdivision

      x = 0.5*(a+b)

      if(x==a .or. x==b) exit   ! No further interval subdivisions possible
      if(abs(x) < 1.d-200) exit ! x next to pole

      ! evaluate value at x

      y = 1. + rho*SUM(z(:)**2/(delta(:)-x))

      if(y==0) then
         ! found exact solution
         exit
      elseif(y>0) then
         b = x
      else
         a = x
      endif

   enddo

   ! Solution:

   dlam = x + dshift
   delta(:) = delta(:) - x

end subroutine solve_secular_equation

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

integer function local_index(idx, my_proc, num_procs, nblk, iflag)

!-------------------------------------------------------------------------------
!  local_index: returns the local index for a given global index
!               If the global index has no local index on the
!               processor my_proc behaviour is defined by iflag
!
!  Parameters
!
!  idx         Global index
!
!  my_proc     Processor row/column for which to calculate the local index
!
!  num_procs   Total number of processors along row/column
!
!  nblk        Blocksize
!
!  iflag       Controls the behaviour if idx is not on local processor
!              iflag< 0 : Return last local index before that row/col
!              iflag==0 : Return 0
!              iflag> 0 : Return next local index after that row/col
!-------------------------------------------------------------------------------

   implicit none

   integer idx, my_proc, num_procs, nblk, iflag

   integer iblk

   iblk = (idx-1)/nblk  ! global block number, 0 based

   if(mod(iblk,num_procs) == my_proc) then

      ! block is local, always return local row/col number

      local_index = (iblk/num_procs)*nblk + mod(idx-1,nblk) + 1

   else

      ! non local block

      if(iflag == 0) then

         local_index = 0

      else

         local_index = (iblk/num_procs)*nblk

         if(mod(iblk,num_procs) > my_proc) local_index = local_index + nblk

         if(iflag>0) local_index = local_index + 1
      endif
   endif

end function local_index

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

subroutine cholesky_real(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  cholesky_real: Cholesky factorization of a real symmetric matrix
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be factorized.
!              Distribution is like in Scalapack.
!              Only upper triangle is needs to be set.
!              On return, the upper triangle contains the Cholesky factor
!              and the lower triangle is set to 0.
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   real*8 a(lda,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
   integer n, nc, i, info
   integer lcs, lce, lrs, lre
   integer tile_size, l_rows_tile, l_cols_tile

   real*8, allocatable:: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   ! 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


   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a

   allocate(tmp1(nblk*nblk))
   allocate(tmp2(nblk,nblk))
   tmp1 = 0
   tmp2 = 0

   allocate(tmatr(l_rows,nblk))
   allocate(tmatc(l_cols,nblk))
   tmatr = 0
   tmatc = 0


   do n = 1, na, nblk

      ! Calculate first local row and column of the still remaining matrix
      ! on the local processor

      l_row1 = local_index(n, my_prow, np_rows, nblk, +1)
      l_col1 = local_index(n, my_pcol, np_cols, nblk, +1)

      l_rowx = local_index(n+nblk, my_prow, np_rows, nblk, +1)
      l_colx = local_index(n+nblk, my_pcol, np_cols, nblk, +1)

      if(n+nblk > na) then

         ! This is the last step, just do a Cholesky-Factorization
         ! of the remaining block

         if(my_prow==prow(n) .and. my_pcol==pcol(n)) then

            call dpotrf('U',na-n+1,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in dpotrf"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

         endif

         exit ! Loop

      endif


      if(my_prow==prow(n)) then

         if(my_pcol==pcol(n)) then

            ! The process owning the upper left remaining block does the
            ! Cholesky-Factorization of this block

            call dpotrf('U',nblk,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in dpotrf"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

            nc = 0
            do i=1,nblk
               tmp1(nc+1:nc+i) = a(l_row1:l_row1+i-1,l_col1+i-1)
               nc = nc+i
            enddo
         endif

         call MPI_Bcast(tmp1,nblk*(nblk+1)/2,MPI_REAL8,pcol(n),mpi_comm_cols,mpierr)

         nc = 0
         do i=1,nblk
            tmp2(1:i,i) = tmp1(nc+1:nc+i)
            nc = nc+i
         enddo

         if(l_cols-l_colx+1>0) &
            call dtrsm('L','U','T','N',nblk,l_cols-l_colx+1,1.d0,tmp2,ubound(tmp2,1),a(l_row1,l_colx),lda)

      endif

      do i=1,nblk

         if(my_prow==prow(n)) tmatc(l_colx:l_cols,i) = a(l_row1+i-1,l_colx:l_cols)
         if(l_cols-l_colx+1>0) &
            call MPI_Bcast(tmatc(l_colx,i),l_cols-l_colx+1,MPI_REAL8,prow(n),mpi_comm_rows,mpierr)

      enddo

      call elpa_transpose_vectors  (tmatc, ubound(tmatc,1), mpi_comm_cols, &
                                    tmatr, ubound(tmatr,1), mpi_comm_rows, &
                                    n, na, nblk, nblk)

      do i=0,(na-1)/tile_size
         lcs = max(l_colx,i*l_cols_tile+1)
         lce = min(l_cols,(i+1)*l_cols_tile)
         lrs = l_rowx
         lre = min(l_rows,(i+1)*l_rows_tile)
         if(lce<lcs .or. lre<lrs) cycle
         call DGEMM('N','T',lre-lrs+1,lce-lcs+1,nblk,-1.d0, &
                    tmatr(lrs,1),ubound(tmatr,1),tmatc(lcs,1),ubound(tmatc,1), &
                    1.d0,a(lrs,lcs),lda)
      enddo

   enddo

   deallocate(tmp1, tmp2, tmatr, tmatc)

   ! Set the lower triangle to 0, it contains garbage (form the above matrix multiplications)

   do i=1,na
      if(my_pcol==pcol(i)) then
         ! column i is on local processor
         l_col1 = local_index(i  , my_pcol, np_cols, nblk, +1) ! local column number
         l_row1 = local_index(i+1, my_prow, np_rows, nblk, +1) ! first row below diagonal
         a(l_row1:l_rows,l_col1) = 0
      endif
   enddo

end subroutine cholesky_real

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

subroutine invert_trm_real(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  invert_trm_real: Inverts a upper triangular matrix
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be inverted.
!              Distribution is like in Scalapack.
!              Only upper triangle is needs to be set.
!              The lower triangle is not referenced.
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   real*8 a(lda,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
   integer n, nc, i, info, ns, nb

   real*8, allocatable:: tmp1(:), tmp2(:,:), tmat1(:,:), tmat2(:,:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a

   allocate(tmp1(nblk*nblk))
   allocate(tmp2(nblk,nblk))
   tmp1 = 0
   tmp2 = 0

   allocate(tmat1(l_rows,nblk))
   allocate(tmat2(nblk,l_cols))
   tmat1 = 0
   tmat2 = 0


   ns = ((na-1)/nblk)*nblk + 1

   do n = ns,1,-nblk

      l_row1 = local_index(n, my_prow, np_rows, nblk, +1)
      l_col1 = local_index(n, my_pcol, np_cols, nblk, +1)

      nb = nblk
      if(na-n+1 < nblk) nb = na-n+1

      l_rowx = local_index(n+nb, my_prow, np_rows, nblk, +1)
      l_colx = local_index(n+nb, my_pcol, np_cols, nblk, +1)


      if(my_prow==prow(n)) then

         if(my_pcol==pcol(n)) then

            call DTRTRI('U','N',nb,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in DTRTRI"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

            nc = 0
            do i=1,nb
               tmp1(nc+1:nc+i) = a(l_row1:l_row1+i-1,l_col1+i-1)
               nc = nc+i
            enddo
         endif

         call MPI_Bcast(tmp1,nb*(nb+1)/2,MPI_REAL8,pcol(n),mpi_comm_cols,mpierr)

         nc = 0
         do i=1,nb
            tmp2(1:i,i) = tmp1(nc+1:nc+i)
            nc = nc+i
         enddo

         if(l_cols-l_colx+1>0) &
            call DTRMM('L','U','N','N',nb,l_cols-l_colx+1,1.d0,tmp2,ubound(tmp2,1),a(l_row1,l_colx),lda)

         if(l_colx<=l_cols)   tmat2(1:nb,l_colx:l_cols) = a(l_row1:l_row1+nb-1,l_colx:l_cols)
         if(my_pcol==pcol(n)) tmat2(1:nb,l_col1:l_col1+nb-1) = tmp2(1:nb,1:nb) ! tmp2 has the lower left triangle 0

      endif

      if(l_row1>1) then
         if(my_pcol==pcol(n)) then
            tmat1(1:l_row1-1,1:nb) = a(1:l_row1-1,l_col1:l_col1+nb-1)
            a(1:l_row1-1,l_col1:l_col1+nb-1) = 0
         endif

         do i=1,nb
            call MPI_Bcast(tmat1(1,i),l_row1-1,MPI_REAL8,pcol(n),mpi_comm_cols,mpierr)
         enddo
      endif

      if(l_cols-l_col1+1>0) &
         call MPI_Bcast(tmat2(1,l_col1),(l_cols-l_col1+1)*nblk,MPI_REAL8,prow(n),mpi_comm_rows,mpierr)

      if(l_row1>1 .and. l_cols-l_col1+1>0) &
         call dgemm('N','N',l_row1-1,l_cols-l_col1+1,nb, -1.d0, &
                    tmat1,ubound(tmat1,1),tmat2(1,l_col1),ubound(tmat2,1), &
                    1.d0, a(1,l_col1),lda)

   enddo

   deallocate(tmp1, tmp2, tmat1, tmat2)

end subroutine invert_trm_real

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

subroutine cholesky_complex(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  cholesky_complex: Cholesky factorization of a complex hermitian matrix
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be factorized.
!              Distribution is like in Scalapack.
!              Only upper triangle is needs to be set.
!              On return, the upper triangle contains the Cholesky factor
!              and the lower triangle is set to 0.
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   complex*16 a(lda,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
   integer n, nc, i, info
   integer lcs, lce, lrs, lre
   integer tile_size, l_rows_tile, l_cols_tile

   complex*16, allocatable:: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   ! 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


   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a

   allocate(tmp1(nblk*nblk))
   allocate(tmp2(nblk,nblk))
   tmp1 = 0
   tmp2 = 0

   allocate(tmatr(l_rows,nblk))
   allocate(tmatc(l_cols,nblk))
   tmatr = 0
   tmatc = 0


   do n = 1, na, nblk

      ! Calculate first local row and column of the still remaining matrix
      ! on the local processor

      l_row1 = local_index(n, my_prow, np_rows, nblk, +1)
      l_col1 = local_index(n, my_pcol, np_cols, nblk, +1)

      l_rowx = local_index(n+nblk, my_prow, np_rows, nblk, +1)
      l_colx = local_index(n+nblk, my_pcol, np_cols, nblk, +1)

      if(n+nblk > na) then

         ! This is the last step, just do a Cholesky-Factorization
         ! of the remaining block

         if(my_prow==prow(n) .and. my_pcol==pcol(n)) then

            call zpotrf('U',na-n+1,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in zpotrf"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

         endif

         exit ! Loop

      endif


      if(my_prow==prow(n)) then

         if(my_pcol==pcol(n)) then

            ! The process owning the upper left remaining block does the
            ! Cholesky-Factorization of this block

            call zpotrf('U',nblk,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in zpotrf"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

            nc = 0
            do i=1,nblk
               tmp1(nc+1:nc+i) = a(l_row1:l_row1+i-1,l_col1+i-1)
               nc = nc+i
            enddo
         endif

         call MPI_Bcast(tmp1,nblk*(nblk+1)/2,MPI_DOUBLE_COMPLEX,pcol(n),mpi_comm_cols,mpierr)

         nc = 0
         do i=1,nblk
            tmp2(1:i,i) = tmp1(nc+1:nc+i)
            nc = nc+i
         enddo

         if(l_cols-l_colx+1>0) &
            call ztrsm('L','U','C','N',nblk,l_cols-l_colx+1,(1.d0,0.d0),tmp2,ubound(tmp2,1),a(l_row1,l_colx),lda)

      endif

      do i=1,nblk

         if(my_prow==prow(n)) tmatc(l_colx:l_cols,i) = conjg(a(l_row1+i-1,l_colx:l_cols))
         if(l_cols-l_colx+1>0) &
            call MPI_Bcast(tmatc(l_colx,i),l_cols-l_colx+1,MPI_DOUBLE_COMPLEX,prow(n),mpi_comm_rows,mpierr)

      enddo

      call elpa_transpose_vectors  (tmatc, 2*ubound(tmatc,1), mpi_comm_cols, &
                                    tmatr, 2*ubound(tmatr,1), mpi_comm_rows, &
                                    2*n-1, 2*na, nblk, 2*nblk)


      do i=0,(na-1)/tile_size
         lcs = max(l_colx,i*l_cols_tile+1)
         lce = min(l_cols,(i+1)*l_cols_tile)
         lrs = l_rowx
         lre = min(l_rows,(i+1)*l_rows_tile)
         if(lce<lcs .or. lre<lrs) cycle
         call ZGEMM('N','C',lre-lrs+1,lce-lcs+1,nblk,(-1.d0,0.d0), &
                    tmatr(lrs,1),ubound(tmatr,1),tmatc(lcs,1),ubound(tmatc,1), &
                    (1.d0,0.d0),a(lrs,lcs),lda)
      enddo

   enddo

   deallocate(tmp1, tmp2, tmatr, tmatc)

   ! Set the lower triangle to 0, it contains garbage (form the above matrix multiplications)

   do i=1,na
      if(my_pcol==pcol(i)) then
         ! column i is on local processor
         l_col1 = local_index(i  , my_pcol, np_cols, nblk, +1) ! local column number
         l_row1 = local_index(i+1, my_prow, np_rows, nblk, +1) ! first row below diagonal
         a(l_row1:l_rows,l_col1) = 0
      endif
   enddo

end subroutine cholesky_complex

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

subroutine invert_trm_complex(na, a, lda, nblk, mpi_comm_rows, mpi_comm_cols)

!-------------------------------------------------------------------------------
!  invert_trm_complex: Inverts a upper triangular matrix
!
!  Parameters
!
!  na          Order of matrix
!
!  a(lda,*)    Distributed matrix which should be inverted.
!              Distribution is like in Scalapack.
!              Only upper triangle is needs to be set.
!              The lower triangle is not referenced.
!
!  lda         Leading dimension of a
!
!  nblk        blocksize of cyclic distribution, must be the same in both directions!
!
!  mpi_comm_rows
!  mpi_comm_cols
!              MPI-Communicators for rows/columns
!
!-------------------------------------------------------------------------------

   implicit none

   integer na, lda, nblk, mpi_comm_rows, mpi_comm_cols
   complex*16 a(lda,*)

   integer my_prow, my_pcol, np_rows, np_cols, mpierr
   integer l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
   integer n, nc, i, info, ns, nb

   complex*16, allocatable:: tmp1(:), tmp2(:,:), tmat1(:,:), tmat2(:,:)

   integer pcol, prow
   pcol(i) = MOD((i-1)/nblk,np_cols) !Processor col for global col number
   prow(i) = MOD((i-1)/nblk,np_rows) !Processor row for global row number


   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)

   l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
   l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a

   allocate(tmp1(nblk*nblk))
   allocate(tmp2(nblk,nblk))
   tmp1 = 0
   tmp2 = 0

   allocate(tmat1(l_rows,nblk))
   allocate(tmat2(nblk,l_cols))
   tmat1 = 0
   tmat2 = 0


   ns = ((na-1)/nblk)*nblk + 1

   do n = ns,1,-nblk

      l_row1 = local_index(n, my_prow, np_rows, nblk, +1)
      l_col1 = local_index(n, my_pcol, np_cols, nblk, +1)

      nb = nblk
      if(na-n+1 < nblk) nb = na-n+1

      l_rowx = local_index(n+nb, my_prow, np_rows, nblk, +1)
      l_colx = local_index(n+nb, my_pcol, np_cols, nblk, +1)


      if(my_prow==prow(n)) then

         if(my_pcol==pcol(n)) then

            call ZTRTRI('U','N',nb,a(l_row1,l_col1),lda,info)
            if(info/=0) then
               print *,"Error in ZTRTRI"
               call MPI_Abort(MPI_COMM_WORLD,1,mpierr)
            endif

            nc = 0
            do i=1,nb
               tmp1(nc+1:nc+i) = a(l_row1:l_row1+i-1,l_col1+i-1)
               nc = nc+i
            enddo
         endif

         call MPI_Bcast(tmp1,nb*(nb+1)/2,MPI_DOUBLE_COMPLEX,pcol(n),mpi_comm_cols,mpierr)

         nc = 0
         do i=1,nb
            tmp2(1:i,i) = tmp1(nc+1:nc+i)
            nc = nc+i
         enddo

         if(l_cols-l_colx+1>0) &
            call ZTRMM('L','U','N','N',nb,l_cols-l_colx+1,(1.d0,0.d0),tmp2,ubound(tmp2,1),a(l_row1,l_colx),lda)

         if(l_colx<=l_cols)   tmat2(1:nb,l_colx:l_cols) = a(l_row1:l_row1+nb-1,l_colx:l_cols)
         if(my_pcol==pcol(n)) tmat2(1:nb,l_col1:l_col1+nb-1) = tmp2(1:nb,1:nb) ! tmp2 has the lower left triangle 0

      endif

      if(l_row1>1) then
         if(my_pcol==pcol(n)) then
            tmat1(1:l_row1-1,1:nb) = a(1:l_row1-1,l_col1:l_col1+nb-1)
            a(1:l_row1-1,l_col1:l_col1+nb-1) = 0
         endif

         do i=1,nb
            call MPI_Bcast(tmat1(1,i),l_row1-1,MPI_DOUBLE_COMPLEX,pcol(n),mpi_comm_cols,mpierr)
         enddo
      endif

      if(l_cols-l_col1+1>0) &
         call MPI_Bcast(tmat2(1,l_col1),(l_cols-l_col1+1)*nblk,MPI_DOUBLE_COMPLEX,prow(n),mpi_comm_rows,mpierr)

      if(l_row1>1 .and. l_cols-l_col1+1>0) &
         call ZGEMM('N','N',l_row1-1,l_cols-l_col1+1,nb, (-1.d0,0.d0), &
                    tmat1,ubound(tmat1,1),tmat2(1,l_col1),ubound(tmat2,1), &
                    (1.d0,0.d0), a(1,l_col1),lda)

   enddo

   deallocate(tmp1, tmp2, tmat1, tmat2)

end subroutine invert_trm_complex

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

integer function least_common_multiple(a, b)

   ! Returns the least common multiple of a and b
   ! There may be more efficient ways to do this, we use the most simple approach

   implicit none
   integer, intent(in) :: a, b

   do least_common_multiple = a, a*(b-1), a
      if(mod(least_common_multiple,b)==0) exit
   enddo
   ! if the loop is left regularly, least_common_multiple = a*b

end function

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

subroutine hh_transform_real(alpha, xnorm_sq, xf, tau)

   ! Similar to LAPACK routine DLARFP, but uses ||x||**2 instead of x(:)
   ! and returns the factor xf by which x has to be scaled.
   ! It also hasn't the special handling for numbers < 1.d-300 or > 1.d150
   ! since this would be expensive for the parallel implementation.

   implicit none
   real*8, intent(inout) :: alpha
   real*8, intent(in)    :: xnorm_sq
   real*8, intent(out)   :: xf, tau

   real*8 BETA

   if( XNORM_SQ==0. ) then

      if( ALPHA>=0. ) then
         TAU = 0.
      else
         TAU = 2.
         ALPHA = -ALPHA
      endif
      XF = 0.

   else

      BETA = SIGN( SQRT( ALPHA**2 + XNORM_SQ ), ALPHA )
      ALPHA = ALPHA + BETA
      IF( BETA<0 ) THEN
         BETA = -BETA
         TAU = -ALPHA / BETA
      ELSE
         ALPHA = XNORM_SQ / ALPHA
         TAU = ALPHA / BETA
         ALPHA = -ALPHA
      END IF
      XF = 1./ALPHA
      ALPHA = BETA

   endif

end subroutine


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

subroutine hh_transform_complex(alpha, xnorm_sq, xf, tau)

   ! Similar to LAPACK routine ZLARFP, but uses ||x||**2 instead of x(:)
   ! and returns the factor xf by which x has to be scaled.
   ! It also hasn't the special handling for numbers < 1.d-300 or > 1.d150
   ! since this would be expensive for the parallel implementation.

   implicit none
   complex*16, intent(inout) :: alpha
   real*8, intent(in)        :: xnorm_sq
   complex*16, intent(out)   :: xf, tau

   real*8 ALPHR, ALPHI, BETA

   ALPHR = DBLE( ALPHA )
   ALPHI = DIMAG( ALPHA )

   if( XNORM_SQ==0. .AND. ALPHI==0. ) then

      if( ALPHR>=0. ) then
         TAU = 0.
      else
         TAU = 2.
         ALPHA = -ALPHA
      endif
      XF = 0.

   else

      BETA = SIGN( SQRT( ALPHR**2 + ALPHI**2 + XNORM_SQ ), ALPHR )
      ALPHA = ALPHA + BETA
      IF( BETA<0 ) THEN
         BETA = -BETA
         TAU = -ALPHA / BETA
      ELSE
         ALPHR = ALPHI * (ALPHI/DBLE( ALPHA ))
         ALPHR = ALPHR + XNORM_SQ/DBLE( ALPHA )
         TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA )
         ALPHA = DCMPLX( -ALPHR, ALPHI )
      END IF
      XF = 1./ALPHA
      ALPHA = BETA

   endif

end subroutine

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

end module ELPA1

! --------------------------------------------------------------------------------------------------
! Please note that the following routines are outside of the module ELPA1
! so that they can be used with real or complex data
! --------------------------------------------------------------------------------------------------

subroutine elpa_transpose_vectors(vmat_s,ld_s,comm_s,vmat_t,ld_t,comm_t,nvs,nvr,nvc,nblk)

!-------------------------------------------------------------------------------
! This routine transposes an array of vectors which are distributed in
! communicator comm_s into its transposed form distributed in communicator comm_t.
! There must be an identical copy of vmat_s in every communicator comm_s.
! After this routine, there is an identical copy of vmat_t in every communicator comm_t.
!
! vmat_s    original array of vectors
! ld_s      leading dimension of vmat_s
! comm_s    communicator over which vmat_s is distributed
! vmat_t    array of vectors in transposed form
! ld_t      leading dimension of vmat_t
! comm_t    communicator over which vmat_t is distributed
! nvs       global index where to start in vmat_s/vmat_t
!           Please note: this is kind of a hint, some values before nvs will be
!           accessed in vmat_s/put into vmat_t
! nvr       global length of vmat_s/vmat_t
! nvc       number of columns in vmat_s/vmat_t
! nblk      block size of block cyclic distribution
!
!-------------------------------------------------------------------------------

   use ELPA1 ! for least_common_multiple

   implicit none

   include 'mpif.h'

   integer, intent(in)   :: ld_s, comm_s, ld_t, comm_t, nvs, nvr, nvc, nblk
   real*8, intent(in)    :: vmat_s(ld_s,nvc)
   real*8, intent(inout) :: vmat_t(ld_t,nvc)

   real*8, allocatable :: aux(:)
   integer myps, mypt, nps, npt
   integer n, lc, k, i, ips, ipt, ns, nl, mpierr
   integer lcm_s_t, nblks_tot, nblks_comm, nblks_skip

   call mpi_comm_rank(comm_s,myps,mpierr)
   call mpi_comm_size(comm_s,nps ,mpierr)
   call mpi_comm_rank(comm_t,mypt,mpierr)
   call mpi_comm_size(comm_t,npt ,mpierr)

   ! The basic idea of this routine is that for every block (in the block cyclic
   ! distribution), the processor within comm_t which owns the diagonal
   ! broadcasts its values of vmat_s to all processors within comm_t.
   ! Of course this has not to be done for every block separately, since
   ! the communictation pattern repeats in the global matrix after
   ! the least common multiple of (nps,npt) blocks

   lcm_s_t   = least_common_multiple(nps,npt) ! least common multiple of nps, npt

   nblks_tot = (nvr+nblk-1)/nblk ! number of blocks corresponding to nvr

   ! Get the number of blocks to be skipped at the begin.
   ! This must be a multiple of lcm_s_t (else it is getting complicated),
   ! thus some elements before nvs will be accessed/set.

   nblks_skip = ((nvs-1)/(nblk*lcm_s_t))*lcm_s_t

   allocate(aux( ((nblks_tot-nblks_skip+lcm_s_t-1)/lcm_s_t) * nblk * nvc ))

   do n = 0, lcm_s_t-1

      ips = mod(n,nps)
      ipt = mod(n,npt)

      if(mypt == ipt) then

         nblks_comm = (nblks_tot-nblks_skip-n+lcm_s_t-1)/lcm_s_t
         if(nblks_comm==0) cycle

         if(myps == ips) then
            k = 0
            do lc=1,nvc
               do i = nblks_skip+n, nblks_tot-1, lcm_s_t
                  ns = (i/nps)*nblk ! local start of block i
                  nl = min(nvr-i*nblk,nblk) ! length
                  aux(k+1:k+nl) = vmat_s(ns+1:ns+nl,lc)
                  k = k+nblk
               enddo
            enddo
         endif

         call MPI_Bcast(aux,nblks_comm*nblk*nvc,MPI_REAL8,ips,comm_s,mpierr)

         k = 0
         do lc=1,nvc
            do i = nblks_skip+n, nblks_tot-1, lcm_s_t
               ns = (i/npt)*nblk ! local start of block i
               nl = min(nvr-i*nblk,nblk) ! length
               vmat_t(ns+1:ns+nl,lc) = aux(k+1:k+nl)
               k = k+nblk
            enddo
         enddo

      endif

   enddo

   deallocate(aux)

end subroutine

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

subroutine elpa_reduce_add_vectors(vmat_s,ld_s,comm_s,vmat_t,ld_t,comm_t,nvr,nvc,nblk)

!-------------------------------------------------------------------------------
! This routine does a reduce of all vectors in vmat_s over the communicator comm_t.
! The result of the reduce is gathered on the processors owning the diagonal
! and added to the array of vectors vmat_t (which is distributed over comm_t).
!
! Opposed to elpa_transpose_vectors, there is NO identical copy of vmat_s
! in the different members within vmat_t (else a reduce wouldn't be necessary).
! After this routine, an allreduce of vmat_t has to be done.
!
! vmat_s    array of vectors to be reduced and added
! ld_s      leading dimension of vmat_s
! comm_s    communicator over which vmat_s is distributed
! vmat_t    array of vectors to which vmat_s is added
! ld_t      leading dimension of vmat_t
! comm_t    communicator over which vmat_t is distributed
! nvr       global length of vmat_s/vmat_t
! nvc       number of columns in vmat_s/vmat_t
! nblk      block size of block cyclic distribution
!
!-------------------------------------------------------------------------------

   use ELPA1 ! for least_common_multiple

   implicit none

   include 'mpif.h'

   integer, intent(in)   :: ld_s, comm_s, ld_t, comm_t, nvr, nvc, nblk
   real*8, intent(in)    :: vmat_s(ld_s,nvc)
   real*8, intent(inout) :: vmat_t(ld_t,nvc)

   real*8, allocatable :: aux1(:), aux2(:)
   integer myps, mypt, nps, npt
   integer n, lc, k, i, ips, ipt, ns, nl, mpierr
   integer lcm_s_t, nblks_tot

   call mpi_comm_rank(comm_s,myps,mpierr)
   call mpi_comm_size(comm_s,nps ,mpierr)
   call mpi_comm_rank(comm_t,mypt,mpierr)
   call mpi_comm_size(comm_t,npt ,mpierr)

   ! Look to elpa_transpose_vectors for the basic idea!

   ! The communictation pattern repeats in the global matrix after
   ! the least common multiple of (nps,npt) blocks

   lcm_s_t   = least_common_multiple(nps,npt) ! least common multiple of nps, npt

   nblks_tot = (nvr+nblk-1)/nblk ! number of blocks corresponding to nvr

   allocate(aux1( ((nblks_tot+lcm_s_t-1)/lcm_s_t) * nblk * nvc ))
   allocate(aux2( ((nblks_tot+lcm_s_t-1)/lcm_s_t) * nblk * nvc ))
   aux1(:) = 0
   aux2(:) = 0

   do n = 0, lcm_s_t-1

      ips = mod(n,nps)
      ipt = mod(n,npt)

      if(myps == ips) then

         k = 0
         do lc=1,nvc
            do i = n, nblks_tot-1, lcm_s_t
               ns = (i/nps)*nblk ! local start of block i
               nl = min(nvr-i*nblk,nblk) ! length
               aux1(k+1:k+nl) = vmat_s(ns+1:ns+nl,lc)
               k = k+nblk
            enddo
         enddo

         if(k>0) call mpi_reduce(aux1,aux2,k,MPI_REAL8,MPI_SUM,ipt,comm_t,mpierr)

         if(mypt == ipt) then
            k = 0
            do lc=1,nvc
               do i = n, nblks_tot-1, lcm_s_t
                  ns = (i/npt)*nblk ! local start of block i
                  nl = min(nvr-i*nblk,nblk) ! length
                  vmat_t(ns+1:ns+nl,lc) = vmat_t(ns+1:ns+nl,lc) + aux2(k+1:k+nl)
                  k = k+nblk
               enddo
            enddo
         endif

      endif

   enddo

   deallocate(aux1)
   deallocate(aux2)

end subroutine

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