elpa2_trans_ev_tridi_to_band_template.F90 88.7 KB
Newer Older
Andreas Marek's avatar
Andreas Marek committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), formerly known as
!      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!    This particular source code file contains additions, changes and
!    enhancements authored by Intel Corporation which is not part of
!    the ELPA consortium.
!
!    More information can be found here:
!    http://elpa.mpcdf.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
!
! ELPA1 -- Faster replacements for ScaLAPACK symmetric eigenvalue routines
!
! Copyright of the original code rests with the authors inside the ELPA
! consortium. The copyright of any additional modifications shall rest
! with their original authors, but shall adhere to the licensing terms
! distributed along with the original code in the file "COPYING".

53
#include "../general/sanity.F90"
Andreas Marek's avatar
Andreas Marek committed
54 55

  subroutine trans_ev_tridi_to_band_&
56 57 58
    &MATH_DATATYPE&
    &_&
    &PRECISION &
59
    (obj, na, nev, nblk, nbw, q, q_dev, ldq, matrixCols,         &
Andreas Marek's avatar
Andreas Marek committed
60
     hh_trans, mpi_comm_rows, mpi_comm_cols, wantDebug, useGPU, max_threads, success, &
61
     kernel)
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90

    !-------------------------------------------------------------------------------
    !  trans_ev_tridi_to_band_real/complex:
    !  Transforms the eigenvectors of a tridiagonal matrix back to the eigenvectors of the band matrix
    !
    !  Parameters
    !
    !  na          Order of matrix a, number of rows of matrix q
    !
    !  nev         Number eigenvectors to compute (= columns of matrix q)
    !
    !  nblk        blocksize of cyclic distribution, must be the same in both directions!
    !
    !  nb          semi bandwith
    !
    !  q           On input: Eigenvectors of tridiagonal matrix
    !              On output: Transformed eigenvectors
    !              Distribution is like in Scalapack.
    !
    !  q_dev       GPU device pointer to q
    !
    !  ldq         Leading dimension of q
    !  matrixCols  local columns of matrix q
    !
    !  mpi_comm_rows
    !  mpi_comm_cols
    !              MPI-Communicators for rows/columns/both
    !
    !-------------------------------------------------------------------------------
91
      use elpa_abstract_impl
92
      use elpa2_workload
Andreas Marek's avatar
Andreas Marek committed
93
      use pack_unpack_cpu
Andreas Marek's avatar
Andreas Marek committed
94
      use pack_unpack_gpu
95
      use compute_hh_trafo
96 97 98
      use cuda_functions
      use precision
      use iso_c_binding
Andreas Marek's avatar
Andreas Marek committed
99
#ifdef WITH_OPENMP
100
      ! use omp_lib
Andreas Marek's avatar
Andreas Marek committed
101
#endif
102
      implicit none
103
#include "../general/precision_kinds.F90"
104
      class(elpa_abstract_impl_t), intent(inout) :: obj
Andreas Marek's avatar
Andreas Marek committed
105
      logical, intent(in)                        :: useGPU
106

Andreas Marek's avatar
Andreas Marek committed
107 108
      integer(kind=ik), intent(in)               :: kernel
      integer(kind=ik), intent(in)               :: na, nev, nblk, nbw, ldq, matrixCols, mpi_comm_rows, mpi_comm_cols
109 110

#ifdef USE_ASSUMED_SIZE
111
      MATH_DATATYPE(kind=rck)                    :: q(ldq,*)
112
#else
113
      MATH_DATATYPE(kind=rck)                    :: q(ldq,matrixCols)
114 115
#endif

116
      MATH_DATATYPE(kind=rck), intent(in)        :: hh_trans(:,:)
Andreas Marek's avatar
Andreas Marek committed
117
      integer(kind=c_intptr_t)                   :: q_dev
Andreas Marek's avatar
Andreas Marek committed
118

Andreas Marek's avatar
Andreas Marek committed
119 120 121 122 123 124
      integer(kind=ik)                           :: np_rows, my_prow, np_cols, my_pcol
      integer(kind=ik)                           :: i, j, ip, sweep, nbuf, l_nev, a_dim2
      integer(kind=ik)                           :: current_n, current_local_n, current_n_start, current_n_end
      integer(kind=ik)                           :: next_n, next_local_n, next_n_start, next_n_end
      integer(kind=ik)                           :: bottom_msg_length, top_msg_length, next_top_msg_length
      integer(kind=ik)                           :: stripe_width, last_stripe_width, stripe_count
