elpa_impl.F90 74 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 actual implementation of the API. Do not use directly! Use the module "elpa"
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

Andreas Marek's avatar
Andreas Marek committed
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
     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
Andreas Marek's avatar
Andreas Marek committed
98
     procedure, private :: elpa_hermitian_multiply_dc           !< for complex valued matrices:   a**H * b
99
     procedure, private :: elpa_hermitian_multiply_fc
100

Andreas Marek's avatar
Andreas Marek committed
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

Andreas Marek's avatar
Andreas Marek committed
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

Andreas Marek's avatar
Andreas Marek committed
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
    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

Andreas Marek's avatar
Andreas Marek committed
486 487
    !>  \brief elpa_solve_d: class method to solve the eigenvalue problem for double real matrices
    !>
488 489
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
Andreas Marek's avatar
Andreas Marek committed
490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
    !>  to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
    !>  with the class method "setup"
    !>
    !>  It is possible to change the behaviour of the method by setting tunable parameters with the
    !>  class method "set"
    !>
    !>  Parameters
    !>
    !>  \param a                                    Distributed matrix for which eigenvalues are to be computed.
    !>                                              Distribution is like in Scalapack.
    !>                                              The full matrix must be set (not only one half like in scalapack).
    !>                                              Destroyed on exit (upper and lower half).
    !>
    !>  \param ev                                   On output: eigenvalues of a, every processor gets the complete set
    !>
    !>  \param q                                    On output: Eigenvectors of a
    !>                                              Distribution is like in Scalapack.
    !>                                              Must be always dimensioned to the full size (corresponding to (na,na))
    !>                                              even if only a part of the eigenvalues is needed.
    !>
    !>  \param error                                integer, optional: returns an error code, which can be queried with elpa_strerr
511
    subroutine elpa_solve_d(self, a, ev, q, error)
512 513
      use elpa2_impl
      use elpa1_impl
514
      use elpa_utilities, only : error_unit
Andreas Marek's avatar
Andreas Marek committed
515
      use iso_c_binding
516
      class(elpa_impl_t)  :: self
Andreas Marek's avatar
Andreas Marek committed
517

518 519 520
#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
521
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
522
#endif
523
      real(kind=c_double) :: ev(self%na)
524

525 526
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
527
      logical             :: success_l
528

529

530
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
531
        success_l = elpa_solve_evp_real_1stage_double_impl(self, a, ev, q)
532

533
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
534
        success_l = elpa_solve_evp_real_2stage_double_impl(self, a, ev, q)
535 536 537 538
      else
        print *,"unknown solver"
        stop
      endif
539

540
      if (present(error)) then
541
        if (success_l) then
542
          error = ELPA_OK
543
        else
544
          error = ELPA_ERROR
545 546 547 548 549 550
        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

551 552
    !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")
553 554 555 556 557 558 559 560 561 562 563
      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])

564
      call elpa_solve_d(self, a, ev, q, error)
565 566
    end subroutine

Andreas Marek's avatar
Andreas Marek committed
567 568 569

    !>  \brief elpa_solve_f: class method to solve the eigenvalue problem for float real matrices
    !>
570 571
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
Andreas Marek's avatar
Andreas Marek committed
572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592
    !>  to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
    !>  with the class method "setup"
    !>
    !>  It is possible to change the behaviour of the method by setting tunable parameters with the
    !>  class method "set"
    !>
    !>  Parameters
    !>
    !>  \param a                                    Distributed matrix for which eigenvalues are to be computed.
    !>                                              Distribution is like in Scalapack.
    !>                                              The full matrix must be set (not only one half like in scalapack).
    !>                                              Destroyed on exit (upper and lower half).
    !>
    !>  \param ev                                   On output: eigenvalues of a, every processor gets the complete set
    !>
    !>  \param q                                    On output: Eigenvectors of a
    !>                                              Distribution is like in Scalapack.
    !>                                              Must be always dimensioned to the full size (corresponding to (na,na))
    !>                                              even if only a part of the eigenvalues is needed.
    !>
    !>  \param error                                integer, optional: returns an error code, which can be queried with elpa_strerr
593
    subroutine elpa_solve_f(self, a, ev, q, error)
594 595
      use elpa2_impl
      use elpa1_impl
596 597
      use elpa_utilities, only : error_unit
      use iso_c_binding
