elpa_impl.F90 45.8 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
!> \brief Fortran module which provides the implementation of the API
51
module elpa_impl
52
  use elpa_abstract_impl
53
  use, intrinsic :: iso_c_binding
54
  implicit none
55

56
57
  private
  public :: elpa_impl_allocate
58

59
!> \brief Definition of the extended elpa_impl_t type
60
  type, extends(elpa_abstract_impl_t) :: elpa_impl_t
Andreas Marek's avatar
Andreas Marek committed
61
   private
62
   type(c_ptr)         :: index = C_NULL_PTR
63

64
   !> \brief methods available with the elpa_impl_t type
65
   contains
66
     !> \brief the puplic methods
67
     ! con-/destructor
68
69
     procedure, public :: setup => elpa_setup                   !< a setup method: implemented in elpa_setup
     procedure, public :: destroy => elpa_destroy               !< a destroy method: implemented in elpa_destroy
70

71
     ! KV store
72
73
74
75
76
77
     procedure, public :: get => elpa_get_integer               !< a get method for integer key/values: implemented in elpa_get_integer
     procedure, public :: get_double => elpa_get_double         !< a get method for double key/values: implemented in elpa_get_double
     procedure, public :: is_set => elpa_is_set                 !< a method to check whether a key/value pair has been set : implemented
                                                                !< in elpa_is_set
     procedure, public :: can_set => elpa_can_set               !< a method to check whether a key/value pair can be set : implemented
                                                                !< in elpa_can_set
78

79
80
81
82
83
84

     ! timer
     procedure, public :: get_time => elpa_get_time
     procedure, public :: print_times => elpa_print_times


85
     !> \brief the private methods
86

87
     procedure, private :: elpa_set_integer                     !< private methods to implement the setting of an integer/double key/value pair
88
     procedure, private :: elpa_set_double
89

90
     procedure, private :: elpa_solve_d               !< private methods to implement the solve step for real/complex
91
                                                                !< double/single matrices
92
93
94
     procedure, private :: elpa_solve_f
     procedure, private :: elpa_solve_dc
     procedure, private :: elpa_solve_fc
95

96
97
98
99
     procedure, private :: elpa_hermitian_multiply_d            !< private methods to implement a "hermitian" multiplication of matrices a and b
     procedure, private :: elpa_hermitian_multiply_f            !< for real valued matrices:   a**T * b
     procedure, private :: elpa_hermitian_multiply_dc            !< for complex valued matrices:   a**H * b
     procedure, private :: elpa_hermitian_multiply_fc
100

101
     procedure, private :: elpa_cholesky_d            !< private methods to implement the cholesky factorisation of
102
                                                                !< real/complex double/single matrices
103
104
105
     procedure, private :: elpa_cholesky_f
     procedure, private :: elpa_cholesky_dc
     procedure, private :: elpa_cholesky_fc
106

107
     procedure, private :: elpa_invert_trm_d          !< private methods to implement the inversion of a triangular
108
                                                                !< real/complex double/single matrix
109
110
111
     procedure, private :: elpa_invert_trm_f
     procedure, private :: elpa_invert_trm_dc
     procedure, private :: elpa_invert_trm_fc
112

113
114
     procedure, private :: elpa_solve_tridi_d         !< private methods to implement the solve step for a real valued
     procedure, private :: elpa_solve_tridi_f         !< double/single tridiagonal matrix
115

116
     procedure, private :: associate_int => elpa_associate_int  !< private method to set some pointers
117

118
  end type elpa_impl_t
119

120
  !> \brief the implementation of the private methods
121
  contains
122
123
124
125
    !> \brief function to allocate an ELPA object
    !> Parameters
    !> \param   error      integer, optional to get an error code
    !> \result  obj        class(elpa_impl_t) allocated ELPA object
126
    function elpa_impl_allocate(error) result(obj)
Andreas Marek's avatar
Andreas Marek committed
127
128
      use precision
      use elpa_utilities, only : error_unit
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
129
      use elpa_generated_fortran_interfaces
Andreas Marek's avatar
Andreas Marek committed
130

131
132
133
134
      type(elpa_impl_t), pointer   :: obj
      integer, optional            :: error

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

Andreas Marek's avatar
Andreas Marek committed
136
      ! check whether init has ever been called
