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,         &
60
     hh_trans, mpi_comm_rows, mpi_comm_cols, wantDebug, useGPU, 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 99
      use cuda_functions
      use precision
      use iso_c_binding
      implicit none
100
#include "../general/precision_kinds.F90"
101
      class(elpa_abstract_impl_t), intent(inout) :: obj
Andreas Marek's avatar
Andreas Marek committed
102
      logical, intent(in)                        :: useGPU
103

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

#ifdef USE_ASSUMED_SIZE
108
      MATH_DATATYPE(kind=rck)                    :: q(ldq,*)
109
#else
110
      MATH_DATATYPE(kind=rck)                    :: q(ldq,matrixCols)
111 112
#endif

113
      MATH_DATATYPE(kind=rck), intent(in)        :: hh_trans(:,:)
Andreas Marek's avatar
Andreas Marek committed
114
      integer(kind=c_intptr_t)                   :: q_dev
Andreas Marek's avatar
Andreas Marek committed
115

Andreas Marek's avatar
Andreas Marek committed
116 117 118 119 120 121
      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
122
#ifdef WITH_OPENMP
Andreas Marek's avatar
Andreas Marek committed
123
      integer(kind=ik)                           :: thread_width, csw, b_off, b_len
124
#endif
Andreas Marek's avatar
Andreas Marek committed
125 126 127
      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
128

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

Andreas Marek's avatar
Andreas Marek committed
137
      type(c_ptr)                                :: aIntern_ptr
138

139 140
      MATH_DATATYPE(kind=rck)   , allocatable    :: row(:)
      MATH_DATATYPE(kind=rck)   , allocatable    :: row_group(:,:)
141 142

#ifdef WITH_OPENMP
143 144
      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(:,:)
145
#else
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 149
#endif

150 151
      integer(kind=c_intptr_t)                   :: aIntern_dev
      integer(kind=c_intptr_t)                   :: bcast_buffer_dev
152
      integer(kind=c_intptr_t)                   :: num
Andreas Marek's avatar
Andreas Marek committed
153
      integer(kind=c_intptr_t)                   :: dev_offset, dev_offset_1
154 155 156 157 158
      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
159

160 161
      integer(kind=ik)                           :: n_times
      integer(kind=ik)                           :: top, chunk, this_chunk
162

163 164
      MATH_DATATYPE(kind=rck), allocatable       :: result_buffer(:,:,:)
      MATH_DATATYPE(kind=rck), allocatable       :: bcast_buffer(:,:)
165

166
      integer(kind=ik)                           :: n_off
167

168 169 170
      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(:)
171 172 173

      ! MPI send/recv tags, arbitrary

174 175 176
      integer(kind=ik), parameter                :: bottom_recv_tag = 111
      integer(kind=ik), parameter                :: top_recv_tag    = 222
      integer(kind=ik), parameter                :: result_recv_tag = 333
177
#ifdef WITH_OPENMP
178 179
      integer(kind=ik)                           :: max_threads, my_thread
      integer(kind=ik)                           :: omp_get_max_threads
180 181 182 183
#endif


      ! Just for measuring the kernel performance
184
      real(kind=c_double)                        :: kernel_time, kernel_time_recv ! MPI_WTIME always needs double
185
      ! long integer
186
      integer(kind=lik)                          :: kernel_flops, kernel_flops_recv
187

188 189 190 191 192
      logical, intent(in)                        :: wantDebug
      logical                                    :: success
      integer(kind=ik)                           :: istat, print_flops
      character(200)                             :: errorMessage
      logical                                    :: successCUDA
193
#ifndef WITH_MPI
194
      integer(kind=ik)                           :: j1
195
#endif
196 197 198 199
      integer(kind=c_intptr_t), parameter        :: size_of_datatype = size_of_&
                                                                     &PRECISION&
                                                                     &_&
                                                                     &MATH_DATATYPE
200

201
      call obj%timer%start("trans_ev_tridi_to_band_&
202
      &MATH_DATATYPE&
Andreas Marek's avatar
Cleanup  
Andreas Marek committed
203
      &" // &
204 205 206
      &PRECISION_SUFFIX &
      )

207
      n_times = 0
208 209 210 211 212 213 214 215 216 217 218 219 220
      if (useGPU) then
        unpack_idx = 0
        row_group_size = 0
      endif

      success = .true.
      kernel_time = 0.0
      kernel_flops = 0

#ifdef WITH_OPENMP
      max_threads = 1
      max_threads = omp_get_max_threads()
#endif
221
      if (wantDebug) call obj%timer%start("mpi_communication")
222 223 224 225
      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)
226
      if (wantDebug) call obj%timer%stop("mpi_communication")
227 228 229 230 231

      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_&