125
#ifdef WITH_OPENMP
Andreas Marek's avatar
Andreas Marek committed
126
      integer(kind=ik)                           :: thread_width, csw, b_off, b_len
127
#endif
Andreas Marek's avatar
Andreas Marek committed
128 129 130
      integer(kind=ik)                           :: num_result_blocks, num_result_buffers, num_bufs_recvd
      integer(kind=ik)                           :: a_off, current_tv_off, max_blk_size
      integer(kind=ik)                           :: mpierr, src, src_offset, dst, offset, nfact, num_blk
131

Andreas Marek's avatar
Andreas Marek committed
132
      logical                                    :: flag
133
#ifdef WITH_OPENMP
134
      MATH_DATATYPE(kind=rck), pointer           :: aIntern(:,:,:,:)
135
#else
136
      MATH_DATATYPE(kind=rck), pointer           :: aIntern(:,:,:)
137
#endif
138
      MATH_DATATYPE(kind=rck)                    :: a_var
139

Andreas Marek's avatar
Andreas Marek committed
140
      type(c_ptr)                                :: aIntern_ptr
141

142 143
      MATH_DATATYPE(kind=rck)   , allocatable    :: row(:)
      MATH_DATATYPE(kind=rck)   , allocatable    :: row_group(:,:)
144 145

#ifdef WITH_OPENMP
146 147
      MATH_DATATYPE(kind=rck), allocatable       :: top_border_send_buffer(:,:), top_border_recv_buffer(:,:)
      MATH_DATATYPE(kind=rck), allocatable       :: bottom_border_send_buffer(:,:), bottom_border_recv_buffer(:,:)
148
#else
149 150
      MATH_DATATYPE(kind=rck), allocatable       :: top_border_send_buffer(:,:,:), top_border_recv_buffer(:,:,:)
      MATH_DATATYPE(kind=rck), allocatable       :: bottom_border_send_buffer(:,:,:), bottom_border_recv_buffer(:,:,:)
151 152
#endif

153 154
      integer(kind=c_intptr_t)                   :: aIntern_dev
      integer(kind=c_intptr_t)                   :: bcast_buffer_dev
155
      integer(kind=c_intptr_t)                   :: num
Andreas Marek's avatar
Andreas Marek committed
156
      integer(kind=c_intptr_t)                   :: dev_offset, dev_offset_1
157 158 159 160 161
      integer(kind=c_intptr_t)                   :: row_dev
      integer(kind=c_intptr_t)                   :: row_group_dev
      integer(kind=c_intptr_t)                   :: hh_tau_dev
      integer(kind=c_intptr_t)                   :: hh_dot_dev
      integer(kind=ik)                           :: row_group_size, unpack_idx
162

163 164
      integer(kind=ik)                           :: n_times
      integer(kind=ik)                           :: top, chunk, this_chunk
165

166 167
      MATH_DATATYPE(kind=rck), allocatable       :: result_buffer(:,:,:)
      MATH_DATATYPE(kind=rck), allocatable       :: bcast_buffer(:,:)
168

169
      integer(kind=ik)                           :: n_off
170

171 172 173
      integer(kind=ik), allocatable              :: result_send_request(:), result_recv_request(:), limits(:)
      integer(kind=ik), allocatable              :: top_send_request(:), bottom_send_request(:)
      integer(kind=ik), allocatable              :: top_recv_request(:), bottom_recv_request(:)
174 175 176

      ! MPI send/recv tags, arbitrary

177 178 179
      integer(kind=ik), parameter                :: bottom_recv_tag = 111
      integer(kind=ik), parameter                :: top_recv_tag    = 222
      integer(kind=ik), parameter                :: result_recv_tag = 333
Andreas Marek's avatar
Andreas Marek committed
180 181 182

      integer(kind=ik), intent(in)               :: max_threads
 
183
#ifdef WITH_OPENMP
Andreas Marek's avatar
Andreas Marek committed
184
      integer(kind=ik)                           :: my_thread
185 186 187 188
#endif


      ! Just for measuring the kernel performance
189
      real(kind=c_double)                        :: kernel_time, kernel_time_recv ! MPI_WTIME always needs double
190
      ! long integer
191
      integer(kind=lik)                          :: kernel_flops, kernel_flops_recv
192

193 194 195 196 197
      logical, intent(in)                        :: wantDebug
      logical                                    :: success
      integer(kind=ik)                           :: istat, print_flops
      character(200)                             :: errorMessage
      logical                                    :: successCUDA
