external/fftpack/fftpack5/mrfti1.F

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