elpa_impl.F90 38 KB
Newer Older
1
2
3
!
!    Copyright 2017, L. Hüdepohl and A. Marek, MPCDF
!
Andreas Marek's avatar
Andreas Marek committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), formerly known as
!      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!    This particular source code file contains additions, changes and
!    enhancements authored by Intel Corporation which is not part of
!    the ELPA consortium.
!
!    More information can be found here:
!    http://elpa.mpcdf.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
48
#include "config-f90.h"
49

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

55
56
  private
  public :: elpa_impl_allocate
57

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

62
   contains
63
     ! con-/destructor
64
     procedure, public :: setup => elpa_setup
65
     procedure, public :: destroy => elpa_destroy
66

67
     ! KV store
68
     procedure, public :: get => elpa_get_integer
69
70
     procedure, public :: get_double => elpa_get_double
     procedure, public :: is_set => elpa_is_set
71
     procedure, public :: can_set => elpa_can_set
72
73

     ! privates:
74

75
     procedure, private :: elpa_set_integer
76
     procedure, private :: elpa_set_double
77

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

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

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

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

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

     procedure, private :: associate_int => elpa_associate_int

103
  end type elpa_impl_t
104
105
106

  contains

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

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

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

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

126
      obj%index = elpa_index_instance_c()
127
128

      ! Associate some important integer pointers for convenience
129
130
131
132
133
134
135
136
      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
137
138
      endif
    end function
Andreas Marek's avatar
Andreas Marek committed
139

140

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

147
      error = ELPA_ERROR
148

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

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

160
161
162
        call self%set("mpi_comm_rows", mpi_comm_rows)
        call self%set("mpi_comm_cols", mpi_comm_cols)

163
        error = ELPA_OK
164
      endif
165

166
167
      if (self%is_set("mpi_comm_rows") == 1 .and. self%is_set("mpi_comm_cols") == 1) then
        error = ELPA_OK
168
      endif
169
170

    end function
171

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

