external/fftpack/fftpack5/dffti1.F

   1 subroutine dffti1 ( n, wa, fac )
   2
   3 !*****************************************************************************80
   4 !
   5 !! DFFTI1 is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    07 February 2006
  12 !
  13 !  Author:
  14 !
  15 !    Original real single precision by Paul Swarztrauber, Richard Valent.
  16 !    Real double precision version by John Burkardt.
  17 !
  18 !  Reference:
  19 !
  20 !    Paul Swarztrauber,
  21 !    Vectorizing the Fast Fourier Transforms,
  22 !    in Parallel Computations,
  23 !    edited by G. Rodrigue,
  24 !    Academic Press, 1982.
  25 !
  26 !    Paul Swarztrauber,
  27 !    Fast Fourier Transform Algorithms for Vector Computers,
  28 !    Parallel Computing, pages 45-63, 1984.
  29 !
  30 !  Parameters:
  31 !
  32 !    Input, integer ( kind = 4 ) N, the number for which factorization and
  33 !    other information is needed.
  34 !
  35 !    Output, real ( kind = 8 ) WA(N), trigonometric information.
  36 !
  37 !    Output, real ( kind = 8 ) FAC(15), factorization information.
  38 !    FAC(1) is N, FAC(2) is NF, the number of factors, and FAC(3:NF+2) are the
  39 !    factors.
  40 !
  41   implicit none
  42
  43   integer ( kind = 4 ) n
  44
  45   real ( kind = 8 ) arg
  46   real ( kind = 8 ) argh
  47   real ( kind = 8 ) argld
  48   real ( kind = 8 ) fac(15)
  49   real ( kind = 8 ) fi
  50   integer ( kind = 4 ) i
  51   integer ( kind = 4 ) ib
  52   integer ( kind = 4 ) ido
  53   integer ( kind = 4 ) ii
  54   integer ( kind = 4 ) ip
  55   integer ( kind = 4 ) ipm
  56   integer ( kind = 4 ) is
  57   integer ( kind = 4 ) j
  58   integer ( kind = 4 ) k1
  59   integer ( kind = 4 ) l1
  60   integer ( kind = 4 ) l2
  61   integer ( kind = 4 ) ld
  62   integer ( kind = 4 ) nf
  63   integer ( kind = 4 ) nfm1
  64   integer ( kind = 4 ) nl
  65   integer ( kind = 4 ) nq
  66   integer ( kind = 4 ) nr
  67   integer ( kind = 4 ) ntry
  68   real ( kind = 8 ) tpi
  69   real ( kind = 8 ) wa(n)
  70
  71   nl = n
  72   nf = 0
  73   j = 0
  74
  75   do while ( 1 < nl )
  76
  77     j = j + 1
  78
  79     if ( j == 1 ) then
  80       ntry = 4
  81     else if ( j == 2 ) then
  82       ntry = 2
  83     else if ( j == 3 ) then
  84       ntry = 3
  85     else if ( j == 4 ) then
  86       ntry = 5
  87     else
  88       ntry = ntry + 2
  89     end if
  90
  91     do
  92
  93       nq = nl / ntry
  94       nr = nl - ntry * nq
  95
  96       if ( nr /= 0 ) then
  97         exit
  98       end if
  99
 100       nf = nf + 1
 101       fac(nf+2) = real ( ntry, kind = 8 )
 102       nl = nq
 103 !
 104 !  If 2 is a factor, make sure it appears first in the list of factors.
 105 !
 106       if ( ntry == 2 ) then
 107         if ( nf /= 1 ) then
 108           do i = 2, nf
 109             ib = nf - i + 2
 110             fac(ib+2) = fac(ib+1)
 111           end do
 112           fac(3) = 2.0D+00
 113         end if
 114       end if
 115
 116     end do
 117
 118   end do
 119
 120   fac(1) = real ( n, kind = 8 )
 121   fac(2) = real ( nf, kind = 8 )
 122   tpi = 8.0D+00 * atan ( 1.0D+00 )
 123   argh = tpi / real ( n, kind = 8 )
 124   is = 0
 125   nfm1 = nf - 1
 126   l1 = 1
 127
 128   do k1 = 1, nfm1
 129     ip = int ( fac(k1+2) )
 130     ld = 0
 131     l2 = l1 * ip
 132     ido = n / l2
 133     ipm = ip - 1
 134     do j = 1, ipm
 135       ld = ld + l1
 136       i = is
 137       argld = real ( ld, kind = 8 ) * argh
 138       fi = 0.0D+00
 139       do ii = 3, ido, 2
 140         i = i + 2
 141         fi = fi + 1.0D+00
 142         arg = fi * argld
 143         wa(i-1) = real ( cos ( arg ), kind = 8 )
 144         wa(i) = real ( sin ( arg ), kind = 8 )
 145       end do
 146       is = is + ido
 147     end do
 148     l1 = l2
 149   end do
 150
 151   return
 152 end