elpa_impl.F90 47.3 KB
Newer Older
1
2
3
!
!    Copyright 2017, L. Hüdepohl and A. Marek, MPCDF
!
Andreas Marek's avatar
Andreas Marek committed
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
!    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.
!
48
#include "config-f90.h"
49

50
51
module elpa_impl
  use elpa_api
52
  use, intrinsic :: iso_c_binding
53
  implicit none
54

55
56
  private
  public :: elpa_impl_allocate
57

58
  type, extends(elpa_t) :: elpa_impl_t
Andreas Marek's avatar
Andreas Marek committed
59
   private
60
   type(c_ptr)         :: index = C_NULL_PTR
61

62
63
64
65
   contains
     procedure, public :: setup => elpa_setup

     procedure, public :: get => elpa_get_integer
66
67
68
69
70
71
     procedure, public :: get_double => elpa_get_double
     procedure, public :: is_set => elpa_is_set
     procedure, public :: print_options => elpa_print_options

     procedure, public :: get_real_kernel => elpa_get_real_kernel
     procedure, public :: get_complex_kernel => elpa_get_complex_kernel
72
73
74
75
76

     procedure, public :: destroy => elpa_destroy

     ! privates:
     procedure, private :: elpa_set_integer
77
     procedure, private :: elpa_set_double
78

79
     procedure, private :: elpa_solve_real_double
80
     procedure, private :: elpa_solve_real_single
81
     procedure, private :: elpa_solve_complex_double
82
     procedure, private :: elpa_solve_complex_single
83

84
85
86
87
     procedure, private :: elpa_multiply_at_b_double
     procedure, private :: elpa_multiply_at_b_single
     procedure, private :: elpa_multiply_ah_b_double
     procedure, private :: elpa_multiply_ah_b_single
88

89
90
91
92
     procedure, private :: elpa_cholesky_double_real
     procedure, private :: elpa_cholesky_single_real
     procedure, private :: elpa_cholesky_double_complex
     procedure, private :: elpa_cholesky_single_complex
93
94
95
96
97

     procedure, private :: elpa_invert_trm_double_real
     procedure, private :: elpa_invert_trm_single_real
     procedure, private :: elpa_invert_trm_double_complex
     procedure, private :: elpa_invert_trm_single_complex
98
99
100

     procedure, private :: elpa_solve_tridi_double_real
     procedure, private :: elpa_solve_tridi_single_real
101
102
103

     procedure, private :: associate_int => elpa_associate_int

104
  end type elpa_impl_t
105
106
107

  contains

108
    function elpa_impl_allocate(error) result(obj)
Andreas Marek's avatar
Andreas Marek committed
109
110
      use precision
      use elpa_utilities, only : error_unit
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
111
      use elpa_generated_fortran_interfaces
Andreas Marek's avatar
Andreas Marek committed
112

113
114
115
116
      type(elpa_impl_t), pointer   :: obj
      integer, optional            :: error

      allocate(obj)
Andreas Marek's avatar
Andreas Marek committed
117

Andreas Marek's avatar
Andreas Marek committed
118
119
      ! check whether init has ever been called
      if (.not.(elpa_initialized())) then
120
        write(error_unit, *) "elpa_allocate(): you must call elpa_init() once before creating instances of ELPA"
121
122
        if(present(error)) then
          error = ELPA_ERROR
123
        endif
Andreas Marek's avatar
Andreas Marek committed
124
125
        return
      endif
Andreas Marek's avatar
Andreas Marek committed
126

127
128
129
      obj%index = elpa_index_instance()

      ! Associate some important integer pointers for convenience
130
131
132
133
134
135
136
137
      obj%na => obj%associate_int("na")
      obj%nev => obj%associate_int("nev")
      obj%local_nrows => obj%associate_int("local_nrows")
      obj%local_ncols => obj%associate_int("local_ncols")
      obj%nblk => obj%associate_int("nblk")

      if(present(error)) then
        error = ELPA_OK
138
139
      endif
    end function
Andreas Marek's avatar
Andreas Marek committed
140

