elpa_impl.F90 84 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

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

70
     ! KV store
71 72 73 74
     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
75

76 77 78 79 80 81

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


82
     !> \brief the private methods
83

84
     procedure, private :: elpa_eigenvectors_d                  !< private methods to implement the solve step for real/complex
85
                                                                !< double/single matrices
86 87 88
     procedure, private :: elpa_eigenvectors_f
     procedure, private :: elpa_eigenvectors_dc
     procedure, private :: elpa_eigenvectors_fc
89

Andreas Marek's avatar
Andreas Marek committed
90 91 92 93 94 95
     procedure, private :: elpa_eigenvalues_d                   !< private methods to implement the solve step for real/complex
                                                                !< double/single matrices; only the eigenvalues are computed
     procedure, private :: elpa_eigenvalues_f
     procedure, private :: elpa_eigenvalues_dc
     procedure, private :: elpa_eigenvalues_fc

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 161 162 163
    !c> /*! \brief C interface for the implementation of the elpa_allocate method
    !c> *
    !c> *  \param  none
    !c> *  \result elpa_t handle
    !c> */
164
    !c> elpa_t elpa_allocate();
165
    function elpa_impl_allocate_c(error) result(ptr) bind(C, name="elpa_allocate")
166 167 168 169 170 171 172 173
      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

174 175 176 177 178
    !c> /*! \brief C interface for the implementation of the elpa_deallocate method
    !c> *
    !c> *  \param  elpa_t  handle of ELPA object to be deallocated
    !c> *  \result void
    !c> */
179
    !c> void elpa_deallocate(elpa_t handle);
180
    subroutine elpa_impl_deallocate_c(handle) bind(C, name="elpa_deallocate")
181 182 183 184 185 186 187 188 189
      type(c_ptr), value :: handle
      type(elpa_impl_t), pointer :: self

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


190 191 192 193
    !> \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
194
    function elpa_setup(self) result(error)
195
      use elpa1_impl, only : elpa_get_communicators_impl
196
      class(elpa_impl_t), intent(inout) :: self
197 198 199
      integer                           :: error
      integer                           :: mpi_comm_parent, mpi_comm_rows, mpi_comm_cols, &
                                           mpierr, process_row, process_col, timings
200

201
#ifdef WITH_MPI
202 203 204 205
      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
206

207 208 209
        call self%get("mpi_comm_parent", mpi_comm_parent)
        call self%get("process_row", process_row)
        call self%get("process_col", process_col)
210
        mpierr = elpa_get_communicators_impl(&
211 212 213
                        mpi_comm_parent, &
                        process_row, &
                        process_col, &
214 215
                        mpi_comm_rows, &
                        mpi_comm_cols)
216

217 218 219
        call self%set("mpi_comm_rows", mpi_comm_rows)
        call self%set("mpi_comm_cols", mpi_comm_cols)

220
        error = ELPA_OK
221
      endif
222

223 224
      if (self%is_set("mpi_comm_rows") == 1 .and. self%is_set("mpi_comm_cols") == 1) then
        error = ELPA_OK
225
      endif
226 227 228
#else
      error = ELPA_OK
#endif
229

230
#ifdef HAVE_DETAILED_TIMINGS
231 232
      call self%get("timings",timings)
      if (timings == 1) then
233 234
        call self%timer%enable()
      endif
235
#endif
236

237
    end function
238

239 240 241 242 243 244
    !c> /*! \brief C interface for the implementation of the elpa_setup method
    !c> *
    !c> *  \param  elpa_t  handle of the ELPA object which describes the problem to
    !c> *                  be set up
    !c> *  \result int     error code, which can be queried with elpa_strerr
    !c> */
245
    !c> int elpa_setup(elpa_t handle);
246
    function elpa_setup_c(handle) result(error) bind(C, name="elpa_setup")
247 248 249 250 251 252 253 254 255
      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


