LCOV - code coverage report
Current view: top level - math - cfft.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 464 466 99.6 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : #ifndef CPP_ESSL
       2             : #define CPP_Singleton
       3             : #endif
       4             : ! ou might want to undefine CPP_singleton to use essl FFT
       5             : MODULE m_cfft
       6             :    IMPLICIT NONE
       7             : CONTAINS
       8             : !     ***************************************************************
       9             : !     multivariate complex fourier transform, computed in place
      10             : !     using mixed-radix fast fourier transform algorithm.
      11             : !     by r. c. singleton, stanford research institute, oct. 1968
      12             : !     arrays a and b originally hold the real and imaginary
      13             : !     components of the data, and return the real and
      14             : !     imaginary components of the resulting fourier coefficients.
      15             : !     multivariate data is indexed according to the fortran
      16             : !     array element successor function, without limit
      17             : !     on the number of implied multiple subscripts.
      18             : !     the subroutine is called once for each variate.
      19             : !     the calls for a multivariate transform may be in any order.
      20             : !     ntot is the total number of complex data values.
      21             : !     n is the dimension of the current variable.
      22             : !     nspan/n is the spacing of consucutive data values
      23             : !     while indexing the current variable.
      24             : !     the sign of isn determines the sign of the complex
      25             : !     exponential, and the magnitude of isn is normally one.
      26             : !     for a single-variate transform,
      27             : !     ntot = n = nspan = (number of complex data values), f.g.
      28             : !     call cft(a,b,n,n,n,1)
      29             : !     a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
      30             : !     is computed by
      31             : !     call cft(a,b,n1*n2*n3,n1,n1,1)
      32             : !     call cft(a,b,n1*n2*n3,n2,n1*n2,1)
      33             : !     call cft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
      34             : !     the data may alternatively be stored in a single complex
      35             : !     array a, then the magnitude of isn changed to two to
      36             : !     give the correct indexing increment and the second parameter
      37             : !     used to pass the initial address for the sequence of
      38             : !     imaginary values, e.g.
      39             : !        real s(2)
      40             : !        equivalence (a,s)
      41             : !        ....
      42             : !        ....
      43             : !        call cft(a,s(2),ntot,n,nspan,2)
      44             : !     arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
      45             : !     are used for temporary storage. if the available storage
      46             : !     is insufficient, the program is terminated by a stop.
      47             : !     maxf must be .ge. the maximum prime factor of n.
      48             : !     maxp must be .gt. the number of prime factors of n.
      49             : !     in addition, if the square-free portion k of n has two or
      50             : !     more prime factors, then maxp must be .ge. k-1.
      51             : !     array storage in nfac for a maximum of 11 factors of n.
      52             : !     if n has more than one square-free factor, the product of the
      53             : !     square-free factors must be .le. 210
      54             : !     *******************************************************************
      55             : !     array storage for maximum prime factor of 199
      56             : !     the following two constants should agree with the array dimensions
      57             : #ifdef CPP_Singleton
      58      862943 :    SUBROUTINE cfft(a, b, ntot, n, nspan, isn)
      59             :       use m_juDFT
      60             :       IMPLICIT NONE
      61             : !     .. Scalar Arguments ..
      62             :       INTEGER :: isn, n, nspan, ntot
      63             : !     ..
      64             : !     .. Array Arguments ..
      65             : !      REAL a(*),b(*)
      66             :       REAL :: a(ntot), b(ntot)
      67             : !     ..
      68             : !     .. Local Scalars ..
      69             :       REAL :: aa, aj, ajm, ajp, ak, akm, akp, bb, bj, bjm, bjp, bk, bkm, bkp, c1, c2, c3, &
      70             :               c72, cd, rad, radf, s1, s120, s2, s3, s72, sd
      71             :       INTEGER :: i, ii, inc, j, jc, jf, jj, k, k1, k2, k3, k4, kk, ks, kspan, kspnn, kt, m, &
      72             :                  maxf, maxp, nn, nt, maxnf
      73             : !     ..
      74             : !     .. Local Arrays ..
      75     1725886 :       REAL, ALLOCATABLE :: at(:), bt(:), ck(:), sk(:)
      76      862943 :       INTEGER, ALLOCATABLE :: nfac(:), np(:)
      77             : !     ..
      78             : !     .. Intrinsic Functions ..
      79             :       INTRINSIC cos, real, mod, sin, sqrt
      80             : !     ..
      81             : !     .. Equivalences ..
      82             :       EQUIVALENCE(i, ii)
      83             : !     ..
      84      862943 :       IF (n < 2) RETURN
      85             : 
      86      862943 :       maxf = MAX(299, n)
      87      862943 :       maxp = MAX(503, n - 1)
      88      862943 :       maxnf = MAX(17, CEILING(sqrt(1.0*n)))
      89      862943 :       ALLOCATE (at(maxf), bt(maxf), ck(maxf), sk(maxf))
      90      862943 :       ALLOCATE (nfac(maxnf), np(maxp))
      91             : 
      92      862943 :       inc = isn
      93             : !     the following constants are rad = 2.*pi , s72 = sin(0.4*pi) ,
      94             : !     c72 = cos(0.4*pi) and s120 = sqrt(0.75)
      95      862943 :       rad = 6.2831853071796
      96      862943 :       s72 = 0.95105651629515
      97      862943 :       c72 = 0.30901699437495
      98      862943 :       s120 = 0.86602540378444
      99      862943 :       IF (isn < 0) THEN
     100       48520 :          s72 = -s72
     101       48520 :          s120 = -s120
     102       48520 :          rad = -rad
     103       48520 :          inc = -inc
     104             :       ENDIF
     105      862943 :       nt = inc*ntot
     106      862943 :       ks = inc*nspan
     107      862943 :       kspan = ks
     108      862943 :       nn = nt - inc
     109      862943 :       jc = ks/n
     110      862943 :       radf = rad*real(jc)*0.5
     111             :       i = 0
     112      862943 :       jf = 0
     113             : !     determine the factors of n
     114      862943 :       m = 0
     115      862943 :       k = n
     116     2246769 :       DO WHILE (k - (k/16)*16 == 0)
     117      691913 :          m = m + 1
     118      691913 :          nfac(m) = 4
     119     1554856 :          k = k/16
     120             :       ENDDO
     121             :       j = 3
     122             :       jj = 9
     123       27830 :       GO TO 50
     124       27830 :    40 m = m + 1
     125       27830 :       nfac(m) = j
     126       47652 :       k = k/jj
     127      910595 :    50 IF (mod(k,jj).EQ.0) GO TO 40
     128      882765 :       j = j + 2
     129      882765 :       jj = j**2
     130      882765 :       IF (jj.LE.k) GO TO 50
     131      862943 :       IF (k.GT.4) GO TO 60
     132      718448 :       kt = m
     133      718448 :       nfac(m+1) = k
     134      718448 :       IF (k.NE.1) m = m + 1
     135      144495 :       GO TO 100
     136      144495 :    60 IF (k- (k/4)*4.NE.0) GO TO 70
     137       98618 :       m = m + 1
     138       98618 :       nfac(m) = 2
     139      144495 :       k = k/4
     140      144495 :    70 kt = m
     141      393377 :       j = 2
     142      393377 :    80 IF (mod(k,j).NE.0) GO TO 90
     143      201501 :       m = m + 1
     144      201501 :       nfac(m) = j
     145      393377 :       k = k/j
     146      393377 : 90    j = ((j + 1)/2)*2 + 1
     147      393377 :       IF (j <= k) GO TO 80
     148      862943 : 100   IF (kt == 0) GO TO 120
     149     2455049 :       DO j=kt,1,-1
     150      818361 :          m = m + 1
     151     1681304 :          nfac(m) = nfac(j)
     152             :       ENDDO
     153             : !     compute fourier transform
     154     1865028 : 120   sd = radf/real(kspan)
     155     1865028 :       cd = 2.0*sin(sd)**2
     156     1865028 :       sd = sin(sd + sd)
     157     1865028 :       kk = 1
     158     1865028 :       i = i + 1
     159     1865028 :       IF (nfac(i) /= 2) GO TO 170
     160             : !     transform for factor of 2 (including rotation factor)
     161      225384 :       kspan = kspan/2
     162     4887391 :       k1 = kspan + 2
     163    79933114 : 130   k2 = kk + kspan
     164    79933114 :       ak = a(k2)
     165    79933114 :       bk = b(k2)
     166    79933114 :       a(k2) = a(kk) - ak
     167    79933114 :       b(k2) = b(kk) - bk
     168    79933114 :       a(kk) = a(kk) + ak
     169    79933114 :       b(kk) = b(kk) + bk
     170    79933114 :       kk = k2 + kspan
     171    79933114 :       IF (kk <= nn) GO TO 130
     172     4887391 :       kk = kk - nn
     173     4887391 :       IF (kk <= jc) GO TO 130
     174      225384 :       IF (kk > kspan) GO TO 360
     175     3210254 : 140   c1 = 1.0 - cd
     176    21254546 :       s1 = sd
     177   101434774 : 150   k2 = kk + kspan
     178   101434774 :       ak = a(kk) - a(k2)
     179   101434774 :       bk = b(kk) - b(k2)
     180   101434774 :       a(kk) = a(kk) + a(k2)
     181   101434774 :       b(kk) = b(kk) + b(k2)
     182   101434774 :       a(k2) = c1*ak - s1*bk
     183   101434774 :       b(k2) = s1*ak + c1*bk
     184   101434774 :       kk = k2 + kspan
     185   101434774 :       IF (kk < nt) GO TO 150
     186    21254546 :       k2 = kk - nt
     187    21254546 :       c1 = -c1
     188    21254546 :       kk = k1 - k2
     189    21254546 :       IF (kk > k2) GO TO 150
     190    11646992 :       ak = c1 - (cd*c1 + sd*s1)
     191    11646992 :       s1 = (sd*c1 - cd*s1) + s1
     192             : !     the following three statements compensate for truncation
     193             : !     error. if rounded arithmetic is used, they may be deleted.
     194             : !     c1=0.5/(ak**2+s1**2)+0.5
     195             : !     s1=c1*s1
     196             : !     c1=c1*ak
     197             : !     next statement should be deleted if non-rounded arithmetic is used
     198    11646992 :       c1 = ak
     199    11646992 :       kk = kk + jc
     200    11646992 :       IF (kk < k2) GO TO 150
     201     3210254 :       k1 = k1 + inc + inc
     202     3210254 :       kk = (k1 - kspan)/2 + jc
     203     3210254 :       IF (kk <= jc + jc) GO TO 140
     204   320558498 :       GO TO 120
     205             : !     transform for factor of 3 (optional code)
     206   320558498 : 160   k1 = kk + kspan
     207   320558498 :       k2 = k1 + kspan
     208   320558498 :       ak = a(kk)
     209   320558498 :       bk = b(kk)
     210   320558498 :       aj = a(k1) + a(k2)
     211   320558498 :       bj = b(k1) + b(k2)
     212   320558498 :       a(kk) = ak + aj
     213   320558498 :       b(kk) = bk + bj
     214   320558498 :       ak = -0.5*aj + ak
     215   320558498 :       bk = -0.5*bj + bk
     216   320558498 :       aj = (a(k1) - a(k2))*s120
     217   320558498 :       bj = (b(k1) - b(k2))*s120
     218   320558498 :       a(k1) = ak - bj
     219   320558498 :       b(k1) = bk + aj
     220   320558498 :       a(k2) = ak + bj
     221   320558498 :       b(k2) = bk - aj
     222   320558498 :       kk = k2 + kspan
     223   320558498 :       IF (kk < nn) GO TO 160
     224    80467106 :       kk = kk - nn
     225    80467106 :       IF (kk <= kspan) GO TO 160
     226     1639644 :       GO TO 320
     227             : !     transform for factor of 4
     228     1639644 : 170   IF (nfac(i) /= 4) GO TO 260
     229     1389927 :       kspnn = kspan
     230    13520440 :       kspan = kspan/4
     231    13520440 : 180   c1 = 1.0
     232    83740300 :       s1 = 0
     233   107386457 : 190   k1 = kk + kspan
     234   107386457 :       k2 = k1 + kspan
     235   107386457 :       k3 = k2 + kspan
     236   107386457 :       akp = a(kk) + a(k2)
     237   107386457 :       akm = a(kk) - a(k2)
     238   107386457 :       ajp = a(k1) + a(k3)
     239   107386457 :       ajm = a(k1) - a(k3)
     240   107386457 :       a(kk) = akp + ajp
     241   107386457 :       ajp = akp - ajp
     242   107386457 :       bkp = b(kk) + b(k2)
     243   107386457 :       bkm = b(kk) - b(k2)
     244   107386457 :       bjp = b(k1) + b(k3)
     245   107386457 :       bjm = b(k1) - b(k3)
     246   107386457 :       b(kk) = bkp + bjp
     247   107386457 :       bjp = bkp - bjp
     248   107386457 :       IF (isn < 0) GO TO 220
     249    87899985 :       akp = akm - bjm
     250    87899985 :       akm = akm + bjm
     251    87899985 :       bkp = bkm + ajm
     252    87899985 :       bkm = bkm - ajm
     253    87899985 :       IF (s1 == 0.0) GO TO 230
     254    45671811 : 200   a(k1) = akp*c1 - bkp*s1
     255    45671811 :       b(k1) = akp*s1 + bkp*c1
     256    45671811 :       a(k2) = ajp*c2 - bjp*s2
     257    45671811 :       b(k2) = ajp*s2 + bjp*c2
     258    45671811 :       a(k3) = akm*c3 - bkm*s3
     259    45671811 :       b(k3) = akm*s3 + bkm*c3
     260    45671811 :       kk = k3 + kspan
     261    45671811 :       IF (kk <= nt) GO TO 190
     262    35546094 : 210   c2 = c1 - (cd*c1 + sd*s1)
     263    35546094 :       s1 = (sd*c1 - cd*s1) + s1
     264             : !     the following three statements compensate for truncation
     265             : !     error. if rounded arithmetic is used, they may be deleted.
     266             : !     c1=0.5/(c2**2+s1**2)+0.5
     267             : !     s1=c1*s1
     268             : !     c1=c1*c2
     269             : !     next statement should be deleted if non-rounded arithmetic is used
     270    35546094 :       c1 = c2
     271    35546094 :       c2 = c1**2 - s1**2
     272    35546094 :       s2 = 2.0*c1*s1
     273    35546094 :       c3 = c2*c1 - s2*s1
     274    35546094 :       s3 = c2*s1 + s2*c1
     275    35546094 :       kk = kk - nt + jc
     276    35546094 :       IF (kk <= kspan) GO TO 190
     277    13520440 :       kk = kk - kspan + inc
     278    13520440 :       IF (kk <= jc) GO TO 180
     279     1389927 :       IF (kspan == jc) GO TO 360
     280    19486472 :       GO TO 120
     281    19486472 : 220   akp = akm + bjm
     282    19486472 :       akm = akm - bjm
     283    19486472 :       bkp = bkm - ajm
     284    19486472 :       bkm = bkm + ajm
     285    19486472 :       IF (s1 /= 0.0) GO TO 200
     286    61714646 : 230   a(k1) = akp
     287    61714646 :       b(k1) = bkp
     288    61714646 :       a(k2) = ajp
     289    61714646 :       b(k2) = bjp
     290    61714646 :       a(k3) = akm
     291    61714646 :       b(k3) = bkm
     292    61714646 :       kk = k3 + kspan
     293    61714646 :       IF (kk <= nt) GO TO 190
     294       12861 :       GO TO 210
     295             : !     transform for factor of 5 (optional code)
     296       12861 : 240   c2 = c72**2 - s72**2
     297     3843730 :       s2 = 2.0*c72*s72
     298    11917378 : 250   k1 = kk + kspan
     299    11917378 :       k2 = k1 + kspan
     300    11917378 :       k3 = k2 + kspan
     301    11917378 :       k4 = k3 + kspan
     302    11917378 :       akp = a(k1) + a(k4)
     303    11917378 :       akm = a(k1) - a(k4)
     304    11917378 :       bkp = b(k1) + b(k4)
     305    11917378 :       bkm = b(k1) - b(k4)
     306    11917378 :       ajp = a(k2) + a(k3)
     307    11917378 :       ajm = a(k2) - a(k3)
     308    11917378 :       bjp = b(k2) + b(k3)
     309    11917378 :       bjm = b(k2) - b(k3)
     310    11917378 :       aa = a(kk)
     311    11917378 :       bb = b(kk)
     312    11917378 :       a(kk) = aa + akp + ajp
     313    11917378 :       b(kk) = bb + bkp + bjp
     314    11917378 :       ak = akp*c72 + ajp*c2 + aa
     315    11917378 :       bk = bkp*c72 + bjp*c2 + bb
     316    11917378 :       aj = akm*s72 + ajm*s2
     317    11917378 :       bj = bkm*s72 + bjm*s2
     318    11917378 :       a(k1) = ak - bj
     319    11917378 :       a(k4) = ak + bj
     320    11917378 :       b(k1) = bk + aj
     321    11917378 :       b(k4) = bk - aj
     322    11917378 :       ak = akp*c2 + ajp*c72 + aa
     323    11917378 :       bk = bkp*c2 + bjp*c72 + bb
     324    11917378 :       aj = akm*s2 - ajm*s72
     325    11917378 :       bj = bkm*s2 - bjm*s72
     326    11917378 :       a(k2) = ak - bj
     327    11917378 :       a(k3) = ak + bj
     328    11917378 :       b(k2) = bk + aj
     329    11917378 :       b(k3) = bk - aj
     330    11917378 :       kk = k4 + kspan
     331    11917378 :       IF (kk < nn) GO TO 250
     332     3843730 :       kk = kk - nn
     333     3843730 :       IF (kk <= kspan) GO TO 250
     334      249717 :       GO TO 320
     335             : !     transform for odd factors
     336      249717 : 260   k = nfac(i)
     337      249717 :       kspnn = kspan
     338      249717 :       kspan = kspan/k
     339      249717 :       IF (k == 3) GO TO 160
     340       48209 :       IF (k == 5) GO TO 240
     341       35348 :       IF (k == jf) GO TO 280
     342       35348 :       jf = k
     343       35348 :       s1 = rad/real(k)
     344       35348 :       c1 = cos(s1)
     345       35348 :       s1 = sin(s1)
     346       35348 :       IF (jf > maxf) GO TO 590
     347       35348 :       ck(jf) = 1.0
     348       35348 :       sk(jf) = 0.0
     349      126942 :       j = 1
     350      126942 : 270   ck(j) = ck(k)*c1 + sk(k)*s1
     351      126942 :       sk(j) = ck(k)*s1 - sk(k)*c1
     352      126942 :       k = k - 1
     353      126942 :       ck(k) = ck(j)
     354      126942 :       sk(k) = -sk(j)
     355      126942 :       j = j + 1
     356      126942 :       IF (j < k) GO TO 270
     357    62940207 : 280   k1 = kk
     358    62940207 :       k2 = kk + kspnn
     359    62940207 :       aa = a(kk)
     360    62940207 :       bb = b(kk)
     361    62940207 :       ak = aa
     362    62940207 :       bk = bb
     363    62940207 :       j = 1
     364   217150281 :       k1 = k1 + kspan
     365   217150281 : 290   k2 = k2 - kspan
     366   217150281 :       j = j + 1
     367   217150281 :       at(j) = a(k1) + a(k2)
     368   217150281 :       ak = at(j) + ak
     369   217150281 :       bt(j) = b(k1) + b(k2)
     370   217150281 :       bk = bt(j) + bk
     371   217150281 :       j = j + 1
     372   217150281 :       at(j) = a(k1) - a(k2)
     373   217150281 :       bt(j) = b(k1) - b(k2)
     374   217150281 :       k1 = k1 + kspan
     375   217150281 :       IF (k1 < k2) GO TO 290
     376    62940207 :       a(kk) = ak
     377    62940207 :       b(kk) = bk
     378    62940207 :       k1 = kk
     379    62940207 :       k2 = kk + kspnn
     380   217150281 :       j = 1
     381   217150281 : 300   k1 = k1 + kspan
     382   217150281 :       k2 = k2 - kspan
     383   217150281 :       jj = j
     384   217150281 :       ak = aa
     385   217150281 :       bk = bb
     386   217150281 :       aj = 0.0
     387   217150281 :       bj = 0.0
     388   796669479 :       k = 1
     389   796669479 : 310   k = k + 1
     390   796669479 :       ak = at(k)*ck(jj) + ak
     391   796669479 :       bk = bt(k)*ck(jj) + bk
     392   796669479 :       k = k + 1
     393   796669479 :       aj = at(k)*sk(jj) + aj
     394   796669479 :       bj = bt(k)*sk(jj) + bj
     395   796669479 :       jj = jj + j
     396   796669479 :       IF (jj > jf) jj = jj - jf
     397   796669479 :       IF (k < jf) GO TO 310
     398   217150281 :       k = jf - j
     399   217150281 :       a(k1) = ak - bj
     400   217150281 :       b(k1) = bk + aj
     401   217150281 :       a(k2) = ak + bj
     402   217150281 :       b(k2) = bk - aj
     403   217150281 :       j = j + 1
     404   217150281 :       IF (j < k) GO TO 300
     405    62940207 :       kk = kk + kspnn
     406    62940207 :       IF (kk <= nn) GO TO 280
     407     5399445 :       kk = kk - nn
     408     5399445 :       IF (kk <= kspan) GO TO 280
     409             : !     multiply by rotation factor (except for factors of 2 and 4)
     410      249717 : 320   IF (i == m) GO TO 360
     411     9786723 :       kk = jc + 1
     412     9786723 : 330   c2 = 1.0 - cd
     413    72382479 :       s1 = sd
     414    72382479 : 340   c1 = c2
     415    72382479 :       s2 = s1
     416   149734830 :       kk = kk + kspan
     417   463515150 : 350   ak = a(kk)
     418   463515150 :       a(kk) = c2*ak - s2*b(kk)
     419   463515150 :       b(kk) = s2*ak + c2*b(kk)
     420   463515150 :       kk = kk + kspnn
     421   463515150 :       IF (kk <= nt) GO TO 350
     422   149734830 :       ak = s1*s2
     423   149734830 :       s2 = s1*c2 + c1*s2
     424   149734830 :       c2 = c1*c2 - ak
     425   149734830 :       kk = kk - nt + kspan
     426   149734830 :       IF (kk <= kspnn) GO TO 350
     427    72382479 :       c2 = c1 - (cd*c1 + sd*s1)
     428    72382479 :       s1 = s1 + (sd*c1 - cd*s1)
     429             : !     the following three statements compensate for truncation
     430             : !     error. if rounded arithmetic is used, they may
     431             : !     be deleted.
     432             : !     c1=0.5/(c2**2+s1**2)+0.5
     433             : !     s1=c1*s1
     434             : !     c2=c1*c2
     435    72382479 :       kk = kk - kspnn + jc
     436    72382479 :       IF (kk <= kspan) GO TO 340
     437     9786723 :       kk = kk - kspan + jc + inc
     438     9786723 :       IF (kk <= jc + jc) GO TO 330
     439      862943 :       GO TO 120
     440             : !     permute the results to normal order---done in two stages
     441             : !     permutation for square factors of n
     442      862943 : 360   np(1) = ks
     443      862943 :       IF (kt == 0) GO TO 450
     444      818327 :       k = kt + kt + 1
     445      818327 :       IF (m < k) k = k - 1
     446      818327 :       j = 1
     447      818361 :       np(k + 1) = jc
     448      818361 : 370   np(j + 1) = np(j)/nfac(j)
     449      818361 :       np(k) = np(k + 1)*nfac(j)
     450      818361 :       j = j + 1
     451      818361 :       k = k - 1
     452      818361 :       IF (j < k) GO TO 370
     453      818327 :       k3 = np(k + 1)
     454      818327 :       kspan = np(2)
     455      818327 :       kk = jc + 1
     456      818327 :       k2 = kspan + 1
     457      818327 :       j = 1
     458      818327 :       IF (n /= ntot) GO TO 410
     459             : !     permutation for single-variate transform (optional code)
     460       17190 : 380   ak = a(kk)
     461       17190 :       a(kk) = a(k2)
     462       17190 :       a(k2) = ak
     463       17190 :       bk = b(kk)
     464       17190 :       b(kk) = b(k2)
     465       17190 :       b(k2) = bk
     466       17190 :       kk = kk + inc
     467       17190 :       k2 = kspan + k2
     468       17190 :       IF (k2 < ks) GO TO 380
     469       18145 : 390   k2 = k2 - np(j)
     470       18145 :       j = j + 1
     471       18145 :       k2 = np(j + 1) + k2
     472       18145 :       IF (k2 > np(j)) GO TO 390
     473       26740 :       j = 1
     474       42975 : 400   IF (kk < k2) GO TO 380
     475       32470 :       kk = kk + inc
     476       32470 :       k2 = kspan + k2
     477       32470 :       IF (k2 < ks) GO TO 400
     478        5730 :       IF (kk < ks) GO TO 390
     479             :       jc = k3
     480    54624968 :       GO TO 450
     481             : !     permutation for multivariate transform
     482   153116006 : 410   k = kk + jc
     483   153116006 : 420   ak = a(kk)
     484   153116006 :       a(kk) = a(k2)
     485   153116006 :       a(k2) = ak
     486   153116006 :       bk = b(kk)
     487   153116006 :       b(kk) = b(k2)
     488   153116006 :       b(k2) = bk
     489   153116006 :       kk = kk + inc
     490   153116006 :       k2 = k2 + inc
     491   153116006 :       IF (kk < k) GO TO 420
     492    54624968 :       kk = kk + ks - jc
     493    54624968 :       k2 = k2 + ks - jc
     494    54624968 :       IF (kk < nt) GO TO 410
     495     4680428 :       k2 = k2 - nt + kspan
     496     4680428 :       kk = kk - nt + jc
     497     4680428 :       IF (k2 < ks) GO TO 410
     498     2938905 : 430   k2 = k2 - np(j)
     499     2938905 :       j = j + 1
     500     2938905 :       k2 = np(j + 1) + k2
     501     2938905 :       IF (k2 > np(j)) GO TO 430
     502     5874631 :       j = 1
     503     8659814 : 440   IF (kk < k2) GO TO 410
     504     6947094 :       kk = kk + jc
     505     6947094 :       k2 = kspan + k2
     506     6947094 :       IF (k2 < ks) GO TO 440
     507     1072463 :       IF (kk < ks) GO TO 430
     508       44616 :       jc = k3
     509      862943 : 450   IF (2*kt + 1 >= m) GO TO 667
     510       47461 :       kspnn = np(kt + 1)
     511             : !     permutation for square-free factors of n
     512       47461 :       j = m - kt
     513      104467 :       nfac(j + 1) = 1
     514      104467 : 460   nfac(j) = nfac(j)*nfac(j + 1)
     515      104467 :       j = j - 1
     516      104467 :       IF (j /= kt) GO TO 460
     517       47461 :       kt = kt + 1
     518       47461 :       nn = nfac(kt) - 1
     519       47461 :       IF (nn > maxp) GO TO 590
     520             :       jj = 0
     521             :       j = 0
     522      438864 :       GO TO 490
     523      438864 : 470   jj = jj - k2
     524      438864 :       k2 = kk
     525      438864 :       k = k + 1
     526     1567379 :       kk = nfac(k)
     527     1567379 : 480   jj = kk + jj
     528     1567379 :       IF (jj >= k2) GO TO 470
     529     1128515 :       np(j) = jj
     530     1175976 : 490   k2 = nfac(kt)
     531     1175976 :       k = kt + 1
     532     1175976 :       kk = nfac(k)
     533     1175976 :       j = j + 1
     534     1175976 :       IF (j <= nn) GO TO 480
     535             : !     determine the permutation cycles of length greater than 1
     536             :       j = 0
     537      836993 :       GO TO 510
     538      836993 : 500   k = kk
     539      836993 :       kk = np(k)
     540      836993 :       np(k) = -kk
     541      836993 :       IF (kk /= j) GO TO 500
     542      872477 :       k3 = kk
     543     1128515 : 510   j = j + 1
     544     1128515 :       kk = np(j)
     545     1128515 :       IF (kk < 0) GO TO 510
     546      291522 :       IF (kk /= j) GO TO 500
     547       82945 :       np(j) = -j
     548       82945 :       IF (j /= nn) GO TO 510
     549       47461 :       maxf = inc*maxf
     550             : !     reorder a and b, following the permutation cycles
     551    42429234 :       GO TO 580
     552    92839725 : 520   j = j - 1
     553    92839725 :       IF (np(j) < 0) GO TO 520
     554       84748 :       jj = jc
     555    42466521 : 530   kspan = jj
     556    42466521 :       IF (jj > maxf) kspan = maxf
     557    42466521 :       jj = jj - kspan
     558    42466521 :       k = np(j)
     559    42466521 :       kk = jc*k + ii + jj
     560    42466521 :       k1 = kk + kspan
     561   118941871 :       k2 = 0
     562   118941871 : 540   k2 = k2 + 1
     563   118941871 :       at(k2) = a(k1)
     564   118941871 :       bt(k2) = b(k1)
     565   118941871 :       k1 = k1 - inc
     566   118941871 :       IF (k1 /= kk) GO TO 540
     567   117669083 : 550   k1 = kk + kspan
     568   117669083 :       k2 = k1 - jc*(k + np(k))
     569   349166103 :       k = -np(k)
     570   349166103 : 560   a(k1) = a(k2)
     571   349166103 :       b(k1) = b(k2)
     572   349166103 :       k1 = k1 - inc
     573   349166103 :       k2 = k2 - inc
     574   349166103 :       IF (k1 /= kk) GO TO 560
     575   117669083 :       kk = k2
     576   117669083 :       IF (k /= j) GO TO 550
     577    42466521 :       k1 = kk + kspan
     578   118941871 :       k2 = 0
     579   118941871 : 570   k2 = k2 + 1
     580   118941871 :       a(k1) = at(k2)
     581   118941871 :       b(k1) = bt(k2)
     582   118941871 :       k1 = k1 - inc
     583   118941871 :       IF (k1 /= kk) GO TO 570
     584    42466521 :       IF (jj /= 0) GO TO 530
     585    42381773 :       IF (j /= 1) GO TO 520
     586     9196912 : 580   j = k3 + 1
     587     9196912 :       nt = nt - kspnn
     588     9196912 :       ii = nt - inc + 1
     589     9196912 :       IF (nt >= 0) GO TO 520
     590           0 :       GOTO 667
     591             : !     error finish, insufficient array storage
     592             : 590   CONTINUE
     593             : !     isn = 0
     594           0 :       WRITE (6, FMT=8000)
     595      815482 :       CALL juDFT_error('array bounds exceeded', calledby='cfft')
     596             : 8000  FORMAT('array bounds exceeded within subroutine cft')
     597             : 667   CONTINUE
     598      862943 :       DEALLOCATE (at, bt, ck, sk, nfac, np)
     599             :    END SUBROUTINE
     600             : 
     601             : #else
     602             :    SUBROUTINE cfft(a, b, ntot, n, nspan, isn)
     603             : 
     604             : !-------------------------------------------------------------*
     605             : ! driver routine for ccfft subroutine instead of cfft on cray *
     606             : !              and dcft, dcft2 and dcft3 essl routines on IBM *
     607             : !-------------------------------------------------------------*
     608             : 
     609             :       IMPLICIT NONE
     610             : 
     611             : ! ... input variables
     612             :       INTEGER :: ntot, n, nspan, isn
     613             :       REAL :: a(ntot), b(ntot)
     614             : 
     615             : ! ... local variables
     616             :       INTEGER :: i, ld1, ld2, n1, n2, n3, dimfft, idim, s(4)
     617             : 
     618             :       LOGICAL :: calc
     619             : 
     620             :       REAL, DIMENSION(:), ALLOCATABLE :: table, aux
     621             :       REAL, DIMENSION(:), ALLOCATABLE :: work, aux1, aux2
     622             :       COMPLEX, DIMENSION(:), ALLOCATABLE :: x
     623             : 
     624             :       INTEGER :: naux, naux1, naux2, lam, la1, la2
     625             :       REAL, PARAMETER :: scale = 1.0
     626             : ! ... save variables
     627             :       SAVE n1, n2, n3
     628             : 
     629             : ! ... data statements
     630             :       DATA s/1, 1, 1, 1/
     631             : 
     632             : ! ... now, what do we have to do ?
     633             : 
     634             :       IF ((ntot == n) .AND. (n == nspan)) THEN
     635             :          !  ...                                          1D-FFT
     636             :          dimfft = 1
     637             :          n1 = n
     638             :          n2 = 2
     639             :          n3 = 2
     640             :          calc = .TRUE.
     641             :       ELSE
     642             :          IF (n == nspan) THEN
     643             :             !  ...                                          2D or 3D first step
     644             :             n1 = n
     645             :             n2 = 0
     646             :             calc = .FALSE.
     647             :          ELSE
     648             :             IF (ntot == nspan) THEN
     649             :                !  ...                                          2D second step or 3D third step
     650             :                IF (n2 == 0) THEN
     651             :                   dimfft = 2
     652             :                   n2 = n
     653             :                   n3 = 1
     654             :                   calc = .TRUE.
     655             :                ELSE
     656             :                   dimfft = 3
     657             :                   n3 = n
     658             :                   calc = .TRUE.
     659             :                ENDIF
     660             :             ELSE
     661             :                !  ...                                          3D second step.
     662             :                n2 = n
     663             :                calc = .FALSE.
     664             :             ENDIF
     665             :          ENDIF
     666             :       ENDIF
     667             : 
     668             :       IF (calc) THEN
     669             : 
     670             :          ! ... build x from a and b
     671             : 
     672             :          ALLOCATE (x(ntot))
     673             :          x = (0.0, 0.0)
     674             :          DO i = 1, ntot
     675             :             x(i) = cmplx(a(i), b(i))
     676             :          ENDDO
     677             :          ! ... do the FFT
     678             : 
     679             : #ifdef CPP_AIX
     680             : 
     681             :          ld1 = n1
     682             :          ld2 = n1*n2
     683             :          IF (dimfft == 1) THEN
     684             :             naux1 = 20000
     685             :             IF (n1 > 2048) naux1 = naux1 + CEILING(2.28*n1)
     686             :             ALLOCATE (aux1(naux1), aux2(naux1))
     687             :          ELSEIF (dimfft == 2) THEN
     688             :             naux1 = 40000 + CEILING(2.28*(n1 + n2))
     689             :             IF (max(n1, n2) <= 2048) naux1 = 40000
     690             :             naux2 = 20000 + CEILING((2*max(n1, n2) + 256)* &
     691             :                                     (2.28 + min(64, n1, n2)))
     692             :             IF (max(n1, n2) < 256) naux2 = 20000
     693             :             ALLOCATE (aux1(naux1), aux2(naux2))
     694             :          ELSE IF (dimfft == 3) THEN
     695             :             IF (max(n2, n3) < 252) THEN
     696             :                IF (n1 <= 2048) THEN
     697             :                   naux = 60000
     698             :                ELSE
     699             :                   naux = 60000 + CEILING(4.56*n1)
     700             :                ENDIF
     701             :             ELSE
     702             :                la1 = CEILING((2*n2 + 256)*(min(64, n1) + 4.56))
     703             :                la2 = CEILING((2*n3 + 256)*(min(64, n1*n2) + 4.56))
     704             :                lam = max(la2, la1)
     705             :                IF ((n2 >= 252) .AND. (n3 < 252)) lam = la1
     706             :                IF ((n2 < 252) .AND. (n3 >= 252)) lam = la2
     707             :                IF (n1 <= 2048) THEN
     708             :                   naux = 60000 + lam
     709             :                ELSE
     710             :                   naux = 60000 + CEILING(4.56*n1) + lam
     711             :                ENDIF
     712             :             ENDIF
     713             :             ALLOCATE (aux(naux))
     714             :          ENDIF
     715             : #else
     716             : 
     717             :          ld1 = n1
     718             :          ld2 = n2
     719             :          s(1) = dimfft
     720             : 
     721             : #ifndef CPP_MPI
     722             :          ! t,j90:
     723             :          idim = 1024*n
     724             :          ALLOCATE (table(16*n + 100), work(idim))
     725             : #else
     726             :          ! t3e:
     727             :          idim = 2*n1*max(n2, 1)*max(n3, 1)
     728             :          ALLOCATE (table(12*(n1 + n2 + n3)), work(idim))
     729             : #endif
     730             : #endif
     731             :          IF (dimfft == 1) THEN
     732             : #ifdef CPP_AIX
     733             :             CALL dcft(-1, x, 1, 1, x, 1, 1, n, 1, -isn, 1.0, &
     734             :                       aux1, naux1, aux2, naux1)
     735             :             CALL dcft(0, x, 1, 1, x, 1, 1, n, 1, -isn, 1.0, &
     736             :                       aux1, naux1, aux2, naux1)
     737             : #else
     738             :             CALL ccfft(0, n, 1.0, x, x, table, work, s)
     739             :             CALL ccfft(isn, n, 1.0, x, x, table, work, s)
     740             : #endif
     741             :          ENDIF
     742             :          IF (dimfft == 2) THEN
     743             : #ifdef CPP_AIX
     744             :             CALL dcft2(-1, x, 1, n1, x, 1, n1, n1, n2, -isn, 1.0, &
     745             :                        aux1, naux1, aux2, naux2)
     746             :             CALL dcft2(0, x, 1, n1, x, 1, n1, n1, n2, -isn, 1.0, &
     747             :                        aux1, naux1, aux2, naux2)
     748             : #else
     749             :             CALL ccfft2d(0, n1, n2, 1.0, x, ld1, x, ld1, table, work, s)
     750             :             CALL ccfft2d(isn, n1, n2, 1.0, x, ld1, x, ld1, table, work, s)
     751             : #endif
     752             :          ENDIF
     753             :          IF (dimfft == 3) THEN
     754             : #ifdef CPP_AIX
     755             :             CALL dcft3(x, ld1, ld2, x, ld1, ld2, n1, n2, n3, &
     756             :                        -isn, scale, aux, naux)
     757             : #else
     758             :             CALL ccfft3d(0, n1, n2, n3, 1.0, x, ld1, ld2, x, ld1, ld2, &
     759             :                          table, work, s)
     760             :             CALL ccfft3d(isn, n1, n2, n3, 1.0, x, ld1, ld2, x, ld1, ld2, &
     761             :                          table, work, s)
     762             : #endif
     763             :          ENDIF
     764             : 
     765             : #ifdef CPP_AIX
     766             :          IF (dimfft == 3) THEN
     767             :             DEALLOCATE (aux)
     768             :          ELSE
     769             :             DEALLOCATE (aux1, aux2)
     770             :          ENDIF
     771             : #else
     772             :          DEALLOCATE (table, work)
     773             : #endif
     774             : 
     775             :          ! ... backup a and b
     776             : 
     777             :          DO i = 1, ntot
     778             :             a(i) = real(x(i))
     779             :             b(i) = aimag(x(i))
     780             :          ENDDO
     781             : 
     782             :          DEALLOCATE (x)
     783             :       ENDIF
     784             : 
     785             :    END subroutine
     786             : #endif
     787             : END module

Generated by: LCOV version 1.13