198
#ifndef WITH_MPI
199
      integer(kind=ik)                           :: j1
200
#endif
201 202 203 204
      integer(kind=c_intptr_t), parameter        :: size_of_datatype = size_of_&
                                                                     &PRECISION&
                                                                     &_&
                                                                     &MATH_DATATYPE
205

206
      call obj%timer%start("trans_ev_tridi_to_band_&
207
      &MATH_DATATYPE&
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
208
      &" // &
209 210 211
      &PRECISION_SUFFIX &
      )

212
      n_times = 0
213 214 215 216 217 218 219 220 221
      if (useGPU) then
        unpack_idx = 0
        row_group_size = 0
      endif

      success = .true.
      kernel_time = 0.0
      kernel_flops = 0

222
      if (wantDebug) call obj%timer%start("mpi_communication")
223 224 225 226
      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)
227
      if (wantDebug) call obj%timer%stop("mpi_communication")
228 229 230 231 232

      if (mod(nbw,nblk)/=0) then
        if (my_prow==0 .and. my_pcol==0) then
          if (wantDebug) then
            write(error_unit,*) 'ELPA2_trans_ev_tridi_to_band_&
233 234
                                &MATH_DATATYPE&
                                &: ERROR: nbw=',nbw,', nblk=',nblk
235
            write(error_unit,*) 'ELPA2_trans_ev_tridi_to_band_&
236 237
                                &MATH_DATATYPE&
                                &: band backtransform works only for nbw==n*nblk'
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
          endif
          success = .false.
          return
        endif
      endif

      nfact = nbw / nblk


      ! local number of eigenvectors
      l_nev = local_index(nev, my_pcol, np_cols, nblk, -1)

      if (l_nev==0) then
#ifdef WITH_OPENMP
        thread_width = 0
#endif
        stripe_width = 0
        stripe_count = 0
        last_stripe_width = 0

      else ! l_nev

#if WITH_OPENMP
        ! Suggested stripe width is 48 since 48*64 real*8 numbers should fit into
        ! every primary cache
263
        ! Suggested stripe width is 48 - should this be reduced for the complex case ???
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294

        if (useGPU) then
          stripe_width = 256 ! Must be a multiple of 4
          stripe_count = (l_nev - 1) / stripe_width + 1
        else ! useGPU
          ! openmp only in non-GPU case
          thread_width = (l_nev-1)/max_threads + 1 ! number of eigenvectors per OMP thread
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
          stripe_width = 48 ! Must be a multiple of 4
#else
          stripe_width = 96 ! Must be a multiple of 8
#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
          stripe_width = 48 ! Must be a multiple of 2
#else
          stripe_width = 48 ! Must be a multiple of 4
#endif
#endif /* COMPLEXCASE */

          stripe_count = (thread_width-1)/stripe_width + 1

          ! Adapt stripe width so that last one doesn't get too small

          stripe_width = (thread_width-1)/stripe_count + 1

#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
295
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
296
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
297
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
298 299 300 301

            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 8 because of AVX-512 memory alignment of 64 bytes
                                                  ! (8 * sizeof(double) == 64)

Andreas Marek's avatar
Retab  
Andreas Marek committed
302
    else
303
            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 4 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
304
                                            ! (4 * sizeof(double) == 32)
305
          endif
306
#else
307
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
308
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
309
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
310 311 312


            stripe_width = ((stripe_width+15)/16)*16 ! Must be a multiple of 16 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
313
                                               ! (16 * sizeof(float) == 64)
314

Andreas Marek's avatar
Retab  
Andreas Marek committed
315
    else
316
            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 8 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
317
                                            ! (8 * sizeof(float) == 32)
318
          endif
319 320 321 322 323
#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
324
          if (kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK1 .or. &
325
              kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK2) then
326

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
327
            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 4 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
328
                                            ! (4 * sizeof(double complex) == 64)
329

Andreas Marek's avatar
Retab  
Andreas Marek committed
330
    else
331 332

            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 2 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
333 334
                                            ! (2 * sizeof(double complex) == 32)
    endif
335
#else
336

337
          if (kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK1 .or. &
338
              kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK2) then
339 340

            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 8 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
341
                                            ! (8 * sizeof(float complex) == 64)
342 343 344

          else
            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 4 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
345 346
                                            ! (4 * sizeof(float complex) == 32)
    endif
347 348 349
#endif
#endif /* COMPLEXCASE */

