Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : MODULE fft_lib
       8             : 
       9             :    USE fft_kinds,                       ONLY: dp
      10             :    USE fft_plan,                        ONLY: fft_plan_type
      11             :    USE fftsg_lib,                       ONLY: fftsg1dm,&
      12             :                                               fftsg3d,&
      13             :                                               fftsg_do_cleanup,&
      14             :                                               fftsg_do_init,&
      15             :                                               fftsg_get_lengths
      16             :    USE fftw3_lib,                       ONLY: &
      17             :         fft_alloc => fftw_alloc, fft_dealloc => fftw_dealloc, fftw31dm, fftw33d, &
      18             :         fftw3_create_plan_1dm, fftw3_create_plan_3d, fftw3_destroy_plan, fftw3_do_cleanup, &
      19             :         fftw3_do_init, fftw3_get_lengths
      20             : #include "../../base/base_uses.f90"
      21             : 
      22             :    IMPLICIT NONE
      23             :    PRIVATE
      24             : 
      25             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fft_lib'
      26             : 
      27             :    PUBLIC :: fft_do_cleanup, fft_do_init, fft_get_lengths, fft_create_plan_3d
      28             :    PUBLIC :: fft_create_plan_1dm, fft_1dm, fft_library, fft_3d, fft_destroy_plan
      29             :    PUBLIC :: fft_alloc, fft_dealloc
      30             : 
      31             : CONTAINS
      32             : ! **************************************************************************************************
      33             : !> \brief Interface to FFT libraries
      34             : !> \param fftlib ...
      35             : !> \return ...
      36             : !> \par History
      37             : !>      IAB 09-Jan-2009 : Modified to use fft_plan_type
      38             : !>                        (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
      39             : !> \author JGH
      40             : ! **************************************************************************************************
      41        8961 :    FUNCTION fft_library(fftlib) RESULT(flib)
      42             : 
      43             :       CHARACTER(len=*), INTENT(IN)                       :: fftlib
      44             :       INTEGER                                            :: flib
      45             : 
      46             :       SELECT CASE (fftlib)
      47             :       CASE DEFAULT
      48          10 :          flib = -1
      49             :       CASE ("FFTSG")
      50          10 :          flib = 1
      51             :       CASE ("FFTW3")
      52        8961 :          flib = 3
      53             :       END SELECT
      54             : 
      55        8961 :    END FUNCTION fft_library
      56             : 
      57             : ! **************************************************************************************************
      58             : !> \brief ...
      59             : !> \param fft_type ...
      60             : !> \param wisdom_file ...
      61             : ! **************************************************************************************************
      62        8961 :    SUBROUTINE fft_do_init(fft_type, wisdom_file)
      63             :       INTEGER, INTENT(IN)                                :: fft_type
      64             :       CHARACTER(LEN=*), INTENT(IN)                       :: wisdom_file
      65             : 
      66        8961 :       SELECT CASE (fft_type)
      67             :       CASE DEFAULT
      68           0 :          CPABORT("fft_do_init")
      69             :       CASE (1)
      70          10 :          CALL fftsg_do_init()
      71             :       CASE (3)
      72        8961 :          CALL fftw3_do_init(wisdom_file)
      73             :       END SELECT
      74             : 
      75        8961 :    END SUBROUTINE
      76             : 
      77             : ! **************************************************************************************************
      78             : !> \brief ...
      79             : !> \param fft_type ...
      80             : !> \param wisdom_file ...
      81             : !> \param ionode ...
      82             : ! **************************************************************************************************
      83        8751 :    SUBROUTINE fft_do_cleanup(fft_type, wisdom_file, ionode)
      84             :       INTEGER, INTENT(IN)                                :: fft_type
      85             :       CHARACTER(LEN=*), INTENT(IN)                       :: wisdom_file
      86             :       LOGICAL, INTENT(IN)                                :: ionode
      87             : 
      88        8751 :       SELECT CASE (fft_type)
      89             :       CASE DEFAULT
      90           0 :          CPABORT("fft_do_cleanup")
      91             :       CASE (1)
      92          10 :          CALL fftsg_do_cleanup()
      93             :       CASE (3)
      94        8751 :          CALL fftw3_do_cleanup(wisdom_file, ionode)
      95             :       END SELECT
      96             : 
      97        8751 :    END SUBROUTINE
      98             : 
      99             : ! **************************************************************************************************
     100             : !> \brief ...
     101             : !> \param fft_type ...
     102             : !> \param DATA ...
     103             : !> \param max_length ...
     104             : ! **************************************************************************************************
     105      107912 :    SUBROUTINE fft_get_lengths(fft_type, DATA, max_length)
     106             : 
     107             :       INTEGER, INTENT(IN)                                :: fft_type
     108             :       INTEGER, DIMENSION(*)                              :: DATA
     109             :       INTEGER, INTENT(INOUT)                             :: max_length
     110             : 
     111      107912 :       SELECT CASE (fft_type)
     112             :       CASE DEFAULT
     113           0 :          CPABORT("fft_get_lengths")
     114             :       CASE (1)
     115      107408 :          CALL fftsg_get_lengths(DATA, max_length)
     116             :       CASE (3)
     117      107912 :          CALL fftw3_get_lengths(DATA, max_length)
     118             :       END SELECT
     119             : 
     120      107912 :    END SUBROUTINE fft_get_lengths
     121             : 
     122             : ! **************************************************************************************************
     123             : 
     124             : ! **************************************************************************************************
     125             : !> \brief ...
     126             : !> \param plan ...
     127             : !> \param fft_type ...
     128             : !> \param fft_in_place ...
     129             : !> \param fsign ...
     130             : !> \param n ...
     131             : !> \param zin ...
     132             : !> \param zout ...
     133             : !> \param plan_style ...
     134             : ! **************************************************************************************************
     135       40756 :    SUBROUTINE fft_create_plan_3d(plan, fft_type, fft_in_place, fsign, n, zin, zout, plan_style)
     136             : 
     137             :       TYPE(fft_plan_type), INTENT(INOUT)                 :: plan
     138             :       INTEGER, INTENT(IN)                                :: fft_type
     139             :       LOGICAL, INTENT(IN)                                :: fft_in_place
     140             :       INTEGER, INTENT(IN)                                :: fsign
     141             :       INTEGER, DIMENSION(3), INTENT(IN)                  :: n
     142             :       COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT)      :: zin, zout
     143             :       INTEGER, INTENT(IN)                                :: plan_style
     144             : 
     145       40756 :       plan%fft_type = fft_type
     146       40756 :       plan%fsign = fsign
     147       40756 :       plan%fft_in_place = fft_in_place
     148      163024 :       plan%n_3d = n
     149       40756 : !$    plan%need_alt_plan = .FALSE.
     150             : 
     151             :       ! Planning only needed for FFTW3
     152       40756 :       IF (fft_type .EQ. 3) THEN
     153       40684 :          CALL fftw3_create_plan_3d(plan, zin, zout, plan_style)
     154       40684 :          plan%valid = .TRUE.
     155             :       END IF
     156             : 
     157       40756 :    END SUBROUTINE fft_create_plan_3d
     158             : 
     159             : !
     160             : ! really ugly, plan is intent out, because plan%fsign is also a status flag
     161             : ! if something goes wrong, plan%fsign is set to zero, and the plan becomes invalid
     162             : !
     163             : ! **************************************************************************************************
     164             : !> \brief ...
     165             : !> \param plan ...
     166             : !> \param scale ...
     167             : !> \param zin ...
     168             : !> \param zout ...
     169             : !> \param stat ...
     170             : ! **************************************************************************************************
     171      614288 :    SUBROUTINE fft_3d(plan, scale, zin, zout, stat)
     172             :       TYPE(fft_plan_type), INTENT(IN)                    :: plan
     173             :       REAL(KIND=dp), INTENT(IN)                          :: scale
     174             :       COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT)      :: zin, zout
     175             :       INTEGER, INTENT(OUT)                               :: stat
     176             : 
     177      614288 :       stat = plan%fsign
     178      614288 :       IF (plan%n_3d(1)*plan%n_3d(2)*plan%n_3d(3) > 0) THEN
     179      614288 :          SELECT CASE (plan%fft_type)
     180             :          CASE DEFAULT
     181           0 :             CPABORT("fft_3d")
     182             :          CASE (1)
     183         886 :             CALL fftsg3d(plan%fft_in_place, stat, scale, plan%n_3d, zin, zout)
     184             :          CASE (3)
     185      614288 :             CALL fftw33d(plan, scale, zin, zout, stat)
     186             :          END SELECT
     187             :       END IF
     188             :       ! stat is set to zero on error, -1,+1 are OK
     189      614288 :       IF (stat .EQ. 0) THEN
     190           0 :          stat = 1
     191             :       ELSE
     192      614288 :          stat = 0
     193             :       END IF
     194             : 
     195      614288 :    END SUBROUTINE fft_3d
     196             : 
     197             : ! **************************************************************************************************
     198             : 
     199             : ! **************************************************************************************************
     200             : !> \brief ...
     201             : !> \param plan ...
     202             : !> \param fft_type ...
     203             : !> \param fsign ...
     204             : !> \param trans ...
     205             : !> \param n ...
     206             : !> \param m ...
     207             : !> \param zin ...
     208             : !> \param zout ...
     209             : !> \param plan_style ...
     210             : ! **************************************************************************************************
     211      149020 :    SUBROUTINE fft_create_plan_1dm(plan, fft_type, fsign, trans, n, m, zin, zout, plan_style)
     212             :       TYPE(fft_plan_type), INTENT(INOUT)                 :: plan
     213             :       INTEGER, INTENT(IN)                                :: fft_type, fsign
     214             :       LOGICAL, INTENT(IN)                                :: trans
     215             :       INTEGER, INTENT(IN)                                :: n, m
     216             :       COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN)         :: zin, zout
     217             :       INTEGER, INTENT(IN)                                :: plan_style
     218             : 
     219      149020 :       plan%fft_type = fft_type
     220      149020 :       plan%fsign = fsign
     221      149020 :       plan%trans = trans
     222      149020 :       plan%n = n
     223      149020 :       plan%m = m
     224      149020 : !$    plan%need_alt_plan = .FALSE.
     225             : 
     226             :       ! Planning only needed for FFTW3
     227      149020 :       IF ((fft_type .EQ. 3) .AND. (n*m .NE. 0)) THEN
     228      148852 :          CALL fftw3_create_plan_1dm(plan, zin, zout, plan_style)
     229      148852 :          plan%valid = .TRUE.
     230             :       ELSE
     231         168 :          plan%valid = .FALSE.
     232             :       END IF
     233             : 
     234      149020 :    END SUBROUTINE fft_create_plan_1dm
     235             : 
     236             : ! **************************************************************************************************
     237             : !> \brief ...
     238             : !> \param plan ...
     239             : ! **************************************************************************************************
     240      263920 :    SUBROUTINE fft_destroy_plan(plan)
     241             :       TYPE(fft_plan_type), INTENT(INOUT)                 :: plan
     242             : 
     243             : ! Planning only needed for FFTW3
     244             : 
     245      263920 :       IF (plan%valid .AND. plan%fft_type .EQ. 3) THEN
     246      189536 :          CALL fftw3_destroy_plan(plan)
     247      189536 :          plan%valid = .FALSE.
     248             :       END IF
     249             : 
     250      263920 :    END SUBROUTINE
     251             : 
     252             : ! **************************************************************************************************
     253             : !> \brief ...
     254             : !> \param plan ...
     255             : !> \param zin ...
     256             : !> \param zout ...
     257             : !> \param scale ...
     258             : !> \param stat ...
     259             : ! **************************************************************************************************
     260     7846906 :    SUBROUTINE fft_1dm(plan, zin, zout, scale, stat)
     261             :       TYPE(fft_plan_type), INTENT(IN)                    :: plan
     262             :       COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT)      :: zin, zout
     263             :       REAL(KIND=dp), INTENT(IN)                          :: scale
     264             :       INTEGER, INTENT(OUT)                               :: stat
     265             : 
     266     7846906 :       stat = plan%fsign
     267     7846906 :       IF (plan%n*plan%m > 0) THEN
     268     7846906 :          SELECT CASE (plan%fft_type)
     269             :          CASE DEFAULT
     270           0 :             CPABORT("fft_1dm")
     271             :          CASE (1)
     272       15522 :             CALL fftsg1dm(stat, plan%trans, plan%n, plan%m, zin, zout, scale)
     273             :          CASE (3)
     274     7846906 :             CALL fftw31dm(plan, zin, zout, scale, stat)
     275             :          END SELECT
     276             :       END IF
     277             :       ! stat is set to zero on error, -1,+1 are OK
     278     7846906 :       IF (stat .EQ. 0) THEN
     279           0 :          stat = 1
     280             :       ELSE
     281     7846906 :          stat = 0
     282             :       END IF
     283             : 
     284     7846906 :    END SUBROUTINE fft_1dm
     285             : 
     286             : END MODULE
     287             :