598
      class(elpa_impl_t)  :: self
599 600 601
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
602
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
603
#endif
604
      real(kind=c_float)  :: ev(self%na)
605

606 607
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
608
      logical             :: success_l
609

610
#ifdef WANT_SINGLE_PRECISION_REAL
611

612
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
613
        success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q)
614

615
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
616
        success_l = elpa_solve_evp_real_2stage_single_impl(self, a, ev, q)
617 618 619 620
      else
        print *,"unknown solver"
        stop
      endif
621

622
      if (present(error)) then
623
        if (success_l) then
624
          error = ELPA_OK
625
        else
626
          error = ELPA_ERROR
627 628 629 630 631
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
632
      print *,"This installation of the ELPA library has not been build with single-precision support"
633
      error = ELPA_ERROR
634 635 636
#endif
    end subroutine

637

638 639
    !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")
640 641 642 643 644 645 646 647 648 649 650
      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])

651
      call elpa_solve_f(self, a, ev, q, error)
652 653 654
    end subroutine


Andreas Marek's avatar
Andreas Marek committed
655 656
    !>  \brief elpa_solve_dc: class method to solve the eigenvalue problem for double complex matrices
    !>
657 658
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
Andreas Marek's avatar
Andreas Marek committed
659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
    !>  to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
    !>  with the class method "setup"
    !>
    !>  It is possible to change the behaviour of the method by setting tunable parameters with the
    !>  class method "set"
    !>
    !>  Parameters
    !>
    !>  \param a                                    Distributed matrix for which eigenvalues are to be computed.
    !>                                              Distribution is like in Scalapack.
    !>                                              The full matrix must be set (not only one half like in scalapack).
    !>                                              Destroyed on exit (upper and lower half).
    !>
    !>  \param ev                                   On output: eigenvalues of a, every processor gets the complete set
    !>
    !>  \param q                                    On output: Eigenvectors of a
    !>                                              Distribution is like in Scalapack.
    !>                                              Must be always dimensioned to the full size (corresponding to (na,na))
    !>                                              even if only a part of the eigenvalues is needed.
    !>
    !>  \param error                                integer, optional: returns an error code, which can be queried with elpa_strerr
680
    subroutine elpa_solve_dc(self, a, ev, q, error)
681 682
      use elpa2_impl
      use elpa1_impl
683 684
      use elpa_utilities, only : error_unit
      use iso_c_binding
685
      class(elpa_impl_t)             :: self
686

687 688 689
#ifdef USE_ASSUMED_SIZE
      complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
690
      complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
691
#endif
692
      real(kind=c_double)            :: ev(self%na)
693

694 695
      integer, optional              :: error
      integer(kind=c_int)            :: error_actual
696
      logical                        :: success_l
697

698
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
699
        success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q)
700

701
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
702
        success_l = elpa_solve_evp_complex_2stage_double_impl(self,  a, ev, q)
703 704 705 706
      else
        print *,"unknown solver"
        stop
      endif
707

708
      if (present(error)) then
709
        if (success_l) then
710
          error = ELPA_OK
711
        else
712
          error = ELPA_ERROR
713 714 715 716 717 718 719
        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


720 721
    !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")
722 723 724 725 726 727 728 729 730 731 732 733
      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])

734
      call elpa_solve_dc(self, a, ev, q, error)
735 736 737
    end subroutine


Andreas Marek's avatar
Andreas Marek committed
738 739
    !>  \brief elpa_solve_fc: class method to solve the eigenvalue problem for float complex matrices
    !>
740 741
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
Andreas Marek's avatar
Andreas Marek committed
742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
    !>  to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
    !>  with the class method "setup"
    !>
    !>  It is possible to change the behaviour of the method by setting tunable parameters with the
    !>  class method "set"
    !>
    !>  Parameters
    !>
    !>  \param a                                    Distributed matrix for which eigenvalues are to be computed.
    !>                                              Distribution is like in Scalapack.
    !>                                              The full matrix must be set (not only one half like in scalapack).
    !>                                              Destroyed on exit (upper and lower half).
    !>
    !>  \param ev                                   On output: eigenvalues of a, every processor gets the complete set
    !>
    !>  \param q                                    On output: Eigenvectors of a
    !>                                              Distribution is like in Scalapack.
    !>                                              Must be always dimensioned to the full size (corresponding to (na,na))
    !>                                              even if only a part of the eigenvalues is needed.
    !>
    !>  \param error                                integer, optional: returns an error code, which can be queried with elpa_strerr