182
      actual_error = elpa_index_set_int_value_c(self%index, name // c_null_char, value)
183

184
185
      if (present(error)) then
        error = actual_error
186

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

      end if
192
193
    end subroutine

194

195
    function elpa_get_integer(self, name, error) result(value)
196
      use iso_c_binding
197
      use elpa_generated_fortran_interfaces
198
      class(elpa_impl_t)             :: self
199
200
      character(*), intent(in)       :: name
      integer(kind=c_int)            :: value
201
      integer, intent(out), optional :: error
202

203
      value = elpa_index_get_int_value_c(self%index, name // c_null_char, error)
204

205
    end function
Andreas Marek's avatar
Andreas Marek committed
206

207
208

    function elpa_is_set(self, name) result(state)
209
210
      use iso_c_binding
      use elpa_generated_fortran_interfaces
211
      class(elpa_impl_t)       :: self
212
      character(*), intent(in) :: name
213
      integer                  :: state
214

215
      state = elpa_index_value_is_set_c(self%index, name // c_null_char)
216
217
    end function

218

219
220
221
222
223
224
225
226
227
228
229
230
231
    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)
232
233
234
      use elpa_generated_fortran_interfaces
      class(elpa_impl_t), intent(in) :: self
      character(kind=c_char, len=*), intent(in) :: option_name
235
236
237
238
      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
239

240
241
242
243
244
245
246
247
      nullify(string)

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

250
251
252
253
      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
254

255
256
257
258
      if (present(error)) then
        error = actual_error
      endif
    end function
259

260
261

    subroutine elpa_set_double(self, name, value, error)
Andreas Marek's avatar
Andreas Marek committed
262
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
263
      use elpa_generated_fortran_interfaces
264
      use elpa_utilities, only : error_unit
265
      class(elpa_impl_t)              :: self
266
      character(*), intent(in)        :: name
267
      real(kind=c_double), intent(in) :: value
268
269
      integer, optional               :: error
      integer                         :: actual_error
Andreas Marek's avatar
Andreas Marek committed
270

271
      actual_error = elpa_index_set_double_value_c(self%index, name // c_null_char, value)
Andreas Marek's avatar
Andreas Marek committed
272

273
274
      if (present(error)) then
        error = actual_error
275

276
      else if (actual_error /= ELPA_OK) then
277
        write(error_unit,'(a,a,es12.5,a)') "ELPA: Error setting option '", name, "' to value ", value, &
278
                " and you did not check for errors!"
279
280
281

      end if
    end subroutine
Andreas Marek's avatar
Andreas Marek committed
282
283


284
    function elpa_get_double(self, name, error) result(value)
Andreas Marek's avatar
Andreas Marek committed
285
      use iso_c_binding
Lorenz Huedepohl's avatar
Lorenz Huedepohl committed
286
      use elpa_generated_fortran_interfaces
287
      class(elpa_impl_t)             :: self
288
      character(*), intent(in)       :: name
289
      real(kind=c_double)            :: value
290
      integer, intent(out), optional :: error
291

292
      value = elpa_index_get_double_value_c(self%index, name // c_null_char, error)
293
294

    end function
Andreas Marek's avatar
Andreas Marek committed
295
296


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

305
306
      type(c_ptr)                    :: value_p

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

314

315
    subroutine elpa_solve_real_double(self, a, ev, q, error)
316
317
      use elpa2_impl
      use elpa1_impl
318
      use elpa_utilities, only : error_unit
319
      use precision
Andreas Marek's avatar
Andreas Marek committed
320
      use iso_c_binding
321
      class(elpa_impl_t)  :: self
Andreas Marek's avatar
Andreas Marek committed
322

323
324
325
#ifdef USE_ASSUMED_SIZE
      real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
326
      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
327
#endif
328
      real(kind=c_double) :: ev(self%na)
329
330

      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
331
332
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
333
      logical             :: success_l
334

335

336
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
337
338
        success_l = elpa_solve_evp_real_1stage_double_impl(self, a, ev, q, time_evp_fwd,     &
                                                           time_evp_solve, time_evp_back)
339

340
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
341
342
343
344
        success_l = elpa_solve_evp_real_2stage_double_impl(self, a, ev, q,  &
                                                           time_evp_fwd,     &
                                                           time_evp_solve, time_evp_back  &
                                                           )
345
346
347
348
      else
        print *,"unknown solver"
        stop
      endif
349

350
      if (present(error)) then
351
        if (success_l) then
352
          error = ELPA_OK
353
        else
354
          error = ELPA_ERROR
355
356
357
358
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
Andreas Marek's avatar
Andreas Marek committed
359

360
361
    end subroutine

362

363
    subroutine elpa_solve_real_single(self, a, ev, q, error)
364
365
      use elpa2_impl
      use elpa1_impl
366
      use elpa_utilities, only : error_unit
367
      use precision
368
      use iso_c_binding
369
      class(elpa_impl_t)  :: self
370
371
372
#ifdef USE_ASSUMED_SIZE
      real(kind=c_float)  :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
373
      real(kind=c_float)  :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
374
#endif
375
      real(kind=c_float)  :: ev(self%na)
376
377

      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
378
379
      integer, optional   :: error
      integer(kind=c_int) :: error_actual
380
      logical             :: success_l
381

382
#ifdef WANT_SINGLE_PRECISION_REAL
383

384
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
385
386
387
        success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q,  &
                                                          time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back)
388

389
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
390
391
392
393
        success_l = elpa_solve_evp_real_2stage_single_impl(self, a, ev, q,  &
                                                          time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back &
                                                          )
394
395
396
397
      else
        print *,"unknown solver"
        stop
      endif
398

399
      if (present(error)) then
400
        if (success_l) then
401
          error = ELPA_OK
402
        else
403
          error = ELPA_ERROR
404
405
406
407
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
408

409
#else
410
      print *,"This installation of the ELPA library has not been build with single-precision support"
411
      error = ELPA_ERROR
412
413
414
#endif
    end subroutine

415

416
    subroutine elpa_solve_complex_double(self, a, ev, q, error)
417
418
      use elpa2_impl
      use elpa1_impl
419
      use elpa_utilities, only : error_unit
420
      use precision
421
      use iso_c_binding
422
      class(elpa_impl_t)             :: self
423

424
425
426
#ifdef USE_ASSUMED_SIZE
      complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *)
#else
427
      complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
428
#endif
429
      real(kind=c_double)            :: ev(self%na)
430

431
432
      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back

433
434
      integer, optional              :: error
      integer(kind=c_int)            :: error_actual
435
      logical                        :: success_l
436

437
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
438
439
440
        success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q,  &
                                                          time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back)
441

442
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
443
444
445
446
        success_l = elpa_solve_evp_complex_2stage_double_impl(self,  a, ev, q,  &
                                                          time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back &
                                                          )
447
448
449
450
      else
        print *,"unknown solver"
        stop
      endif
451

452
      if (present(error)) then
453
        if (success_l) then
454
          error = ELPA_OK
455
        else
456
          error = ELPA_ERROR
457
458
459
460
461
462
463
464
        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


465
    subroutine elpa_solve_complex_single(self, a, ev, q, error)
466
467
      use elpa2_impl
      use elpa1_impl
468
469
470
      use elpa_utilities, only : error_unit

      use iso_c_binding
471
      use precision
472
      class(elpa_impl_t)            :: self
473
474
475
476
477
478
#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)
479

480
      real(kind=c_double) :: time_evp_fwd, time_evp_solve, time_evp_back
481
482
      integer, optional             :: error
      integer(kind=c_int)           :: error_actual
483
      logical                       :: success_l
484
485

#ifdef WANT_SINGLE_PRECISION_COMPLEX
486

487
      if (self%get("solver") .eq. ELPA_SOLVER_1STAGE) then
488
489
490
        success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q,  &
                                                          time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back)