137
      if ( elpa_initialized() .ne. ELPA_OK) then
138
        write(error_unit, *) "elpa_allocate(): you must call elpa_init() once before creating instances of ELPA"
139
140
        if(present(error)) then
          error = ELPA_ERROR
141
        endif
Andreas Marek's avatar
Andreas Marek committed
142
143
        return
      endif
Andreas Marek's avatar
Andreas Marek committed
144

145
      obj%index = elpa_index_instance_c()
146
147

      ! Associate some important integer pointers for convenience
148
149
150
151
152
153
154
155
      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
156
157
      endif
    end function
Andreas Marek's avatar
Andreas Marek committed
158

159
160

    !c> elpa_t elpa_allocate();
161
    function elpa_impl_allocate_c(error) result(ptr) bind(C, name="elpa_allocate")
162
163
164
165
166
167
168
169
170
171
      integer(kind=c_int) :: error
      type(c_ptr) :: ptr
      type(elpa_impl_t), pointer :: obj

      obj => elpa_impl_allocate(error)
      ptr = c_loc(obj)
    end function


    !c> void elpa_deallocate(elpa_t handle);
172
    subroutine elpa_impl_deallocate_c(handle) bind(C, name="elpa_deallocate")
173
174
175
176
177
178
179
180
181
      type(c_ptr), value :: handle
      type(elpa_impl_t), pointer :: self

      call c_f_pointer(handle, self)
      call self%destroy()
      deallocate(self)
    end subroutine


182
183
184
185
    !> \brief function to setup an ELPA object and to store the MPI communicators internally
    !> Parameters
    !> \param   self       class(elpa_impl_t), the allocated ELPA object
    !> \result  error      integer, the error code
186
    function elpa_setup(self) result(error)
187
      use elpa1_impl, only : elpa_get_communicators_impl
188
      class(elpa_impl_t), intent(inout) :: self
189
      integer :: error, error2
190
      integer :: mpi_comm_rows, mpi_comm_cols, mpierr
191

192
#ifdef WITH_MPI
193
194
195
196
      error = ELPA_ERROR
      if (self%is_set("mpi_comm_parent") == 1 .and. &
          self%is_set("process_row") == 1 .and. &
          self%is_set("process_col") == 1) then
197

198
199
200
201
202
203
        mpierr = elpa_get_communicators_impl(&
                        self%get("mpi_comm_parent"), &
                        self%get("process_row"), &
                        self%get("process_col"), &
                        mpi_comm_rows, &
                        mpi_comm_cols)
204

205
206
207
        call self%set("mpi_comm_rows", mpi_comm_rows)
        call self%set("mpi_comm_cols", mpi_comm_cols)

208
        error = ELPA_OK
209
      endif
210

211
212
      if (self%is_set("mpi_comm_rows") == 1 .and. self%is_set("mpi_comm_cols") == 1) then
        error = ELPA_OK
213
      endif
214
215
216
#else
      error = ELPA_OK
#endif
217

218
219
220
221
      if (self%get("timings") == 1) then
        call self%timer%enable()
      endif

222
    end function
223

224
225

    !c> int elpa_setup(elpa_t handle);
226
    function elpa_setup_c(handle) result(error) bind(C, name="elpa_setup")
227
228
229
230
231
232
233
234
235
      type(c_ptr), intent(in), value :: handle
      type(elpa_impl_t), pointer :: self
      integer(kind=c_int) :: error

      call c_f_pointer(handle, self)
      error = self%setup()
    end function


236
237
238
239
240
241
    !> \brief subroutine to set an integer key/value pair
    !> Parameters
    !> \param   self       class(elpa_impl_t) the allocated ELPA object
    !> \param   name       string, the key
    !> \param   value      integer, the value to be set
    !> \result  error      integer, the error code
242
    subroutine elpa_set_integer(self, name, value, error)
243
      use iso_c_binding
244
245
      use elpa_generated_fortran_interfaces
      use elpa_utilities, only : error_unit
246
      class(elpa_impl_t)              :: self
247
248
      character(*), intent(in)        :: name
      integer(kind=c_int), intent(in) :: value
249
250
      integer, optional               :: error
      integer                         :: actual_error
251