350
#if REALCASE == 1
351
        last_stripe_width = l_nev - (stripe_count-1)*stripe_width
352 353 354 355 356
#endif
#if COMPLEXCASE == 1
! only needed in no OMP case check thsis
! last_stripe_width = l_nev - (stripe_count-1)*stripe_width
#endif
357 358

        endif ! useGPU
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
359

360
#else /* WITH_OPENMP */
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
361

362 363
        ! Suggested stripe width is 48 since 48*64 real*8 numbers should fit into
        ! every primary cache
364
        ! Suggested stripe width is 48 - should this be reduced for the complex case ???
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394

        if (useGPU) then
          stripe_width = 256 ! Must be a multiple of 4
          stripe_count = (l_nev - 1) / stripe_width + 1

        else ! useGPU
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
          stripe_width = 48 ! Must be a multiple of 4
#else
          stripe_width = 96 ! Must be a multiple of 8
#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
          stripe_width = 48 ! Must be a multiple of 2
#else
          stripe_width = 48 ! Must be a multiple of 4
#endif
#endif /* COMPLEXCASE */

          stripe_count = (l_nev-1)/stripe_width + 1

          ! Adapt stripe width so that last one doesn't get too small

          stripe_width = (l_nev-1)/stripe_count + 1

#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
395
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
396
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
397
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
398 399 400 401

            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 8 because of AVX-512 memory alignment of 64 bytes
                                                  ! (8 * sizeof(double) == 64)

Andreas Marek's avatar
Retab  
Andreas Marek committed
402
    else
403
            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 4 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
404
                                            ! (4 * sizeof(double) == 32)
405
          endif
406
#else
407
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
408
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
409
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
410 411 412


            stripe_width = ((stripe_width+15)/16)*16 ! Must be a multiple of 16 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
413
                                               ! (16 * sizeof(float) == 64)
414

Andreas Marek's avatar
Retab  
Andreas Marek committed
415
    else
416
            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 8 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
417
                                            ! (8 * sizeof(float) == 32)
418
          endif
419 420 421 422 423
#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
424

425
          if (kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK1 .or. &
426
              kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK2) then
427

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
428
            stripe_width = ((stripe_width+7)/8)*8 ! Must be a multiple of 4 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
429
                                            ! (4 * sizeof(double complex) == 64)
430

Andreas Marek's avatar
Retab  
Andreas Marek committed
431
    else
432 433

            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 2 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
434 435
                                            ! (2 * sizeof(double complex) == 32)
    endif
436
#else
437

438
          if (kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK1 .or. &
439
              kernel .eq. ELPA_2STAGE_COMPLEX_AVX512_BLOCK2) then
440

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
441
            stripe_width = ((stripe_width+15)/16)*16 ! Must be a multiple of 8 because of AVX-512 memory alignment of 64 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
442
                                            ! (8 * sizeof(float complex) == 64)
443 444 445

          else
            stripe_width = ((stripe_width+3)/4)*4 ! Must be a multiple of 4 because of AVX/SSE memory alignment of 32 bytes
Andreas Marek's avatar
Retab  
Andreas Marek committed
446 447
                                            ! (4 * sizeof(float complex) == 32)
    endif
448 449 450 451 452 453 454 455 456 457 458 459 460 461
#endif
#endif /* COMPLEXCASE */
        endif ! useGPU

        last_stripe_width = l_nev - (stripe_count-1)*stripe_width

#endif /* WITH_OPENMP */
      endif ! l_nev

      ! Determine the matrix distribution at the beginning

      allocate(limits(0:np_rows), stat=istat, errmsg=errorMessage)
      if (istat .ne. 0) then
        print *,"trans_ev_tridi_to_band_&
462 463
                &MATH_DATATYPE&
                &: error when allocating limits"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
464
        stop 1
465
      endif
466
      call determine_workload(obj,na, nbw, np_rows, limits)
467 468 469 470 471 472

      max_blk_size = maxval(limits(1:np_rows) - limits(0:np_rows-1))

      a_dim2 = max_blk_size + nbw

      if (useGPU) then
473 474
        num =  (stripe_width*a_dim2*stripe_count)* size_of_datatype
        successCUDA = cuda_malloc(aIntern_dev, stripe_width*a_dim2*stripe_count* size_of_datatype)
475 476 477 478
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMalloc"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
479
          stop 1
480 481 482 483 484 485 486
        endif

        successCUDA = cuda_memset(aIntern_dev , 0, num)
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMemset"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
487
          stop 1