256 257 258 259 260 261 262 263 264
    !c> /*! \brief C interface for the implementation of the elpa_set_integer method
    !c> *  This method is available to the user as C generic elpa_set method
    !c> *
    !c> *  \param  handle  handle of the ELPA object for which a key/value pair should be set
    !c> *  \param  name    the name of the key
    !c> *  \param  value   the value to be set for the key
    !c> *  \param  error   on return the error code, which can be queried with elpa_strerr()
    !c> *  \result void
    !c> */
265
    !c> void elpa_set_integer(elpa_t handle, const char *name, int value, int *error);
266
    subroutine elpa_set_integer_c(handle, name_p, value, error) bind(C, name="elpa_set_integer")
267 268 269 270 271 272 273 274 275 276 277 278 279
      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


280 281 282 283 284 285 286 287 288
    !c> /*! \brief C interface for the implementation of the elpa_get_integer method
    !c> *  This method is available to the user as C generic elpa_get method
    !c> *
    !c> *  \param  handle  handle of the ELPA object for which a key/value pair should be queried
    !c> *  \param  name    the name of the key
    !c> *  \param  value   the value to be obtain for the key
    !c> *  \param  error   on return the error code, which can be queried with elpa_strerr()
    !c> *  \result void
    !c> */
289 290
    !c> void elpa_get_integer(elpa_t handle, const char *name, int *value, int *error);
    subroutine elpa_get_integer_c(handle, name_p, value, error) bind(C, name="elpa_get_integer")
Andreas Marek's avatar
Andreas Marek committed
291 292 293 294 295 296 297 298 299
      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)
300 301
      call elpa_get_integer(self, name, value, error)
    end subroutine
Andreas Marek's avatar
Andreas Marek committed
302 303


304 305 306 307 308
    !> \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
309
    function elpa_is_set(self, name) result(state)
310 311
      use iso_c_binding
      use elpa_generated_fortran_interfaces
312
      class(elpa_impl_t)       :: self
313
      character(*), intent(in) :: name
314
      integer                  :: state
315

316
      state = elpa_index_value_is_set_c(self%index, name // c_null_char)
317 318
    end function

319 320 321 322 323 324
    !> \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
325 326 327 328 329 330 331 332 333 334 335 336 337
    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)
338 339 340
      use elpa_generated_fortran_interfaces
      class(elpa_impl_t), intent(in) :: self
      character(kind=c_char, len=*), intent(in) :: option_name
341 342 343 344
      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
345

346 347
      nullify(string)

348
      call self%get(option_name, val, actual_error)
349 350 351 352 353
      if (actual_error /= ELPA_OK) then
        if (present(error)) then
          error = actual_error
        endif
        return
354 355
      endif

356 357 358 359
      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
360

361 362 363 364
      if (present(error)) then
        error = actual_error
      endif
    end function
365

Andreas Marek's avatar
Andreas Marek committed
366

367 368 369 370 371 372 373 374 375
    !c> /*! \brief C interface for the implementation of the elpa_set_double method
    !c> *  This method is available to the user as C generic elpa_set method
    !c> *
    !c> *  \param  handle  handle of the ELPA object for which a key/value pair should be set
    !c> *  \param  name    the name of the key
    !c> *  \param  value   the value to be set for the key
    !c> *  \param  error   on return the error code, which can be queried with elpa_strerr()
    !c> *  \result void
    !c> */
376
    !c> void elpa_set_double(elpa_t handle, const char *name, double value, int *error);
377
    subroutine elpa_set_double_c(handle, name_p, value, error) bind(C, name="elpa_set_double")
378 379 380 381 382 383 384 385 386 387 388 389
      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

390

391
    !c> /*! \brief C interface for the implementation of the elpa_get_double method
392 393 394 395 396 397 398 399
    !c> *  This method is available to the user as C generic elpa_get method
    !c> *
    !c> *  \param  handle  handle of the ELPA object for which a key/value pair should be queried
    !c> *  \param  name    the name of the key
    !c> *  \param  value   the value to be obtain for the key
    !c> *  \param  error   on return the error code, which can be queried with elpa_strerr()
    !c> *  \result void
    !c> */