252
      actual_error = elpa_index_set_int_value_c(self%index, name // c_null_char, value, 0)
253

254
255
      if (present(error)) then
        error = actual_error
256

257
      else if (actual_error /= ELPA_OK) then
258
259
        write(error_unit,'(a,i0,a)') "ELPA: Error setting option '" // name // "' to value ", value, &
                " (got: " // elpa_strerr(actual_error) // ") and you did not check for errors!"
260
      end if
261
262
    end subroutine

263
264

    !c> void elpa_set_integer(elpa_t handle, const char *name, int value, int *error);
265
    subroutine elpa_set_integer_c(handle, name_p, value, error) bind(C, name="elpa_set_integer")
266
267
268
269
270
271
272
273
274
275
276
277
278
      type(c_ptr), intent(in), value :: handle
      type(elpa_impl_t), pointer :: self
      type(c_ptr), intent(in), value :: name_p
      character(len=elpa_strlen_c(name_p)), pointer :: name
      integer(kind=c_int), intent(in), value :: value
      integer(kind=c_int), optional, intent(in) :: error

      call c_f_pointer(handle, self)
      call c_f_pointer(name_p, name)
      call elpa_set_integer(self, name, value, error)
    end subroutine


279
280
281
282
283
284
    !> \brief function to get an integer key/value pair
    !> Parameters
    !> \param   self       class(elpa_impl_t) the allocated ELPA object
    !> \param   name       string, the key
    !> \param   error      integer, optional, to store an error code
    !> \result  value      integer, the value of the key/vaue pair
285
    function elpa_get_integer(self, name, error) result(value)
286
      use iso_c_binding
287
      use elpa_generated_fortran_interfaces
288
      use elpa_utilities, only : error_unit
289
      class(elpa_impl_t)             :: self
290
291
      character(*), intent(in)       :: name
      integer(kind=c_int)            :: value
292
      integer, intent(out), optional :: error
293
      integer                        :: actual_error
294

295
296
297
298
299
300
301
      value = elpa_index_get_int_value_c(self%index, name // c_null_char, actual_error)
      if (present(error)) then
        error = actual_error
      else if (actual_error /= ELPA_OK) then
        write(error_unit,'(a)') "ELPA: Error getting option '" // name // "'" // &
                " (got: " // elpa_strerr(actual_error) // ") and you did not check for errors!"
      end if
302
    end function
Andreas Marek's avatar
Andreas Marek committed
303

Andreas Marek's avatar
Andreas Marek committed
304
305

    !c> int elpa_get_integer(elpa_t handle, const char *name, int *error);
306
    function elpa_get_integer_c(handle, name_p, error) result(value) bind(C, name="elpa_get_integer")
Andreas Marek's avatar
Andreas Marek committed
307
308
309
310
311
312
313
314
315
316
317
318
319
      type(c_ptr), intent(in), value :: handle
      type(elpa_impl_t), pointer :: self
      type(c_ptr), intent(in), value :: name_p
      character(len=elpa_strlen_c(name_p)), pointer :: name
      integer(kind=c_int)  :: value
      integer(kind=c_int), optional, intent(inout) :: error

      call c_f_pointer(handle, self)
      call c_f_pointer(name_p, name)
      value = elpa_get_integer(self, name, error)
    end function


320
321
322
323
324
    !> \brief function to check whether a key/value pair is set
    !> Parameters
    !> \param   self       class(elpa_impl_t) the allocated ELPA object
    !> \param   name       string, the key
    !> \result  state      integer, the state of the key/value pair
325
    function elpa_is_set(self, name) result(state)
326
327
      use iso_c_binding
      use elpa_generated_fortran_interfaces
328
      class(elpa_impl_t)       :: self
329
      character(*), intent(in) :: name
330
      integer                  :: state
331

332
      state = elpa_index_value_is_set_c(self%index, name // c_null_char)
333
334
    end function

335
336
337
338
339
340
    !> \brief function to check whether a key/value pair can be set
    !> Parameters
    !> \param   self       class(elpa_impl_t) the allocated ELPA object
    !> \param   name       string, the key
    !> \param   value      integer, value
    !> \result  error      integer, error code
341
342
343
344
345
346
347
348
349
350
351
352
353
    function elpa_can_set(self, name, value) result(error)
      use iso_c_binding
      use elpa_generated_fortran_interfaces
      class(elpa_impl_t)       :: self
      character(*), intent(in) :: name
      integer(kind=c_int), intent(in) :: value
      integer                  :: error

      error = elpa_index_int_is_valid_c(self%index, name // c_null_char, value)
    end function


    function elpa_value_to_string(self, option_name, error) result(string)
354
355
356
      use elpa_generated_fortran_interfaces
      class(elpa_impl_t), intent(in) :: self
      character(kind=c_char, len=*), intent(in) :: option_name
357
358
359
360
      type(c_ptr) :: ptr
      integer, intent(out), optional :: error
      integer :: val, actual_error
      character(kind=c_char, len=elpa_index_int_value_to_strlen_c(self%index, option_name // C_NULL_CHAR)), pointer :: string
361

362
363
364
365
366
367
368
369
      nullify(string)

      val = self%get(option_name, actual_error)
      if (actual_error /= ELPA_OK) then
        if (present(error)) then
          error = actual_error
        endif
        return
370
371
      endif

372
373
374
375
      actual_error = elpa_int_value_to_string_c(option_name // C_NULL_CHAR, val, ptr)
      if (c_associated(ptr)) then
        call c_f_pointer(ptr, string)
      endif
376

377
378
379
380
      if (present(error)) then
        error = actual_error
      endif
    end function
381

382
383

    subroutine elpa_set_double(self, name, value, error)
Andreas Marek's avatar
Andreas Marek committed
384
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
385
      use elpa_generated_fortran_interfaces
386
      use elpa_utilities, only : error_unit
387
      class(elpa_impl_t)              :: self
388
      character(*), intent(in)        :: name
389
      real(kind=c_double), intent(in) :: value
390
391
      integer, optional               :: error
      integer                         :: actual_error
Andreas Marek's avatar
Andreas Marek committed
392

393
      actual_error = elpa_index_set_double_value_c(self%index, name // c_null_char, value, 0)
Andreas Marek's avatar
Andreas Marek committed
394

395
396
397
      if (present(error)) then
        error = actual_error
      else if (actual_error /= ELPA_OK) then
398
399
        write(error_unit,'(a,es12.5,a)') "ELPA: Error setting option '" // name // "' to value ", value, &
                " (got: " // elpa_strerr(actual_error) // ") and you did not check for errors!"
400
401
      end if
    end subroutine
Andreas Marek's avatar
Andreas Marek committed
402
403


404
    !c> void elpa_set_double(elpa_t handle, const char *name, double value, int *error);
405
    subroutine elpa_set_double_c(handle, name_p, value, error) bind(C, name="elpa_set_double")
406
407
408
409
410
411
412
413
414
415
416
417
418
      type(c_ptr), intent(in), value :: handle
      type(elpa_impl_t), pointer :: self
      type(c_ptr), intent(in), value :: name_p
      character(len=elpa_strlen_c(name_p)), pointer :: name
      real(kind=c_double), intent(in), value :: value
      integer(kind=c_int), optional, intent(in) :: error

      call c_f_pointer(handle, self)
      call c_f_pointer(name_p, name)
      call elpa_set_double(self, name, value, error)
    end subroutine


419
    function elpa_get_double(self, name, error) result(value)
Andreas Marek's avatar
Andreas Marek committed
420
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
421
      use elpa_generated_fortran_interfaces
422
      use elpa_utilities, only : error_unit
423
      class(elpa_impl_t)             :: self
424
      character(*), intent(in)       :: name
425
      real(kind=c_double)            :: value
426
      integer, intent(out), optional :: error
427
      integer                        :: actual_error
428

429
430
431
432
433
434
435
      value = elpa_index_get_double_value_c(self%index, name // c_null_char, actual_error)
      if (present(error)) then
        error = actual_error
      else if (actual_error /= ELPA_OK) then
        write(error_unit,'(a)') "ELPA: Error getting option '" // name // "'" // &
                " (got: " // elpa_strerr(actual_error) // ") and you did not check for errors!"
      end if
436
    end function
Andreas Marek's avatar
Andreas Marek committed
437

Andreas Marek's avatar
Andreas Marek committed
438
    !c> int elpa_get_double(elpa_t handle, const char *name, int *error);
439
    function elpa_get_double_c(handle, name_p, error) result(value) bind(C, name="elpa_get_double")
Andreas Marek's avatar
Andreas Marek committed
440
441
442
443
444
445
446
447
448
449
450
451
452
      type(c_ptr), intent(in), value :: handle
      type(elpa_impl_t), pointer :: self
      type(c_ptr), intent(in), value :: name_p
      character(len=elpa_strlen_c(name_p)), pointer :: name
      real(kind=c_double)  :: value
      integer(kind=c_int), optional, intent(inout) :: error

      call c_f_pointer(handle, self)
      call c_f_pointer(name_p, name)
      value = elpa_get_double(self, name, error)
    end function


453
    function elpa_associate_int(self, name) result(value)
Andreas Marek's avatar
Andreas Marek committed
454
      use iso_c_binding
455
      use elpa_generated_fortran_interfaces
456
457
      use elpa_utilities, only : error_unit
      class(elpa_impl_t)             :: self
458
459
      character(*), intent(in)       :: name
      integer(kind=c_int), pointer   :: value
Andreas Marek's avatar
Andreas Marek committed
460

461
462
      type(c_ptr)                    :: value_p

463
      value_p = elpa_index_get_int_loc_c(self%index, name // c_null_char)
464
465
466
      if (.not. c_associated(value_p)) then
        write(error_unit, '(a,a,a)') "ELPA: Warning, received NULL pointer for entry '", name, "'"
      endif
467
468
      call c_f_pointer(value_p, value)
    end function
Andreas Marek's avatar
Andreas Marek committed
469

470

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
    function elpa_get_time(self, name1, name2, name3, name4, name5, name6) result(s)
      class(elpa_impl_t), intent(in) :: self
      ! this is clunky, but what can you do..
      character(len=*), intent(in), optional :: name1, name2, name3, name4, name5, name6
      real(kind=c_double) :: s

      s = self%timer%get(name1, name2, name3, name4, name5, name6)
    end function


    subroutine elpa_print_times(self)
      class(elpa_impl_t), intent(in) :: self
      call self%timer%print()
    end subroutine


487
    subroutine elpa_solve_d(self, a, ev, q, error)
488
489
      use elpa2_impl
      use elpa1_impl
490
      use elpa_utilities, only : error_unit
491
      use precision
Andreas Marek's avatar
Andreas Marek committed
492
      use iso_c_binding
493
      class(elpa_impl_t)  :: self
Andreas Marek's avatar
Andreas Marek committed
494

495
496
497
#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
498
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
499
#endif
500
      real(kind=c_double) :: ev(self%na)
501

502
503
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
504
      logical             :: success_l
505

506

507
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
508
        success_l = elpa_solve_evp_real_1stage_double_impl(self, a, ev, q)
509

510
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
511
        success_l = elpa_solve_evp_real_2stage_double_impl(self, a, ev, q)
512
513
514
515
      else
        print *,"unknown solver"
        stop
      endif
516

517
      if (present(error)) then
518
        if (success_l) then
519
          error = ELPA_OK
520
        else
521
          error = ELPA_ERROR
522
523
524
525
526
527
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
    end subroutine

528
529
    !c> void elpa_solve_d(elpa_t handle, double *a, double *ev, double *q, int *error);
    subroutine elpa_solve_d_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_d")
530
531
532
533
534
535
536
537
538
539
540
      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
      integer(kind=c_int), optional, intent(in) :: error

      real(kind=c_double), pointer :: a(:, :), q(:, :), ev(:)
      type(elpa_impl_t), pointer  :: self

      call c_f_pointer(handle, self)
      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
      call c_f_pointer(ev_p, ev, [self%na])
      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])

541
      call elpa_solve_d(self, a, ev, q, error)
542
543
    end subroutine

544
    subroutine elpa_solve_f(self, a, ev, q, error)
545
546
      use elpa2_impl
      use elpa1_impl
547
      use elpa_utilities, only : error_unit
548
      use precision
549
      use iso_c_binding
550
      class(elpa_impl_t)  :: self
551
552
553
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
554
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
555
#endif
556
      real(kind=c_float)  :: ev(self%na)
557

558
559
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
560
      logical             :: success_l
561

562
#ifdef WANT_SINGLE_PRECISION_REAL
563

564
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
565
        success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q)
566

567
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
568
        success_l = elpa_solve_evp_real_2stage_single_impl(self, a, ev, q)
569
570
571
572
      else
        print *,"unknown solver"
        stop
      endif
573

574
      if (present(error)) then
575
        if (success_l) then
576
          error = ELPA_OK
577
        else
578
          error = ELPA_ERROR
579
580
581
582
583
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
584
      print *,"This installation of the ELPA library has not been build with single-precision support"
585
      error = ELPA_ERROR
586
587
588
#endif
    end subroutine

589

590
591
    !c> void elpa_solve_f(elpa_t handle, float *a, float *ev, float *q, int *error);
    subroutine elpa_solve_f_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_f")
592
593
594
595
596
597
598
599
600
601
602
      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
      integer(kind=c_int), optional, intent(in) :: error

      real(kind=c_float), pointer :: a(:, :), q(:, :), ev(:)
      type(elpa_impl_t), pointer  :: self

      call c_f_pointer(handle, self)
      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
      call c_f_pointer(ev_p, ev, [self%na])
      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])

603
      call elpa_solve_f(self, a, ev, q, error)
604
605
606
    end subroutine


607
    subroutine elpa_solve_dc(self, a, ev, q, error)
608
609
      use elpa2_impl
      use elpa1_impl
610
      use elpa_utilities, only : error_unit
611
      use precision
612
      use iso_c_binding
613
      class(elpa_impl_t)             :: self
614

615
616
617
#ifdef USE_ASSUMED_SIZE
      complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
618
      complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
619
#endif
620
      real(kind=c_double)            :: ev(self%na)
621

622
623
      integer, optional              :: error
      integer(kind=c_int)            :: error_actual
624
      logical                        :: success_l
625

626
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
627
        success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q)
628

629
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
630
        success_l = elpa_solve_evp_complex_2stage_double_impl(self,  a, ev, q)
631
632
633
634
      else
        print *,"unknown solver"
        stop
      endif
635

636
      if (present(error)) then
637
        if (success_l) then
638
          error = ELPA_OK
639
        else
640
          error = ELPA_ERROR
641
642
643
644
645
646
647
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
    end subroutine


648
649
    !c> void elpa_solve_dc(elpa_t handle, double complex *a, double *ev, double complex *q, int *error);
    subroutine elpa_solve_dc_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_dc")
650
651
652
653
654
655
656
657
658
659
660
661
      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
      integer(kind=c_int), optional, intent(in) :: error

      complex(kind=c_double_complex), pointer :: a(:, :), q(:, :)
      real(kind=c_double), pointer :: ev(:)
      type(elpa_impl_t), pointer  :: self

      call c_f_pointer(handle, self)
      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
      call c_f_pointer(ev_p, ev, [self%na])
      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])

662
      call elpa_solve_dc(self, a, ev, q, error)
663
664
665
    end subroutine


666
    subroutine elpa_solve_fc(self, a, ev, q, error)
667
668
      use elpa2_impl
      use elpa1_impl
669
670
671
      use elpa_utilities, only : error_unit

      use iso_c_binding
672
      use precision
673
      class(elpa_impl_t)            :: self
674
675
676
677
678
679
#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)
680

681
682
      integer, optional             :: error
      integer(kind=c_int)           :: error_actual
683
      logical                       :: success_l
684
685

#ifdef WANT_SINGLE_PRECISION_COMPLEX
686

687
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
688
        success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q)
689

690
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
691
        success_l = elpa_solve_evp_complex_2stage_single_impl(self,  a, ev, q)
692
693
694
695
      else
        print *,"unknown solver"
        stop
      endif
696

697
      if (present(error)) then
698
        if (success_l) then
699
          error = ELPA_OK
700
        else
701
          error = ELPA_ERROR
702
703
704
705
706
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
707
      print *,"This installation of the ELPA library has not been build with single-precision support"
708
      error = ELPA_ERROR
709
710
711
#endif
    end subroutine

712

713
714
    !c> void elpa_solve_fc(elpa_t handle, float complex *a, float *ev, float complex *q, int *error);
    subroutine elpa_solve_fc_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_fc")
715
716
717
718
719
720
721
722
723
724
725
726
      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
      integer(kind=c_int), optional, intent(in) :: error

      complex(kind=c_float_complex), pointer :: a(:, :), q(:, :)
      real(kind=c_float), pointer :: ev(:)
      type(elpa_impl_t), pointer  :: self

      call c_f_pointer(handle, self)
      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
      call c_f_pointer(ev_p, ev, [self%na])
      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])