488 489
        endif

490
        num =  (l_nev)* size_of_datatype
491 492 493 494 495
        successCUDA = cuda_malloc( row_dev,num)
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMalloc "
Andreas Marek's avatar
Andreas Marek committed
496
          stop 1
497 498 499 500 501 502 503
        endif

        successCUDA = cuda_memset(row_dev , 0, num)
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMemset "
Andreas Marek's avatar
Andreas Marek committed
504
          stop 1
505 506 507 508 509 510 511 512
        endif

        ! "row_group" and "row_group_dev" are needed for GPU optimizations
        allocate(row_group(l_nev, nblk), stat=istat, errmsg=errorMessage)
        if (istat .ne. 0) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error when allocating row_group"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
513
          stop 1
514 515
        endif

516
        row_group(:, :) = 0.0_rck
517
        num =  (l_nev*nblk)* size_of_datatype
518 519 520 521 522
        successCUDA = cuda_malloc(row_group_dev, num)
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMalloc"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
523
          stop 1
524 525 526 527 528 529 530
        endif

        successCUDA = cuda_memset(row_group_dev , 0, num)
        if (.not.(successCUDA)) then
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
          &: error in cudaMemset"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
531
          stop 1
532 533 534 535 536 537 538 539 540 541
        endif

      else ! GPUs are not used

#if 0
! realcase or complexcase
!DEC$ ATTRIBUTES ALIGN: 64:: aIntern
#endif

#ifdef WITH_OPENMP
542
        if (posix_memalign(aIntern_ptr, 64_c_intptr_t, stripe_width*a_dim2*stripe_count*max_threads*     &
543
               C_SIZEOF(a_var)) /= 0) then
544 545
          print *,"trans_ev_tridi_to_band_&
          &MATH_DATATYPE&
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
546
          &: error when allocating aIntern"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
547
          stop 1
548 549 550 551 552 553
        endif

        call c_f_pointer(aIntern_ptr, aIntern, [stripe_width,a_dim2,stripe_count,max_threads])
        ! allocate(aIntern(stripe_width,a_dim2,stripe_count,max_threads), stat=istat, errmsg=errorMessage)

        ! aIntern(:,:,:,:) should be set to 0 in a parallel region, not here!
554 555 556

#else /* WITH_OPENMP */

557
        if (posix_memalign(aIntern_ptr, 64_c_intptr_t, stripe_width*a_dim2*stripe_count*  &
558
            C_SIZEOF(a_var)) /= 0) then
559
          print *,"trans_ev_tridi_to_band_real: error when allocating aIntern"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
560
          stop 1
561 562 563 564 565
        endif

        call c_f_pointer(aIntern_ptr, aIntern,[stripe_width,a_dim2,stripe_count] )
        !allocate(aIntern(stripe_width,a_dim2,stripe_count), stat=istat, errmsg=errorMessage)

566
        aIntern(:,:,:) = 0.0_rck
567 568 569 570 571 572 573 574
#endif /* WITH_OPENMP */
      endif !useGPU

      allocate(row(l_nev), stat=istat, errmsg=errorMessage)
      if (istat .ne. 0) then
        print *,"trans_ev_tridi_to_band_&
        &MATH_DATATYPE&
        &: error when allocating row"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
575
        stop 1
576 577
      endif

578
      row(:) = 0.0_rck
579 580 581 582 583 584 585 586 587 588 589 590

      ! Copy q from a block cyclic distribution into a distribution with contiguous rows,
      ! and transpose the matrix using stripes of given stripe_width for cache blocking.

      ! The peculiar way it is done below is due to the fact that the last row should be
      ! ready first since it is the first one to start below

#ifdef WITH_OPENMP
      ! Please note about the OMP usage below:
      ! This is not for speed, but because we want the matrix a in the memory and
      ! in the cache of the correct thread (if possible)