491

492
      else if (self%get("solver") .eq. ELPA_SOLVER_2STAGE) then
493
494
495
496
        success_l = elpa_solve_evp_complex_2stage_single_impl(self,  a, ev, q,  &
                                                           time_evp_fwd,     &
                                                          time_evp_solve, time_evp_back, &
                                                          )
497
498
499
500
      else
        print *,"unknown solver"
        stop
      endif
501

502
      if (present(error)) then
503
        if (success_l) then
504
          error = ELPA_OK
505
        else
506
          error = ELPA_ERROR
507
508
509
510
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
      endif
511

512
#else
513
      print *,"This installation of the ELPA library has not been build with single-precision support"
514
      error = ELPA_ERROR
515
516
517
#endif
    end subroutine

518

519
    subroutine elpa_multiply_at_b_double (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
520
                                          c, ldc, ldcCols, error)
521
      use iso_c_binding
522
      use elpa1_auxiliary_impl
523
      use precision
524
      class(elpa_impl_t)              :: self
525
      character*1                     :: uplo_a, uplo_c
526
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
527
528
529
#ifdef USE_ASSUMED_SIZE
      real(kind=rk8)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
530
      real(kind=rk8)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
531
#endif
532
      integer, optional               :: error
533
534
      logical                         :: success_l

535
      success_l = elpa_mult_at_b_real_double_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
536
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
537
      if (present(error)) then
538
        if (success_l) then
539
          error = ELPA_OK
540
        else
541
          error = ELPA_ERROR
542
543
544
545
546
547
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
    end subroutine

548

549
    subroutine elpa_multiply_at_b_single (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
550
                                          c, ldc, ldcCols, error)
551
      use iso_c_binding
552
      use elpa1_auxiliary_impl
553
      use precision
554
      class(elpa_impl_t)              :: self
555
      character*1                     :: uplo_a, uplo_c
556
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
557
558
559
#ifdef USE_ASSUMED_SIZE
      real(kind=rk4)                  :: a(lda,*), b(ldb,*), c(ldc,*)
#else
560
      real(kind=rk4)                  :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
561
#endif
562
      integer, optional               :: error
563
564
      logical                         :: success_l
#ifdef WANT_SINGLE_PRECISION_REAL
565
      success_l = elpa_mult_at_b_real_single_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
566
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
567
      if (present(error)) then
568
        if (success_l) then
569
          error = ELPA_OK
570
        else
571
          error = ELPA_ERROR
572
573
574
575
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
576
577
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
578
      error = ELPA_ERROR
579
580
581
#endif
    end subroutine

582

583
    subroutine elpa_multiply_ah_b_double (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
584
                                          c, ldc, ldcCols, error)
585
      use iso_c_binding
586
      use elpa1_auxiliary_impl
587
      use precision
588
      class(elpa_impl_t)              :: self
589
      character*1                     :: uplo_a, uplo_c
590
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
591
592
593
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck8)               :: a(lda,*), b(ldb,*), c(ldc,*)
#else
594
      complex(kind=ck8)               :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
595
#endif
596
      integer, optional               :: error
597
598
      logical                         :: success_l

599
      success_l = elpa_mult_ah_b_complex_double_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
600
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
601
      if (present(error)) then
602
        if (success_l) then
603
          error = ELPA_OK
604
        else
605
          error = ELPA_ERROR
606
607
608
609
610
611
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
    end subroutine

612

613
    subroutine elpa_multiply_ah_b_single (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
614
                                          c, ldc, ldcCols, error)
615
      use iso_c_binding
616
      use elpa1_auxiliary_impl
617
      use precision
618
      class(elpa_impl_t)              :: self
619
      character*1                     :: uplo_a, uplo_c
620
      integer(kind=ik), intent(in)    :: na, lda, ldaCols, ldb, ldbCols, ldc, ldcCols, ncb