727
      call elpa_solve_fc(self, a, ev, q, error)
728
729
730
    end subroutine


731
    subroutine elpa_hermitian_multiply_d (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
732
                                          c, ldc, ldcCols, error)
733
      use iso_c_binding
734
      use elpa1_auxiliary_impl
735
      use precision
736
      class(elpa_impl_t)              :: self
737
      character*1                     :: uplo_a, uplo_c
738
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
739
740
741
#ifdef USE_ASSUMED_SIZE
      real(kind=rk8)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
742
      real(kind=rk8)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
743
#endif
744
      integer, optional               :: error
745
746
      logical                         :: success_l

747
748
      success_l = elpa_mult_at_b_real_double_impl(self, uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                                  c, ldc, ldcCols)
749
      if (present(error)) then
750
        if (success_l) then
751
          error = ELPA_OK
752
        else
753
          error = ELPA_ERROR
754
755
        endif
      else if (.not. success_l) then
756
        write(error_unit,'(a)') "ELPA: Error in hermitian_multiply() and you did not check for errors!"
757
758
759
      endif
    end subroutine

760

761
    subroutine elpa_hermitian_multiply_f (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
762
                                          c, ldc, ldcCols, error)
763
      use iso_c_binding
764
      use elpa1_auxiliary_impl