141

142
    function elpa_setup(self) result(success)
143
      use elpa1_impl, only : elpa_get_communicators_impl
144
145
      class(elpa_impl_t), intent(inout) :: self
      integer :: success, success2
146
      integer :: mpi_comm_rows, mpi_comm_cols, mpierr
147

148
      success = ELPA_ERROR
149

150
151
152
      if (self%is_set("mpi_comm_parent") == ELPA_OK .and. &
          self%is_set("process_row") == ELPA_OK .and. &
          self%is_set("process_col") == ELPA_OK) then
153

154
155
156
157
158
159
        mpierr = elpa_get_communicators_impl(&
                        self%get("mpi_comm_parent"), &
                        self%get("process_row"), &
                        self%get("process_col"), &
                        mpi_comm_rows, &
                        mpi_comm_cols)
160

161
162
163
164
165
        call self%set("mpi_comm_rows", mpi_comm_rows)
        call self%set("mpi_comm_cols", mpi_comm_cols)

        success = ELPA_OK
      endif
166

167
168
169
      if (self%is_set("mpi_comm_rows") == ELPA_OK .and. self%is_set("mpi_comm_cols") == ELPA_OK) then
        success = ELPA_OK
      endif
170
171

    end function
172
173

    subroutine elpa_set_integer(self, name, value, success)
174
      use iso_c_binding
175
176
      use elpa_generated_fortran_interfaces
      use elpa_utilities, only : error_unit
177
      class(elpa_impl_t)              :: self
178
179
180
181
      character(*), intent(in)        :: name
      integer(kind=c_int), intent(in) :: value
      integer, optional               :: success
      integer                         :: actual_success
182