400 401
    !c> void elpa_get_double(elpa_t handle, const char *name, double *value, int *error);
    subroutine elpa_get_double_c(handle, name_p, value, error) bind(C, name="elpa_get_double")
Andreas Marek's avatar
Andreas Marek committed
402 403 404 405 406 407 408 409 410
      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)
411 412
      call elpa_get_double(self, name, value, error)
    end subroutine
Andreas Marek's avatar
Andreas Marek committed
413 414


415
    function elpa_associate_int(self, name) result(value)
Andreas Marek's avatar
Andreas Marek committed
416
      use iso_c_binding
417
      use elpa_generated_fortran_interfaces
418 419
      use elpa_utilities, only : error_unit
      class(elpa_impl_t)             :: self
420 421
      character(*), intent(in)       :: name
      integer(kind=c_int), pointer   :: value
Andreas Marek's avatar
Andreas Marek committed
422

423 424
      type(c_ptr)                    :: value_p

425
      value_p = elpa_index_get_int_loc_c(self%index, name // c_null_char)
426 427 428
      if (.not. c_associated(value_p)) then
        write(error_unit, '(a,a,a)') "ELPA: Warning, received NULL pointer for entry '", name, "'"
      endif
429 430
      call c_f_pointer(value_p, value)
    end function
Andreas Marek's avatar
Andreas Marek committed
431

432

433 434 435 436 437 438
    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

439
#ifdef HAVE_DETAILED_TIMINGS
440
      s = self%timer%get(name1, name2, name3, name4, name5, name6)
441 442 443
#else
      s = -1.0
#endif
444 445 446 447 448
    end function


    subroutine elpa_print_times(self)
      class(elpa_impl_t), intent(in) :: self
449
#ifdef HAVE_DETAILED_TIMINGS
450
      call self%timer%print()
451
#endif
452 453
    end subroutine

454
    !>  \brief elpa_eigenvectors_d: class method to solve the eigenvalue problem for double real matrices
Andreas Marek's avatar
Andreas Marek committed
455
    !>
456 457
    !>  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
458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478
    !>  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
479
    subroutine elpa_eigenvectors_d(self, a, ev, q, error)
480 481
      use elpa2_impl
      use elpa1_impl
482
      use elpa_utilities, only : error_unit
Andreas Marek's avatar
Andreas Marek committed
483
      use iso_c_binding
484
      class(elpa_impl_t)  :: self
Andreas Marek's avatar
Andreas Marek committed
485

486 487 488
#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
489
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
490
#endif
491
      real(kind=c_double) :: ev(self%na)
492

493
      integer, optional   :: error
494
      integer(kind=c_int) :: solver
495
      logical             :: success_l
496

497

498 499
      call self%get("solver", solver)
      if (solver .eq. ELPA_SOLVER_1STAGE) then
500
        success_l = elpa_solve_evp_real_1stage_double_impl(self, a, ev, q)
501

502
      else if (solver .eq. ELPA_SOLVER_2STAGE) then
503
        success_l = elpa_solve_evp_real_2stage_double_impl(self, a, ev, q)
504 505 506 507
      else
        print *,"unknown solver"
        stop
      endif
508

509
      if (present(error)) then
510
        if (success_l) then
511
          error = ELPA_OK
512
        else
513
          error = ELPA_ERROR
514 515 516 517 518 519
        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

520 521
    !c> void elpa_eigenvectors_d(elpa_t handle, double *a, double *ev, double *q, int *error);
    subroutine elpa_eigenvectors_d_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_eigenvectors_d")
522 523 524 525 526 527 528 529 530 531 532
      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])

533
      call elpa_eigenvectors_d(self, a, ev, q, error)
534 535
    end subroutine

Andreas Marek's avatar
Andreas Marek committed
536

537
    !>  \brief elpa_eigenvectors_f: class method to solve the eigenvalue problem for float real matrices
Andreas Marek's avatar
Andreas Marek committed
538
    !>
539 540
    !>  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
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
    !>  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
562
    subroutine elpa_eigenvectors_f(self, a, ev, q, error)