765
      use precision
766
      class(elpa_impl_t)              :: self
767
      character*1                     :: uplo_a, uplo_c
768
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
769
770
771
#ifdef USE_ASSUMED_SIZE
      real(kind=rk4)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
772
      real(kind=rk4)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
773
#endif
774
      integer, optional               :: error
775
776
      logical                         :: success_l
#ifdef WANT_SINGLE_PRECISION_REAL
777
778
      success_l = elpa_mult_at_b_real_single_impl(self, uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                                  c, ldc, ldcCols)
779
      if (present(error)) then
780
        if (success_l) then
781
          error = ELPA_OK
782
        else
783
          error = ELPA_ERROR
784
785
        endif
      else if (.not. success_l) then
786
        write(error_unit,'(a)') "ELPA: Error in hermitian_multiply() and you did not check for errors!"
787
      endif
788
789
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
790
      error = ELPA_ERROR
791
792
793
#endif
    end subroutine

794

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

811
812
      success_l = elpa_mult_ah_b_complex_double_impl(self, uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                                     c, ldc, ldcCols)
813
      if (present(error)) then
814
        if (success_l) then
815
          error = ELPA_OK
816
        else
817
          error = ELPA_ERROR
818
819
        endif
      else if (.not. success_l) then
820
        write(error_unit,'(a)') "ELPA: Error in hermitian_multiply() and you did not check for errors!"