232 233
                                &MATH_DATATYPE&
                                &: ERROR: nbw=',nbw,', nblk=',nblk
234
            write(error_unit,*) 'ELPA2_trans_ev_tridi_to_band_&
235 236
                                &MATH_DATATYPE&
                                &: band backtransform works only for nbw==n*nblk'
237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
          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
262
        ! Suggested stripe width is 48 - should this be reduced for the complex case ???
263 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

        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
294
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
295
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
296
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
297 298 299 300

            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
301
    else
302
            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
303
                                            ! (4 * sizeof(double) == 32)
304
          endif
305
#else
306
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
307
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
308
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
309 310 311


            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
312
                                               ! (16 * sizeof(float) == 64)
313

Andreas Marek's avatar
Retab  
Andreas Marek committed
314
    else
315
            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
316
                                            ! (8 * sizeof(float) == 32)
317
          endif
318 319 320 321 322
#endif
#endif /* REALCASE */

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

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
326
            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
327
                                            ! (4 * sizeof(double complex) == 64)
328

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

            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
332 333
                                            ! (2 * sizeof(double complex) == 32)
    endif
334
#else
335

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

            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
340
                                            ! (8 * sizeof(float complex) == 64)
341 342 343

          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
344 345
                                            ! (4 * sizeof(float complex) == 32)
    endif
346 347 348
#endif
#endif /* COMPLEXCASE */

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

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

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

361 362
        ! Suggested stripe width is 48 since 48*64 real*8 numbers should fit into
        ! every primary cache
363
        ! Suggested stripe width is 48 - should this be reduced for the complex case ???
364 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

        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
394
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
395
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
396
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
397 398 399 400

            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
401
    else
402
            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
403
                                            ! (4 * sizeof(double) == 32)
404
          endif
405
#else
406
          if (kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK2 .or. &
407
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK4 .or. &
408
              kernel .eq. ELPA_2STAGE_REAL_AVX512_BLOCK6) then
409 410 411


            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
412
                                               ! (16 * sizeof(float) == 64)
413

Andreas Marek's avatar
Retab  
Andreas Marek committed
414
    else
415
            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
416
                                            ! (8 * sizeof(float) == 32)
417
          endif
418 419 420 421 422
#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
423

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

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
427
            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
428
                                            ! (4 * sizeof(double complex) == 64)
429

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

            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
433 434
                                            ! (2 * sizeof(double complex) == 32)
    endif
435
#else
436

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

Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
440
            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
441
                                            ! (8 * sizeof(float complex) == 64)
442 443 444

          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
445 446
                                            ! (4 * sizeof(float complex) == 32)
    endif
447 448 449 450 451 452 453 454 455 456 457 458 459 460
#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_&
461 462
                &MATH_DATATYPE&
                &: error when allocating limits"//errorMessage
Andreas Marek's avatar
Andreas Marek committed
463
        stop 1
464
      endif
465
      call determine_workload(obj,na, nbw, np_rows, limits)
466 467 468 469 470 471

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

      a_dim2 = max_blk_size + nbw

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

489
        num =  (l_nev)* size_of_datatype
490 491 492 493 494
        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
495
          stop 1
496 497 498 499 500 501 502
        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
503
          stop 1
504 505 506 507 508 509 510 511
        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
512
          stop 1
513 514
        endif

515
        row_group(:, :) = 0.0_rck
516
        num =  (l_nev*nblk)* size_of_datatype
517 518 519 520 521
        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
522
          stop 1
523 524 525 526 527 528 529
        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
530
          stop 1
531 532 533 534 535 536 537 538 539 540
        endif

      else ! GPUs are not used

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

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

#else /* WITH_OPENMP */

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

565
        aIntern(:,:,:) = 0.0_rck
566 567 568 569 570 571 572 573
#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
574
        stop 1
575 576
      endif

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

      ! 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)

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

