ifft235.f 4.35 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

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 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
      MODULE m_ifft
      use m_juDFT
      CONTAINS
      INTEGER FUNCTION  ifft235 (iofile,ksfft,n,gmaxp)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C   this function checks wether n can be expressed as :
C   n =  (2**P) * (3**Q) * (5**R) to match withe the MFFT - routines
C   used in this program. If close to n there is a number n' which
c   is solely expressable as n'= 2**p then this number sis choosen
c   and is outputted by the function. If n is not expressable as
C   n =  (2**P) * (3**Q) * (5**R) a number n' is choosen close to ns
c   is selected which fullfilles these requirements.
c
c
c                            Stefan Bl"ugel , kfa, Oct. 1993
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IMPLICIT NONE
c
c----> declaration part
c
      INTEGER iofile,ksfft,n
      REAL    gmaxp
c
c----> local variabels
c
      INTEGER np,nn,nl,iexp,is,itrial
      REAL rl
      REAL fac_1,two,fac_2,fac_3,one,fac_4
      PARAMETER ( two=2.0e0,fac_1=1.001e0,fac_2=1.7e0,fac_3=1.95e0 )
      PARAMETER ( one=1.0e0,fac_4=0.03e0 )
      INTRINSIC abs,log,mod,real
c
c
      ifft235=0
      IF ( n.eq.0 ) THEN
         WRITE (iofile,*) 'n should not be zero'
          CALL juDFT_error("n should not be zero",calledby="ifft235")
      ENDIF
c
c====> RADIX 2 CALCULATION ONLY
c
      IF ( ksfft .EQ. 0 ) THEN
          rl   = log(real(n))/log(two)
          nl   = rl
          IF ( abs( n - 2 ** nl ) .GT. abs(n - 2 ** (nl+1) ) ) THEN
               np = nl + 1
          ELSE
               np = nl
          ENDIF
          IF ( ( np .eq. nl ) .AND. ( gmaxp * real(nl)/rl .lt. fac_1 ) )
     >         np = nl + 1
          ifft235 = 2 ** np
          RETURN
       ENDIF
c
c====> RADIX 2 , 3, 5 CALCULATION
c
c
c----> since gmaxp is large enough, try also smaller n
c
      IF ( gmaxp .GT. fac_2 .AND. gmaxp .LT. fac_3 ) THEN
           is = -1
      ELSE
           is =  1
      ENDIF
      np   = n
      nn   = n
c
      DO 20 itrial = 1 , 200
c
c----> check whether there is a number close by, which is 2**NL
c      ( FFT very fast )
c
          rl   = log(real(nn))/log(two)
          nl   = rl
          IF (( abs( real(nl)/rl - one ) .LT. fac_4 ).AND.
     +        ( 2**nl.GE.n ) ) THEN
             IF (gmaxp * real(nl)/rl .GT. one) ifft235 = 2 ** nl
             RETURN
          ELSE IF (( abs( real(nl+1)/rl - one ) .lt. fac_4 ).AND.
     +             ( 2**(nl+1).GE.n ) ) THEN
             ifft235 = 2 ** (nl + 1)
             RETURN
          ELSE
c
c----> no , no binary number is arround, check wether number can be
c      divided by 2,3 or 5
c
             DO 10 iexp = 1 , nl+1
                IF (mod(nn,2) .EQ. 0) THEN
                   nn = nn / 2
                ELSE IF (mod(nn,3) .EQ. 0 ) THEN
                   nn = nn / 3
                ELSE IF (mod(nn,5) .EQ. 0 ) THEN
                   nn = nn / 5
                ENDIF
                IF ( nn.eq.2 .OR. nn.eq.3 .OR. nn.eq.5 ) THEN
c---> o.k.
                   ifft235 = np
                   RETURN
                ENDIF
  10         ENDDO
c
c----> no , n  cannot be expressed as (2**P) * (3**Q) * (5**R)
c      change number
             nn = n + (is ** (itrial)) * ((itrial + 1)/2)
             np = n + (is ** (itrial)) * ((itrial + 1)/2)
          ENDIF
  20    ENDDO
c
      WRITE (iofile,'('' to few trials '', i3)') itrial
      STOP
c
      END FUNCTION ifft235
C-----------------------------------------------------------------------
      INTEGER FUNCTION i2357(ii)

!
! simple setup to determine fft-length for ESSL calls
!
      IMPLICIT NONE
      INTEGER, INTENT (IN)  :: ii
      INTEGER i,j,k,h,m,n,nn

      nn = 12582912
      DO m = 0,1 
        DO k = 0,1 
          DO j = 0,1 
            DO i = 0,1 
              DO h = 1,25
                n = 2**h * 3**i * 5**j * 7**k * 11**m
                IF ( (n.GT.ii).AND.(n.LT.nn) ) nn = n
              ENDDO
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      i2357 = nn

      END FUNCTION i2357
 
      END MODULE m_ifft