821
822
823
      endif
    end subroutine

824

825
    subroutine elpa_hermitian_multiply_fc (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
826
                                          c, ldc, ldcCols, error)
827
      use iso_c_binding
828
      use elpa1_auxiliary_impl
829
      use precision
830
      class(elpa_impl_t)              :: self
831
      character*1                     :: uplo_a, uplo_c
832
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
833
834
835
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck4)               :: a(lda,*), b(ldb,*), c(ldc,*)
#else
836
      complex(kind=ck4)               :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
837
#endif
838
      integer, optional               :: error
839
840
841
      logical                         :: success_l

#ifdef WANT_SINGLE_PRECISION_COMPLEX
842
843
      success_l = elpa_mult_ah_b_complex_single_impl(self, uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                                     c, ldc, ldcCols)
844
      if (present(error)) then
845
        if (success_l) then
846
          error = ELPA_OK
847
        else
848
          error = ELPA_ERROR
849
850
        endif
      else if (.not. success_l) then
851
852
        write(error_unit,'(a)') "ELPA: Error in hermitian_multiply() and you did not check for errors!"
      endif
853
854
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
855
      error = ELPA_ERROR
856
857
858
#endif
    end subroutine

859

860
    subroutine elpa_cholesky_d (self, a, error)
861
      use iso_c_binding
862
      use elpa1_auxiliary_impl