597
      call obj%timer%stop("OpenMP parallel" // PRECISION_SUFFIX)
598 599 600 601 602 603 604 605 606 607 608 609 610
#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
611
              if (wantDebug) call obj%timer%start("mpi_communication")
612
              call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
613
                            src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
614
              if (wantDebug) call obj%timer%stop("mpi_communication")
615 616 617 618 619 620 621

#else /* WITH_MPI */

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

#endif /* WITH_MPI */

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

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

633 634 635
              enddo
!$omp end parallel do

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

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

#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
662
                if (wantDebug) call obj%timer%start("mpi_communication")
663
                call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
664
                              src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
665
                if (wantDebug) call obj%timer%stop("mpi_communication")
666 667 668 669 670 671 672

#else /* WITH_MPI */

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

#endif /* WITH_MPI */

673
                call unpack_row_&
674 675 676
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
677
                                (obj,aIntern, row,i-limits(ip), stripe_count, stripe_width, last_stripe_width)
678 679 680 681 682 683 684 685 686 687
              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
688 689
                 ! An unpacking of the current row group may occur before queuing the next row
                 call unpack_and_prepare_row_group_&
690 691 692
                      &MATH_DATATYPE&
                      &_gpu_&
                      &PRECISION &
Andreas Marek's avatar
Retab  
Andreas Marek committed
693
                  ( &
Andreas Marek's avatar
Andreas Marek committed
694 695 696 697 698
                               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.)

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

702 703 704 705 706 707 708 709 710
!#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
711 712 713 714 715 716 717

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

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

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

728 729 730
              enddo
!$omp end parallel do

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

#else /* WITH_OPENMP */

              if (useGPU) then

              else
738 739 740 741
                call unpack_row_&
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
742
                                (obj,aIntern, row,i-limits(ip),  stripe_count, stripe_width, last_stripe_width)
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758
              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
759
                if (wantDebug) call obj%timer%start("mpi_communication")
760
                call MPI_Send(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
761
                              dst, 0, mpi_comm_rows, mpierr)
762
                if (wantDebug) call obj%timer%stop("mpi_communication")
763 764 765 766 767 768
#endif /* WITH_MPI */
              endif
            enddo
          enddo

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

770 771 772 773 774 775 776 777
          ! 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
778
              if (wantDebug) call obj%timer%start("mpi_communication")
779
              call MPI_Send(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
780
                            ip, 0, mpi_comm_rows, mpierr)
781
              if (wantDebug) call obj%timer%stop("mpi_communication")
782 783 784 785 786 787 788 789 790 791 792
#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
793
              if (wantDebug) call obj%timer%start("mpi_communication")
794
              call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
795
                            src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
796
              if (wantDebug) call obj%timer%stop("mpi_communication")
797 798 799 800 801 802
#else /* WITH_MPI */

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

#endif /* WITH_MPI */

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

#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
818
                call unpack_and_prepare_row_group_&
819 820 821 822
                     &MATH_DATATYPE&
                     &_gpu_&
                     &PRECISION&
                     &( &
Andreas Marek's avatar
Andreas Marek committed
823 824 825 826
                  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.)
827 828

#ifdef WITH_MPI
829
               if (wantDebug) call obj%timer%start("mpi_communication")
830
               call MPI_Recv(row_group(:, row_group_size), l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
831
                             src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
832
               if (wantDebug) call obj%timer%stop("mpi_communication")
833 834 835 836 837 838 839
#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
840
                if (wantDebug) call obj%timer%start("mpi_communication")
841
                call MPI_Recv(row, l_nev, MPI_MATH_DATATYPE_PRECISION_EXPL, &
842
                              src, 0, mpi_comm_rows, MPI_STATUS_IGNORE, mpierr)
843
                if (wantDebug) call obj%timer%stop("mpi_communication")
844 845 846 847 848
#else /* WITH_MPI */

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

#endif
849
                call unpack_row_&
850 851 852
                     &MATH_DATATYPE&
                     &_cpu_&
                     &PRECISION &
853
                                (obj,aIntern, row,i-limits(my_prow), stripe_count, stripe_width, last_stripe_width)
854 855 856 857 858 859 860 861 862 863 864
              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
865
        call unpack_and_prepare_row_group_&
866 867 868 869
             &MATH_DATATYPE&
             &_gpu_&
             &PRECISION&
             &( &
Andreas Marek's avatar
Andreas Marek committed
870 871 872 873 874
          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.)

875 876 877 878 879 880
        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
881
           stop 1
882 883 884 885 886 887 888 889 890 891 892 893 894
         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
895
        stop 1
896 897 898 899 900 901 902
      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
903
        stop 1
904 905 906 907 908 909 910
      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
911
        stop 1
912 913 914 915 916 917 918 919 920
      endif

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

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

      if (my_prow > 0 .and. l_nev>0) then ! note: row 0 always sends
        do j = 1, min(num_result_buffers, num_result_blocks)
925
          call MPI_Irecv(result_buffer(1,1,j), l_nev*nblk, MPI_MATH_DATATYPE_PRECISION_EXPL,     &
926 927 928
                         0, result_recv_tag, mpi_comm_rows, result_recv_request(j), mpierr)
        enddo
      endif
929
      if (wantDebug) call obj%timer%stop("mpi_communication")
930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945
#else /* WITH_MPI */

      ! carefull the "recv" has to be done at the corresponding wait or send
      ! result_buffer(1: l_nev*nblk,1,j) =result_buffer(1:l_nev*nblk,1,nbuf)

#endif /* WITH_MPI */

      num_bufs_recvd = 0 ! No buffers received yet

      ! Initialize top/bottom requests

      allocate(top_send_request(stripe_count), stat<