563 564
      use elpa2_impl
      use elpa1_impl
565 566
      use elpa_utilities, only : error_unit
      use iso_c_binding
567
      class(elpa_impl_t)  :: self
568 569 570
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
571
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
572
#endif
573
      real(kind=c_float)  :: ev(self%na)
574

575
      integer, optional   :: error
576
      integer(kind=c_int) :: solver
577
#ifdef WANT_SINGLE_PRECISION_REAL
578
      logical             :: success_l
579

580 581
      call self%get("solver",solver)
      if (solver .eq. ELPA_SOLVER_1STAGE) then
582
        success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q)
583

584
      else if (solver .eq. ELPA_SOLVER_2STAGE) then
585
        success_l = elpa_solve_evp_real_2stage_single_impl(self, a, ev, q)
586 587 588 589
      else
        print *,"unknown solver"
        stop
      endif
590

591
      if (present(error)) then
592
        if (success_l) then
593
          error = ELPA_OK
594
        else
595
          error = ELPA_ERROR
596 597 598 599 600
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
601
      print *,"This installation of the ELPA library has not been build with single-precision support"
602
      error = ELPA_ERROR
603 604 605
#endif
    end subroutine

606

607 608
    !c> void elpa_eigenvectors_f(elpa_t handle, float *a, float *ev, float *q, int *error);
    subroutine elpa_eigenvectors_f_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_eigenvectors_f")
609 610 611 612 613 614 615 616 617 618 619
      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])

620
      call elpa_eigenvectors_f(self, a, ev, q, error)
621 622 623
    end subroutine


624
    !>  \brief elpa_eigenvectors_dc: class method to solve the eigenvalue problem for double complex matrices
Andreas Marek's avatar
Andreas Marek committed
625
    !>
626 627
    !>  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
628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648
    !>  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
649
    subroutine elpa_eigenvectors_dc(self, a, ev, q, error)
650 651
      use elpa2_impl
      use elpa1_impl
652 653
      use elpa_utilities, only : error_unit
      use iso_c_binding
654
      class(elpa_impl_t)             :: self
655

656 657 658
#ifdef USE_ASSUMED_SIZE
      complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
659
      complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
660
#endif
661
      real(kind=c_double)            :: ev(self%na)
662

663
      integer, optional              :: error
664
      integer(kind=c_int)            :: solver
665
      logical                        :: success_l
666

667 668
      call self%get("solver", solver)
      if (solver .eq. ELPA_SOLVER_1STAGE) then
669
        success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q)
670

671
      else if (solver .eq. ELPA_SOLVER_2STAGE) then
672
        success_l = elpa_solve_evp_complex_2stage_double_impl(self,  a, ev, q)
673 674 675 676
      else
        print *,"unknown solver"
        stop
      endif
677

678
      if (present(error)) then
679
        if (success_l) then
680
          error = ELPA_OK
681
        else
682
          error = ELPA_ERROR
683 684 685 686 687 688 689
        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


690 691
    !c> void elpa_eigenvectors_dc(elpa_t handle, double complex *a, double *ev, double complex *q, int *error);
    subroutine elpa_eigenvectors_dc_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_eigenvectors_dc")
692 693 694 695 696 697 698 699 700 701 702 703
      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])

704
      call elpa_eigenvectors_dc(self, a, ev, q, error)
705 706 707
    end subroutine


708
    !>  \brief elpa_eigenvectors_fc: class method to solve the eigenvalue problem for float complex matrices
Andreas Marek's avatar
Andreas Marek committed
709
    !>
710 711
    !>  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
712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732
    !>  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
733
    subroutine elpa_eigenvectors_fc(self, a, ev, q, error)
734 735
      use elpa2_impl
      use elpa1_impl
736 737 738
      use elpa_utilities, only : error_unit

      use iso_c_binding
739
      class(elpa_impl_t)            :: self