591
      call obj%timer%start("OpenMP parallel" // PRECISION_SUFFIX)
592 593
      !$omp parallel do private(my_thread), schedule(static, 1)
      do my_thread = 1, max_threads
594
        aIntern(:,:,:,my_thread) = 0.0_rck ! if possible, do first touch allocation!
595 596 597
      enddo
      !$omp end parallel do

598
      call obj%timer%stop("OpenMP parallel" // PRECISION_SUFFIX)
599 600 601 602 603 604 605 606 607 608 609 610 611
#endif /* WITH_OPENMP */

      do ip = np_rows-1, 0, -1
        if (my_prow == ip) then
          ! Receive my rows which have not yet been received
          src_offset = local_index(limits(ip), my_prow, np_rows, nblk, -1)
          do i=limits(ip)+1,limits(ip+1)
            src = mod((i-1)/nblk, np_rows)

            if (src < my_prow) then
#ifdef WITH_OPENMP

#ifdef WITH_MPI
612
              if (wantDebug) call obj%timer%start("mpi_communication")
613
              call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
614
                            src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
615
              if (wantDebug) call obj%timer%stop("mpi_communication")
616 617 618 619 620 621 622

#else /* WITH_MPI */

!              row(1:l_nev) = row(1:l_nev)

#endif /* WITH_MPI */

623
              call obj%timer%start("OpenMP parallel" // PRECISION_SUFFIX)
624 625 626

!$omp parallel do private(my_thread), schedule(static, 1)
              do my_thread = 1, max_threads
627 628 629 630
                call unpack_row_&
                     &MATH_DATATYPE&
                     &_cpu_openmp_&
                     &PRECISION &
631
                                  (obj,aIntern, row, i-limits(ip), my_thread, stripe_count, &
632 633
                                   thread_width, stripe_width, l_nev)

634 635 636
              enddo
!$omp end parallel do

637
              call obj%timer%stop("OpenMP parallel" // PRECISION_SUFFIX)
638 639 640 641

#else /* WITH_OPENMP */
              if (useGPU) then
                ! An unpacking of the current row group may occur before queuing the next row
642 643 644 645
                call unpack_and_prepare_row_group_&
                &MATH_DATATYPE&
                &_gpu_&
                &PRECISION &
Andreas Marek's avatar
Retab  
Andreas Marek committed
646 647
                              ( &
                              row_group, row_group_dev, aIntern_dev, stripe_count, &
Andreas Marek's avatar
Andreas Marek committed
648 649 650
                                          stripe_width, last_stripe_width, a_dim2, l_nev,&
                                          row_group_size, nblk, unpack_idx, &
                                           i - limits(ip), .false.)
651
#ifdef WITH_MPI
652
                if (wantDebug) call obj%timer%start("mpi_communication")
653
                call MPI_Recv(row_group(:, row_group_size), l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
654
                              src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
655
                if (wantDebug) call obj%timer%stop("mpi_communication")
656 657 658 659 660 661 662

#else /* WITH_MPI */
                row_group(1:l_nev, row_group_size) = row(1:l_nev) ! is this correct?
#endif /* WITH_MPI */

              else ! useGPU
#ifdef WITH_MPI
663
                if (wantDebug) call obj%timer%start("mpi_communication")
664
                call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
665
                              src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
666
                if (wantDebug) call obj%timer%stop("mpi_communication")
667 668 669 670 671 672 673

#else /* WITH_MPI */

!                row(1:l_nev) = row(1:l_nev)

#endif /* WITH_MPI */

674
                call unpack_row_&
675 676 677
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
678
                                (obj,aIntern, row,i-limits(ip), stripe_count, stripe_width, last_stripe_width)
679 680 681 682 683 684 685 686 687 688
              endif ! useGPU
#endif /* WITH_OPENMP */

            elseif (src == my_prow) then

              src_offset = src_offset+1

              if (useGPU) then
#ifndef WITH_OPENMP

Andreas Marek's avatar
Andreas Marek committed
689 690
                 ! An unpacking of the current row group may occur before queuing the next row
                 call unpack_and_prepare_row_group_&
691 692 693
                      &MATH_DATATYPE&
                      &_gpu_&
                      &PRECISION &
Andreas Marek's avatar
Retab  
Andreas Marek committed
694
                  ( &
Andreas Marek's avatar
Andreas Marek committed
695 696 697 698 699
                               row_group, row_group_dev, aIntern_dev, stripe_count, &
                               stripe_width, last_stripe_width, a_dim2, l_nev,&
                               row_group_size, nblk, unpack_idx, &
                               i - limits(ip), .false.)

700
                row_group(:, row_group_size) = q(src_offset, 1:l_nev)
701
#else /* WITH_OPENMP */
Andreas Marek's avatar
Andreas Marek committed
702

703 704 705 706 707 708 709 710 711
!#if COMPLEXCASE == 1
!! why is an cuda call in the openmp region?
!                call unpack_and_prepare_row_group_complex_gpu_&
!                     &PRECISION&
!                     &(row_group, row_group_dev, aIntern_dev, stripe_count, stripe_width, &
!                      last_stripe_width, a_dim2, l_nev, row_group_size, nblk,      &
!                      unpack_idx, i - limits(ip),.false.)
!                      row_group(:, row_group_size) = q(src_offset, 1:l_nev)
!#endif
712 713 714 715 716 717 718

#endif /* not OpenMP */
              else
                row(:) = q(src_offset, 1:l_nev)
              endif

#ifdef WITH_OPENMP
719
              call obj%timer%start("OpenMP parallel" // PRECISION_SUFFIX)
720 721 722

!$omp parallel do private(my_thread), schedule(static, 1)
              do my_thread = 1, max_threads
723 724 725 726
                call unpack_row_&
                     &MATH_DATATYPE&
                     &_cpu_openmp_&
                     &PRECISION &
727
                                   (obj,aIntern, row, i-limits(ip), my_thread, stripe_count, thread_width, stripe_width, l_nev)
728

729 730 731
              enddo
!$omp end parallel do

732
              call obj%timer%stop("OpenMP parallel" // PRECISION_SUFFIX)
733 734 735 736 737 738

#else /* WITH_OPENMP */

              if (useGPU) then

              else
739 740 741 742
                call unpack_row_&
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
743
                                (obj,aIntern, row,i-limits(ip),  stripe_count, stripe_width, last_stripe_width)
744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759
              endif

#endif /* WITH_OPENMP */

            endif
          enddo

          ! Send all rows which have not yet been send
          src_offset = 0
          do dst = 0, ip-1
            do i=limits(dst)+1,limits(dst+1)
              if (mod((i-1)/nblk, np_rows) == my_prow) then
                src_offset = src_offset+1
                row(:) = q(src_offset, 1:l_nev)

#ifdef WITH_MPI
760
                if (wantDebug) call obj%timer%start("mpi_communication")
761
                call MPI_Send(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
762
                              dst, 0, mpi_comm_rows, mpierr)
763
                if (wantDebug) call obj%timer%stop("mpi_communication")
764 765 766 767 768 769
#endif /* WITH_MPI */
              endif
            enddo
          enddo

        else if (my_prow < ip) then
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
770

771 772 773 774 775 776 777 778
          ! Send all rows going to PE ip
          src_offset = local_index(limits(ip), my_prow, np_rows, nblk, -1)
          do i=limits(ip)+1,limits(ip+1)
            src = mod((i-1)/nblk, np_rows)
            if (src == my_prow) then
              src_offset = src_offset+1
              row(:) = q(src_offset, 1:l_nev)
#ifdef WITH_MPI
779
              if (wantDebug) call obj%timer%start("mpi_communication")
780
              call MPI_Send(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
781
                            ip, 0, mpi_comm_rows, mpierr)
782
              if (wantDebug) call obj%timer%stop("mpi_communication")
783 784 785 786 787 788 789 790 791 792 793
#endif /* WITH_MPI */
            endif
          enddo

          ! Receive all rows from PE ip
          do i=limits(my_prow)+1,limits(my_prow+1)
            src = mod((i-1)/nblk, np_rows)
            if (src == ip) then
#ifdef WITH_OPENMP

#ifdef WITH_MPI
794
              if (wantDebug) call obj%timer%start("mpi_communication")
795
              call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
796
                            src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
797
              if (wantDebug) call obj%timer%stop("mpi_communication")
798 799 800 801 802 803
#else /* WITH_MPI */

!              row(1:l_nev) = row(1:l_nev)

#endif /* WITH_MPI */

804
              call obj%timer%start("OpenMP parallel" // PRECISION_SUFFIX)
805 806
!$omp parallel do private(my_thread), schedule(static, 1)
              do my_thread = 1, max_threads
807 808 809 810
                call unpack_row_&
                     &MATH_DATATYPE&
                     &_cpu_openmp_&
                     &PRECISION &
811
                                 (obj,aIntern, row, i-limits(my_prow), my_thread, stripe_count, thread_width, stripe_width, l_nev)
812 813
              enddo
!$omp end parallel do
814
              call obj%timer%stop("OpenMP parallel" // PRECISION_SUFFIX)
815 816 817 818

#else /* WITH_OPENMP */
              if (useGPU) then
                ! An unpacking of the current row group may occur before queuing the next row
Andreas Marek's avatar
Andreas Marek committed
819
                call unpack_and_prepare_row_group_&
820 821 822 823
                     &MATH_DATATYPE&
                     &_gpu_&
                     &PRECISION&
                     &( &
Andreas Marek's avatar
Andreas Marek committed
824 825 826 827
                  row_group, row_group_dev, aIntern_dev, stripe_count,  &
                  stripe_width, last_stripe_width, a_dim2, l_nev,       &
                  row_group_size, nblk, unpack_idx,                     &
                  i - limits(my_prow), .false.)
828 829

#ifdef WITH_MPI
830
               if (wantDebug) call obj%timer%start("mpi_communication")
831
               call MPI_Recv(row_group(:, row_group_size), l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
832
                             src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
833
               if (wantDebug) call obj%timer%stop("mpi_communication")
834 835 836 837 838 839 840
#else /* WITH_MPI */

               row_group(1:l_nev,row_group_size) = row(1:l_nev) ! is this correct ?
#endif /* WITH_MPI */

              else ! useGPU
#ifdef WITH_MPI
841
                if (wantDebug) call obj%timer%start("mpi_communication")
842
                call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
843
                              src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
844
                if (wantDebug) call obj%timer%stop("mpi_communication")
845 846 847 848 849
#else /* WITH_MPI */

!                row(1:l_nev) = row(1:l_nev)

#endif
850
                call unpack_row_&
851 852 853
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
854
                                (obj,aIntern, row,i-limits(my_prow), stripe_count, stripe_width, last_stripe_width)
855 856 857 858 859 860 861 862 863 864 865
              endif ! useGPU

#endif /* WITH_OPENMP */

            endif
          enddo
        endif
      enddo

      if (useGPU) then
        ! Force an unpacking of all remaining rows that haven't been unpacked yet
Andreas Marek's avatar
Andreas Marek committed
866
        call unpack_and_prepare_row_group_&
867 868 869 870
             &MATH_DATATYPE&
             &_gpu_&
             &PRECISION&
             &( &
Andreas Marek's avatar
Andreas Marek committed
871 872 873 874 875
          row_group, row_group_dev, aIntern_dev, stripe_count, &
          stripe_width, last_stripe_width, &
          a_dim2, l_nev, row_group_size, nblk, unpack_idx,     &
          -1, .true.)

876 877 878 879 880 881
        successCUDA = cuda_devicesynchronize()

         if (.not.(successCUDA)) then
           print *,"trans_ev_tridi_to_band_&
           &MATH_DATATYPE&
           &: error in cudaDeviceSynchronize"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
882
           stop 1
883 884 885 886 887 888 889 890 891 892 893 894 895
         endif
      endif

      ! Set up result buffer queue

      num_result_blocks = ((na-1)/nblk + np_rows - my_prow) / np_rows

      num_result_buffers = 4*nfact
      allocate(result_buffer(l_nev,nblk,num_result_buffers), stat=istat, errmsg=errorMessage)
      if (istat .ne. 0) then
        print *,"trans_ev_tridi_to_band_&
        &MATH_DATATYPE&
        &: error when allocating result_buffer"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
896
        stop 1
897 898 899 900 901 902 903
      endif

      allocate(result_send_request(num_result_buffers), stat=istat, errmsg=errorMessage)
      if (istat .ne. 0) then
        print *,"trans_ev_tridi_to_band_&
        &MATH_DATATYPE&
        &: error when allocating result_send_request"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
904
        stop 1
905 906 907 908 909 910 911
      endif

      allocate(result_recv_request(num_result_buffers), stat=istat, errmsg=errorMessage)
      if (istat .ne. 0) then
        print *,"trans_ev_tridi_to_band_&
        &MATH_DATATYPE&
        &: error when allocating result_recv_request"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
912
        stop 1
913 914 915 916 917 918 919 920 921
      endif

#ifdef WITH_MPI
      result_send_request(:) = MPI_REQUEST_NULL
      result_recv_request(:) = MPI_REQUEST_NULL
#endif

      ! Queue up buffers
#ifdef WITH_MPI
922
      if (wantDebug) call obj%timer%start("mpi_communication")
923 924 925

      if (my_prow > 0 .and. l_nev>0) then ! note: row 0 always sends
        do j = 1, min(num_result_buffers, num_result_blocks)
926
          call MPI_Irecv(result_buffer(1,1,j), l_nev*nblk, MPI_MATH_DATATYPE_PRECISION_EXPL,     &
927 928 929
                         0, result_recv_tag, mpi_comm_rows, result_recv_request(j), mpierr)
        enddo
      endif
930
      if (wantDebug) call obj%timer%stop("mpi_communication")