183
      actual_success = elpa_index_set_int_value(self%index, name // c_null_char, value)
184

185
186
187
188
      if (present(success)) then
        success = actual_success

      else if (actual_success /= ELPA_OK) then
189
190
        write(error_unit,'(a,a,a,i0,a)') "ELPA: Error setting option '", name, "' to value ", value, &
                " and you did not check for errors!"
191
192

      end if
193
194
    end subroutine

195

196
197
198
199
200
201
202
203
204
205
206
207
208
    subroutine elpa_set_integer_by_string(self, name, value, success)
      use iso_c_binding
      use elpa_generated_fortran_interfaces
      use elpa_utilities, only : error_unit
      class(elpa_impl_t)              :: self
      character(*), intent(in)        :: name
      character(*), intent(in)        :: value
      integer, optional               :: success
      integer                         :: actual_success

    end subroutine


209
    function elpa_get_integer(self, name, success) result(value)
210
      use iso_c_binding
211
      use elpa_generated_fortran_interfaces
212
      class(elpa_impl_t)             :: self
213
214
215
      character(*), intent(in)       :: name
      integer(kind=c_int)            :: value
      integer, intent(out), optional :: success
216

217
      value = elpa_index_get_int_value(self%index, name // c_null_char, success)
218

219
    end function
Andreas Marek's avatar
Andreas Marek committed
220

221
222

    function elpa_is_set(self, name) result(state)
223
224
      use iso_c_binding
      use elpa_generated_fortran_interfaces
225
      class(elpa_impl_t)       :: self
226
      character(*), intent(in) :: name
227
      integer                  :: state
228

229
      state = elpa_index_value_is_set(self%index, name // c_null_char)
230
231
    end function

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

    subroutine elpa_print_options(self, option_name, unit)
      use elpa_utilities, only : output_unit
      use elpa_generated_fortran_interfaces
      class(elpa_impl_t), intent(in) :: self
      character(kind=c_char, len=*), intent(in) :: option_name
      integer, intent(in), optional :: unit
      integer :: unit_actual, i, val

      if (present(unit)) then
        unit_actual = unit
      else
        unit_actual = output_unit
      endif

      do i = 0, elpa_index_cardinality_c(self%index, option_name // C_NULL_CHAR)
        val = elpa_index_enumerate_c(self%index, option_name // C_NULL_CHAR, i)
        if (elpa_index_int_is_valid_c(self%index, option_name // C_NULL_CHAR, val) == ELPA_OK) then
          write(unit_actual, '(a)') " " // elpa_c_string(elpa_index_int_value_to_string_helper_c(option_name // C_NULL_CHAR, val))
        end if
      end do
    end subroutine


256
    subroutine elpa_set_double(self, name, value, success)
Andreas Marek's avatar
Andreas Marek committed
257
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
258
      use elpa_generated_fortran_interfaces
259
      use elpa_utilities, only : error_unit
260
      class(elpa_impl_t)              :: self
261
      character(*), intent(in)        :: name
262
      real(kind=c_double), intent(in) :: value
263
264
      integer, optional               :: success
      integer                         :: actual_success
Andreas Marek's avatar
Andreas Marek committed
265

266
      actual_success = elpa_index_set_double_value(self%index, name // c_null_char, value)
Andreas Marek's avatar
Andreas Marek committed
267

268
269
270
271
      if (present(success)) then
        success = actual_success

      else if (actual_success /= ELPA_OK) then
272
        write(error_unit,'(a,a,es12.5,a)') "ELPA: Error setting option '", name, "' to value ", value, &
273
                " and you did not check for errors!"
274
275
276

      end if
    end subroutine
Andreas Marek's avatar
Andreas Marek committed
277
278


279
    function elpa_get_double(self, name, success) result(value)
Andreas Marek's avatar
Andreas Marek committed
280
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
281
      use elpa_generated_fortran_interfaces
282
      class(elpa_impl_t)             :: self
283
      character(*), intent(in)       :: name
284
      real(kind=c_double)            :: value
285
286
287
      integer, intent(out), optional :: success
      type(c_ptr) :: c_success_ptr

288
      value = elpa_index_get_double_value(self%index, name // c_null_char, success)
289
290

    end function
Andreas Marek's avatar
Andreas Marek committed
291
292


293
    function elpa_associate_int(self, name) result(value)
Andreas Marek's avatar
Andreas Marek committed
294
      use iso_c_binding
295
      use elpa_generated_fortran_interfaces
296
297
      use elpa_utilities, only : error_unit
      class(elpa_impl_t)             :: self
298
299
      character(*), intent(in)       :: name
      integer(kind=c_int), pointer   :: value
Andreas Marek's avatar
Andreas Marek committed
300

301
302
303
      type(c_ptr)                    :: value_p

      value_p = elpa_index_get_int_loc(self%index, name // c_null_char)
304
305
306
      if (.not. c_associated(value_p)) then
        write(error_unit, '(a,a,a)') "ELPA: Warning, received NULL pointer for entry '", name, "'"
      endif
307
308
      call c_f_pointer(value_p, value)
    end function
Andreas Marek's avatar
Andreas Marek committed
309

310
311

    subroutine elpa_solve_real_double(self, a, ev, q, success)
312
313
      use elpa2_impl
      use elpa1_impl
314
      use elpa_utilities, only : error_unit
315
      use precision
Andreas Marek's avatar
Andreas Marek committed
316
      use iso_c_binding
317
      class(elpa_impl_t)  :: self
Andreas Marek's avatar
Andreas Marek committed
318

319
320
321
#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
322
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
323
#endif
324
      real(kind=c_double) :: ev(self%na)
325
326

      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
327
328
      integer, optional   :: success
      integer(kind=c_int) :: success_internal
329
      logical             :: success_l, summary_timings
330

331
      logical             :: useGPU, useQR
332
      integer(kind=c_int) :: kernel
333

334
335
336
337
338
339
340
341
342
343
344
345
      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

        summary_timings = .true.
      else
        summary_timings = .false.
      endif


346
347
348
349
350
351
352
353
354
355
356
      if (self%get("gpu",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry gpu"
          stop
        endif

        useGPU = .true.
      else
        useGPU = .false.
      endif

357
      useQR = self%get("qr") == 1
358

359
360
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
        success_l = elpa_solve_evp_real_1stage_double_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
361
362
363
364
                                                           self%local_nrows,  self%nblk, self%local_ncols, &
                                                           self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
                                                           self%get("mpi_comm_parent"), useGPU, time_evp_fwd,     &
                                                           time_evp_solve, time_evp_back, summary_timings)
365

366
367
368
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
        kernel = self%get_real_kernel()
        success_l = elpa_solve_evp_real_2stage_double_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
369
370
371
372
373
                                                           self%local_nrows,  self%nblk, self%local_ncols, &
                                                           self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
                                                           self%get("mpi_comm_parent"), time_evp_fwd,     &
                                                           time_evp_solve, time_evp_back, summary_timings, useGPU, &
                                                           kernel, useQR)
374
375
376
377
      else
        print *,"unknown solver"
        stop
      endif
378
379
380
381
382
383
384
385
386
387

      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
Andreas Marek's avatar
Andreas Marek committed
388

389
390
391
392
393
394
      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

395
396
397
        call self%set("time_evp_fwd", time_evp_fwd)
        call self%set("time_evp_solve", time_evp_solve)
        call self%set("time_evp_back", time_evp_back)
398
399
      else

400
401
402
        call self%set("time_evp_fwd", -1.0_rk8)
        call self%set("time_evp_solve", -1.0_rk8)
        call self%set("time_evp_back", -1.0_rk8)
403
      endif
404
405
    end subroutine

406

407
    subroutine elpa_solve_real_single(self, a, ev, q, success)
408
409
      use elpa2_impl
      use elpa1_impl
410
      use elpa_utilities, only : error_unit
411
      use precision
412
      use iso_c_binding
413
      class(elpa_impl_t)  :: self
414
415
416
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
417
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
418
#endif
419
      real(kind=c_float)  :: ev(self%na)
420
421

      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
422
423
      integer, optional   :: success
      integer(kind=c_int) :: success_internal
424
      logical             :: success_l, summary_timings
425

426
      logical             :: useGPU, useQR
427
      integer(kind=c_int) :: kernel
428

429
#ifdef WANT_SINGLE_PRECISION_REAL
430
431
432
433
434
435
436
437
438
439
      if (self%get("timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

        summary_timings = .true.
      else
        summary_timings = .false.
      endif
440
441
442
443
444
445
446
447
448
449
450
451

      if (self%get("gpu",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry gpu"
          stop
        endif

        useGPU = .true.
      else
        useGPU = .false.
      endif

452
      useQR = self%get("qr") == 1
453

454
455
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
        success_l = elpa_solve_evp_real_1stage_single_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
456
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
457
458
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"), useGPU, time_evp_fwd,     &
459
                                                          time_evp_solve, time_evp_back, summary_timings)
460

461
462
463
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
        kernel = self%get_real_kernel()
        success_l = elpa_solve_evp_real_2stage_single_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
464
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
465
466
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"), time_evp_fwd,     &
467
                                                          time_evp_solve, time_evp_back, summary_timings, useGPU, &
468
                                                          kernel, useQR)
469
470
471
472
      else
        print *,"unknown solver"
        stop
      endif
473
474
475
476
477
478
479
480
481
482

      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
483
484
485
486
487
488
489
490


      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

491
492
493
        call self%set("time_evp_fwd", time_evp_fwd)
        call self%set("time_evp_solve", time_evp_solve)
        call self%set("time_evp_back", time_evp_back)
494
495
      else

496
497
498
        call self%set("time_evp_fwd", -1.0_rk8)
        call self%set("time_evp_solve", -1.0_rk8)
        call self%set("time_evp_back", -1.0_rk8)
499
      endif
500
#else
501
      print *,"This installation of the ELPA library has not been build with single-precision support"
502
503
504
505
      success = ELPA_ERROR
#endif
    end subroutine

506
507

    subroutine elpa_solve_complex_double(self, a, ev, q, success)
508
509
      use elpa2_impl
      use elpa1_impl
510
      use elpa_utilities, only : error_unit
511
      use precision
512
      use iso_c_binding
513
      class(elpa_impl_t)             :: self
514

515
516
517
#ifdef USE_ASSUMED_SIZE
      complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
518
      complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
519
#endif
520
      real(kind=c_double)            :: ev(self%na)
521

522
523
      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back

524
525
      integer, optional              :: success
      integer(kind=c_int)            :: success_internal
526
      logical                        :: success_l, summary_timings
527

528
      logical                        :: useGPU
529
      integer(kind=c_int) :: kernel
530
531
532
533
534
535
536
537
538
539
      if (self%get("timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

        summary_timings = .true.
      else
        summary_timings = .false.
      endif
540
541
542
543
544
545
546
547
548
549
550
551

      if (self%get("gpu",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry gpu"
          stop
        endif

        useGPU = .true.
      else
        useGPU = .false.
      endif

552
553
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
        success_l = elpa_solve_evp_complex_1stage_double_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
554
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
555
556
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"), useGPU, time_evp_fwd,     &
557
                                                          time_evp_solve, time_evp_back, summary_timings)
558

559
560
561
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
        kernel = self%get_complex_kernel()
        success_l = elpa_solve_evp_complex_2stage_double_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
562
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
563
564
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"), time_evp_fwd,     &
565
                                                          time_evp_solve, time_evp_back, summary_timings, useGPU, &
566
                                                          kernel)
567
568
569
570
      else
        print *,"unknown solver"
        stop
      endif
571
572
573
574
575
576
577
578
579
580
581

      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif

582
583
584
585
586
587
      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

588
589
590
        call self%set("time_evp_fwd", time_evp_fwd)
        call self%set("time_evp_solve", time_evp_solve)
        call self%set("time_evp_back", time_evp_back)
591
592
      else

593
594
595
        call self%set("time_evp_fwd", -1.0_rk8)
        call self%set("time_evp_solve", -1.0_rk8)
        call self%set("time_evp_back", -1.0_rk8)
596
      endif
597
598
599
    end subroutine


600
    subroutine elpa_solve_complex_single(self, a, ev, q, success)
601
602
      use elpa2_impl
      use elpa1_impl
603
604
605
      use elpa_utilities, only : error_unit

      use iso_c_binding
606
      use precision
607
      class(elpa_impl_t)            :: self
608
609
610
611
612
613
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck4)             :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
      complex(kind=ck4)             :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
#endif
      real(kind=rk4)                :: ev(self%na)
614

615
      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
616
617
      integer, optional             :: success
      integer(kind=c_int)           :: success_internal
618
      logical                       :: success_l, summary_timings
619

620
      logical                       :: useGPU
621
      integer(kind=c_int) :: kernel
622
#ifdef WANT_SINGLE_PRECISION_COMPLEX
623
624
625
626
627
628
629
630
631
632
633
634

      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

        summary_timings = .true.
      else
        summary_timings = .false.
      endif

635
636
637
638
639
640
641
642
643
644
645
      if (self%get("gpu",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry gpu"
          stop
        endif

        useGPU = .true.
      else
        useGPU = .false.
      endif

646
647
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
        success_l = elpa_solve_evp_complex_1stage_single_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
648
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
649
650
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"), useGPU, time_evp_fwd,     &
651
                                                          time_evp_solve, time_evp_back, summary_timings)
652

653
654
655
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
        kernel = self%get_complex_kernel()
        success_l = elpa_solve_evp_complex_2stage_single_impl(self, self%na, self%nev, a, self%local_nrows, ev, q,  &
656
                                                          self%local_nrows,  self%nblk, self%local_ncols, &
657
658
                                                          self%get("mpi_comm_rows"), self%get("mpi_comm_cols"),         &
                                                          self%get("mpi_comm_parent"),  time_evp_fwd,     &
659
                                                          time_evp_solve, time_evp_back, summary_timings, useGPU, &
660
                                                          kernel)
661
662
663
664
      else
        print *,"unknown solver"
        stop
      endif
665
666
667
668
669
670
671
672
673
674

      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
675
676
677
678
679
680
681

      if (self%get("summary_timings",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry summary timings"
          stop
        endif

682
683
684
        call self%set("time_evp_fwd", time_evp_fwd)
        call self%set("time_evp_solve", time_evp_solve)
        call self%set("time_evp_back", time_evp_back)
685
686
      else

687
688
689
        call self%set("time_evp_fwd", -1.0_rk8)
        call self%set("time_evp_solve", -1.0_rk8)
        call self%set("time_evp_back", -1.0_rk8)
690
691
      endif

692
#else
693
      print *,"This installation of the ELPA library has not been build with single-precision support"
694
695
696
697
      success = ELPA_ERROR
#endif
    end subroutine

698

699
700
701
    subroutine elpa_multiply_at_b_double (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                          c, ldc, ldcCols, success)
      use iso_c_binding
702
      use elpa1_auxiliary_impl
703
      use precision
704
      class(elpa_impl_t)              :: self
705
      character*1                     :: uplo_a, uplo_c
706
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
707
708
709
#ifdef USE_ASSUMED_SIZE
      real(kind=rk8)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
710
      real(kind=rk8)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
711
#endif
712
713
714
      integer, optional               :: success
      logical                         :: success_l

715
      success_l = elpa_mult_at_b_real_double_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
716
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
717
718
719
720
721
722
723
724
725
726
727
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
    end subroutine

728

729
730
731
    subroutine elpa_multiply_at_b_single (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                          c, ldc, ldcCols, success)
      use iso_c_binding
732
      use elpa1_auxiliary_impl
733
      use precision
734
      class(elpa_impl_t)              :: self
735
      character*1                     :: uplo_a, uplo_c
736
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
737
738
739
#ifdef USE_ASSUMED_SIZE
      real(kind=rk4)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
740
      real(kind=rk4)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
741
#endif
742
743
744
      integer, optional               :: success
      logical                         :: success_l
#ifdef WANT_SINGLE_PRECISION_REAL
745
      success_l = elpa_mult_at_b_real_single_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
746
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
747
748
749
750
751
752
753
754
755
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
756
757
758
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
      success = ELPA_ERROR
759
760
761
#endif
    end subroutine

762

763
764
765
    subroutine elpa_multiply_ah_b_double (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                          c, ldc, ldcCols, success)
      use iso_c_binding
766
      use elpa1_auxiliary_impl
767
      use precision
768
      class(elpa_impl_t)              :: self
769
      character*1                     :: uplo_a, uplo_c
770
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
771
772
773
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck8)               :: a(lda,*), b(ldb,*), c(ldc,*)
#else
774
      complex(kind=ck8)               :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
775
#endif
776
777
778
      integer, optional               :: success
      logical                         :: success_l

779
      success_l = elpa_mult_ah_b_complex_double_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
780
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
781
782
783
784
785
786
787
788
789
790
791
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
    end subroutine

792

793
794
795
    subroutine elpa_multiply_ah_b_single (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                          c, ldc, ldcCols, success)
      use iso_c_binding
796
      use elpa1_auxiliary_impl
797
      use precision
798
      class(elpa_impl_t)              :: self
799
      character*1                     :: uplo_a, uplo_c
800
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
801
802
803
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck4)               :: a(lda,*), b(ldb,*), c(ldc,*)
#else
804
      complex(kind=ck4)               :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
805
#endif
806
807
808
809
      integer, optional               :: success
      logical                         :: success_l

#ifdef WANT_SINGLE_PRECISION_COMPLEX
810
      success_l = elpa_mult_ah_b_complex_single_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
811
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
812
813
814
815
816
817
818
819
820
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
821
822
823
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
      success = ELPA_ERROR
824
825
826
#endif
    end subroutine

827

828
    subroutine elpa_cholesky_double_real (self, a, success)
829
      use iso_c_binding
830
      use elpa1_auxiliary_impl
831
      use precision
832
      class(elpa_impl_t)              :: self
833
834
835
#ifdef USE_ASSUMED_SIZE
      real(kind=rk8)                  :: a(self%local_nrows,*)
#else
836
      real(kind=rk8)                  :: a(self%local_nrows,self%local_ncols)
837
#endif
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
      integer, optional               :: success
      logical                         :: success_l
      integer(kind=c_int)             :: success_internal
      logical                         :: wantDebugIntern

      if (self%get("wantDebug",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry wantDebug"
          stop
        endif

        wantDebugIntern = .true.
      else
        wantDebugIntern = .false.
      endif

854

855
      success_l = elpa_cholesky_real_double_impl (self%na, a, self%local_nrows, self%nblk, &
856
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
857
858
859
860
861
862
863
864
865
866
867
868
                                                 wantDebugIntern)
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine

869

870
    subroutine elpa_cholesky_single_real(self, a, success)
871
      use iso_c_binding
872
      use elpa1_auxiliary_impl
873
      use precision
874
      class(elpa_impl_t)              :: self
875
876
877
#ifdef USE_ASSUMED_SIZE
      real(kind=rk4)                  :: a(self%local_nrows,*)
#else
878
      real(kind=rk4)                  :: a(self%local_nrows,self%local_ncols)
879
#endif
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
      integer, optional               :: success
      logical                         :: success_l
      integer(kind=c_int)             :: success_internal
      logical                         :: wantDebugIntern

      if (self%get("wantDebug",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry wantDebug"
          stop
        endif

        wantDebugIntern = .true.
      else
        wantDebugIntern = .false.
      endif

#if WANT_SINGLE_PRECISION_REAL
897
      success_l = elpa_cholesky_real_single_impl (self%na, a, self%local_nrows, self%nblk, &
898
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
899
                                                 wantDebugIntern)
900
901
902
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
      success = ELPA_ERROR
903
904
905
906
907
908
909
910
911
912
913
914
#endif
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine

915

916
    subroutine elpa_cholesky_double_complex (self, a, success)
917
      use iso_c_binding
918
      use elpa1_auxiliary_impl
919
      use precision
920
      class(elpa_impl_t)              :: self
921
922
923
924
925
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck8)               :: a(self%local_nrows,*)
#else
      complex(kind=ck8)               :: a(self%local_nrows,self%local_ncols)
#endif
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
      integer, optional               :: success
      logical                         :: success_l
      integer(kind=c_int)             :: success_internal
      logical                         :: wantDebugIntern

      if (self%get("wantDebug",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry wantDebug"
          stop
        endif

        wantDebugIntern = .true.
      else
        wantDebugIntern = .false.
      endif

942
      success_l = elpa_cholesky_complex_double_impl (self%na, a, self%local_nrows, self%nblk, &
943
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
944
945
946
947
948
949
950
951
952
953
954
955
                                                 wantDebugIntern)
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine

956

957
    subroutine elpa_cholesky_single_complex (self, a, success)
958
      use iso_c_binding
959
      use elpa1_auxiliary_impl
960
      use precision
961
      class(elpa_impl_t)              :: self
962
963
964
965
966
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck4)               :: a(self%local_nrows,*)
#else
      complex(kind=ck4)               :: a(self%local_nrows,self%local_ncols)
#endif
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
      integer, optional               :: success
      logical                         :: success_l
      integer(kind=c_int)             :: success_internal
      logical                         :: wantDebugIntern

      if (self%get("wantDebug",success_internal) .eq. 1) then
        if (success_internal .ne. ELPA_OK) then
          print *,"Could not querry wantDebug"
          stop
        endif

        wantDebugIntern = .true.
      else
        wantDebugIntern = .false.
      endif
982
#if WANT_SINGLE_PRECISION_COMPLEX
983
      success_l = elpa_cholesky_complex_single_impl (self%na, a, self%local_nrows, self%nblk, &
984
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
985
                                                 wantDebugIntern)
986
987
988
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
      success = ELPA_ERROR
989
990
991
992
993
994
995
996
997
998
999
#endif
      if (present(success)) then
        if (success_l) then
          success = ELPA_OK
        else
          success = ELPA_ERROR
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine
1000

1001

1002
1003
    subroutine elpa_invert_trm_double_real (self, a, success)