740
#ifdef USE_ASSUMED_SIZE
Andreas Marek's avatar
Andreas Marek committed
741
      complex(kind=c_float_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
742
#else
Andreas Marek's avatar
Andreas Marek committed
743
      complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
744
#endif
Andreas Marek's avatar
Andreas Marek committed
745
      real(kind=c_float)            :: ev(self%na)
746

747
      integer, optional             :: error
748
      integer(kind=c_int)           :: solver
749
#ifdef WANT_SINGLE_PRECISION_COMPLEX
750
      logical                       :: success_l
751

752 753
      call self%get("solver", solver)
      if (solver .eq. ELPA_SOLVER_1STAGE) then
754
        success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q)
755

756
      else if (solver .eq. ELPA_SOLVER_2STAGE) then
757
        success_l = elpa_solve_evp_complex_2stage_single_impl(self,  a, ev, q)
758 759 760 761
      else
        print *,"unknown solver"
        stop
      endif
762

763
      if (present(error)) then
764
        if (success_l) then
765
          error = ELPA_OK
766
        else
767
          error = ELPA_ERROR
768 769 770 771 772
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
#else
773
      print *,"This installation of the ELPA library has not been build with single-precision support"
774
      error = ELPA_ERROR
775 776 777
#endif
    end subroutine

778

779 780
    !c> void elpa_eigenvectors_fc(elpa_t handle, float complex *a, float *ev, float complex *q, int *error);
    subroutine elpa_eigenvectors_fc_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_eigenvectors_fc")
781 782 783 784 785 786 787 788 789 790 791 792
      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])

793
      call elpa_eigenvectors_fc(self, a, ev, q, error)
794 795
    end subroutine

Andreas Marek's avatar
Andreas Marek committed
796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 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 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911



    !>  \brief elpa_eigenvalues_d: class method to solve the eigenvalue problem for double real matrices
    !>
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
    !>  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 error                                integer, optional: returns an error code, which can be queried with elpa_strerr
    subroutine elpa_eigenvalues_d(self, a, ev, error)
      use elpa2_impl
      use elpa1_impl
      use elpa_utilities, only : error_unit
      use iso_c_binding
      class(elpa_impl_t)  :: self

#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *)
#else
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols)
#endif
      real(kind=c_double) :: ev(self%na)

      integer, optional   :: error
      integer(kind=c_int) :: solver
      logical             :: success_l


      call self%get("solver", solver)
      if (solver .eq. ELPA_SOLVER_1STAGE) then
        success_l = elpa_solve_evp_real_1stage_double_impl(self, a, ev)

      else if (solver .eq. ELPA_SOLVER_2STAGE) then
        success_l = elpa_solve_evp_real_2stage_double_impl(self, a, ev)
      else
        print *,"unknown solver"
        stop
      endif

      if (present(error)) then
        if (success_l) then
          error = ELPA_OK
        else
          error = 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
    end subroutine

    !c> void elpa_eigenvalues_d(elpa_t handle, double *a, double *ev, int *error);
    subroutine elpa_eigenvalues_d_c(handle, a_p, ev_p, error) bind(C, name="elpa_eigenvalues_d")
      type(c_ptr), intent(in), value :: handle, a_p, ev_p
      integer(kind=c_int), optional, intent(in) :: error

      real(kind=c_double), pointer :: a(:, :), 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 elpa_eigenvalues_d(self, a, ev, error)
    end subroutine


    !>  \brief elpa_eigenvectors_f: class method to solve the eigenvalue problem for float real matrices
    !>
    !>  The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
    !>  blocksize, the number of eigenvectors
    !>  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 error                                integer, optional: returns an error code, which can be queried with elpa_strerr
    subroutine elpa_eigenvalues_f(self, a, ev, error)
      use elpa2_impl
      use elpa1_impl
      use elpa_utilities, only : error_unit
      use iso_c_binding
      class(elpa_impl_t)  :: self
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *)
#else
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols)
#endif
      real(kind=c_float)  :: ev(self%na)

      integer, optional   :: error
      integer(kind=c_int) :: solver
#ifdef WANT_SINGLE_PRECISION_REAL
912
      logical             :: success_l
Andreas Marek's avatar
Andreas Marek committed
913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 <