621
622
623
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck4)               :: a(lda,*), b(ldb,*), c(ldc,*)
#else
624
      complex(kind=ck4)               :: a(lda,ldaCols), b(ldb,ldbCols), c(ldc,ldcCols)
625
#endif
626
      integer, optional               :: error
627
628
629
      logical                         :: success_l

#ifdef WANT_SINGLE_PRECISION_COMPLEX
630
      success_l = elpa_mult_ah_b_complex_single_impl(uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, self%nblk, &
631
                              self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), c, ldc, ldcCols)
632
      if (present(error)) then
633
        if (success_l) then
634
          error = ELPA_OK
635
        else
636
          error = ELPA_ERROR
637
638
639
640
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in multiply_a_b() and you did not check for errors!"
      endif
641
642
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
643
      error = ELPA_ERROR
644
645
646
#endif
    end subroutine

647

648
    subroutine elpa_cholesky_double_real (self, a, error)
649
      use iso_c_binding
650
      use elpa1_auxiliary_impl
651
      use precision
652
      class(elpa_impl_t)              :: self
653
654
655
#ifdef USE_ASSUMED_SIZE
      real(kind=rk8)                  :: a(self%local_nrows,*)
#else
656
      real(kind=rk8)                  :: a(self%local_nrows,self%local_ncols)
657
#endif
658
      integer, optional               :: error
659
      logical                         :: success_l
660
      integer(kind=c_int)             :: error_actual
661
662
      logical                         :: wantDebugIntern

663
664
      if (self%get("wantDebug",error_actual) .eq. 1) then
        if (error_actual .ne. ELPA_OK) then
665
666
667
668
669
670
671
672
673
          print *,"Could not querry wantDebug"
          stop
        endif

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

674

675
      success_l = elpa_cholesky_real_double_impl (self%na, a, self%local_nrows, self%nblk, &
676
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
677
                                                 wantDebugIntern)
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
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine

689

690
    subroutine elpa_cholesky_single_real(self, a, error)
691
      use iso_c_binding
692
      use elpa1_auxiliary_impl
693
      use precision
694
      class(elpa_impl_t)              :: self
695
696
697
#ifdef USE_ASSUMED_SIZE
      real(kind=rk4)                  :: a(self%local_nrows,*)
#else
698
      real(kind=rk4)                  :: a(self%local_nrows,self%local_ncols)
699
#endif
700
      integer, optional               :: error
701
      logical                         :: success_l
702
      integer(kind=c_int)             :: error_actual
703
704
      logical                         :: wantDebugIntern

705
706
      if (self%get("wantDebug",error_actual) .eq. 1) then
        if (error_actual .ne. ELPA_OK) then
707
708
709
710
711
712
713
714
715
716
          print *,"Could not querry wantDebug"
          stop
        endif

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

#if WANT_SINGLE_PRECISION_REAL
717
      success_l = elpa_cholesky_real_single_impl (self%na, a, self%local_nrows, self%nblk, &
718
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
719
                                                 wantDebugIntern)
720
721
#else
      print *,"This installation of the ELPA library has not been build with single-precision support"
722
      error = ELPA_ERROR
723
#endif
724
      if (present(error)) then
725
        if (success_l) then
726
          error = ELPA_OK
727
        else
728
          error = ELPA_ERROR
729
730
731
732
733
734
        endif
      else if (.not. success_l) then
        write(error_unit,'(a)') "ELPA: Error in cholesky() and you did not check for errors!"
      endif
    end subroutine

735

736
    subroutine elpa_cholesky_double_complex (self, a, error)
737
      use iso_c_binding
738
      use elpa1_auxiliary_impl
739
      use precision
740
      class(elpa_impl_t)              :: self
741
742
743
744
745
#ifdef USE_ASSUMED_SIZE
      complex(kind=ck8)               :: a(self%local_nrows,*)
#else
      complex(kind=ck8)               :: a(self%local_nrows,self%local_ncols)
#endif
746
      integer, optional               :: error
747
      logical                         :: success_l
748
      integer(kind=c_int)             :: error_actual
749
750
      logical                         :: wantDebugIntern

751
752
      if (self%get("wantDebug",error_actual) .eq. 1) then
        if (error_actual .ne. ELPA_OK) then
753
754
755
756
757
758
759
760
761
          print *,"Could not querry wantDebug"
          stop
        endif

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

762
      success_l = elpa_cholesky_complex_double_impl (self%na, a, self%local_nrows, self%nblk, &
763
                                                 self%local_ncols, self%get("mpi_comm_rows"), self%get("mpi_comm_cols"), &
764
                                                 wantDebugIntern)
765
      if (present(error)) then
766
        if (success_l) then
767
          error = ELPA_OK