763
    subroutine elpa_solve_fc(self, a, ev, q, error)
764 765
      use elpa2_impl
      use elpa1_impl
766 767 768
      use elpa_utilities, only : error_unit

      use iso_c_binding
769
      class(elpa_impl_t)            :: self
770
#ifdef USE_ASSUMED_SIZE
Andreas Marek's avatar
Andreas Marek committed
771
      complex(kind=c_float_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
772
#else
Andreas Marek's avatar
Andreas Marek committed
773
      complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
774
#endif
Andreas Marek's avatar
Andreas Marek committed
775
      real(kind=c_float)            :: ev(self%na)
776

777 778
      integer, optional             :: error
      integer(kind=c_int)           :: error_actual
779
      logical                       :: success_l
780 781

#ifdef WANT_SINGLE_PRECISION_COMPLEX
782

783
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
784
        success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q)
785

786
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
787
        success_l = elpa_solve_evp_complex_2stage_single_impl(self,  a, ev, q)
788 789 790 791
      else
        print *,"unknown solver"
        stop
      endif
792

793
      if (present(error)) then
794
        if (success_l) then
795
          error = ELPA_OK
796
        else
797
          error = ELPA_ERROR
798 799 800 801 802
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
803
      print *,"This installation of the ELPA library has not been build with single-precision support"
804
      error = ELPA_ERROR
805 806 807
#endif
    end subroutine

808

809 810
    !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")
811 812 813 814 815 816 817 818 819 820 821 822
      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])

823
      call elpa_solve_fc(self, a, ev, q, error)
824 825
    end subroutine

Andreas Marek's avatar
Andreas Marek committed
826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862
    !> \brief  elpa_hermitian_multiply_d: class method to perform C : = A**T * B for double real matrices
    !>         where   A is a square matrix (na,na) which is optionally upper or lower triangular
    !>                 B is a (na,ncb) matrix
    !>                 C is a (na,ncb) matrix where optionally only the upper or lower
    !>                   triangle may be computed
    !>
    !> the MPI commicators and the block-cyclic distribution block size are already known to the type.
    !> Thus the class method "setup" must be called BEFORE this method is used
    !>
    !> \details
    !>
    !> \param  uplo_a               'U' if A is upper triangular
    !>                              'L' if A is lower triangular
    !>                              anything else if A is a full matrix
    !>                              Please note: This pertains to the original A (as set in the calling program)
    !>                                           whereas the transpose of A is used for calculations
    !>                              If uplo_a is 'U' or 'L', the other triangle is not used at all,
    !>                              i.e. it may contain arbitrary numbers
    !> \param uplo_c                'U' if only the upper diagonal part of C is needed
    !>                              'L' if only the upper diagonal part of C is needed
    !>                              anything else if the full matrix C is needed
    !>                              Please note: Even when uplo_c is 'U' or 'L', the other triangle may be
    !>                                            written to a certain extent, i.e. one shouldn't rely on the content there!
    !> \param na                    Number of rows/columns of global matrix A, number of rows of global matrices B and C
    !> \param ncb                   Number of columns  of global matrices B and C
    !> \param a                     matrix a
    !> \param nrows_a               number of rows of local (sub) matrix a
    !> \param ncols_a               number of columns of local (sub) matrix a
    !> \param b                     matrix b
    !> \param nrows_b               number of rows of local (sub) matrix b
    !> \param ncols_b               number of columns of local (sub) matrix b
    !> \param c                     matrix c
    !> \param nrows_c               number of rows of local (sub) matrix c
    !> \param ncols_c               number of columns of local (sub) matrix c
    !> \param error                 optional argument, error code which can be queried with elpa_strerr
    subroutine elpa_hermitian_multiply_d (self,uplo_a, uplo_c, na, ncb, a, nrows_a, ncols_a, b, nrows_b, ncols_b, &
                                          c, nrows_c, ncols_c, error)
863
      use iso_c_binding
864
      use elpa1_auxiliary_impl
865
      class(elpa_impl_t)              :: self
866
      character*1                     :: uplo_a, uplo_c