external/fftpack/fftpack5/mcstf1.F

   1 subroutine mcstf1 ( lot, jump, n, inc, x, wsave, dsum, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! MCSTF1 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   implicit none
  35
  36   integer ( kind = 4 ) inc
  37
  38   real ( kind = 8 ) dsum(*)
  39   integer ( kind = 4 ) i
  40   integer ( kind = 4 ) ier
  41   integer ( kind = 4 ) ier1
  42   integer ( kind = 4 ) jump
  43   integer ( kind = 4 ) k
  44   integer ( kind = 4 ) kc
  45   integer ( kind = 4 ) lenx
  46   integer ( kind = 4 ) lj
  47   integer ( kind = 4 ) lnsv
  48   integer ( kind = 4 ) lnwk
  49   integer ( kind = 4 ) lot
  50   integer ( kind = 4 ) m
  51   integer ( kind = 4 ) m1
  52   integer ( kind = 4 ) modn
  53   integer ( kind = 4 ) n
  54   integer ( kind = 4 ) nm1
  55   integer ( kind = 4 ) np1
  56   integer ( kind = 4 ) ns2
  57   real ( kind = 4 ) snm1
  58   real ( kind = 4 ) t1
  59   real ( kind = 4 ) t2
  60   real ( kind = 4 ) tx2
  61   real ( kind = 4 ) work(*)
  62   real ( kind = 4 ) wsave(*)
  63   real ( kind = 4 ) x(inc,*)
  64   real ( kind = 4 ) x1h
  65   real ( kind = 4 ) x1p3
  66   real ( kind = 4 ) xi
  67
  68   ier = 0
  69   nm1 = n - 1
  70   np1 = n + 1
  71   ns2 = n / 2
  72   lj = ( lot - 1 ) * jump + 1
  73
  74   if ( n < 2 ) then
  75     return
  76   end if
  77
  78   if ( n == 2 ) then
  79
  80     do m = 1, lj, jump
  81       x1h = x(m,1) + x(m,2)
  82       x(m,2) = 0.5E+00 * ( x(m,1) - x(m,2) )
  83       x(m,1) = 0.5E+00 * x1h
  84     end do
  85
  86     return
  87
  88   end if
  89
  90   if ( n == 3 ) then
  91
  92     do m = 1, lj, jump
  93       x1p3 = x(m,1) + x(m,3)
  94       tx2 = x(m,2) + x(m,2)
  95       x(m,2) = 0.5E+00 * ( x(m,1) - x(m,3) )
  96       x(m,1) = 0.25E+00 * ( x1p3 + tx2 )
  97       x(m,3) = 0.25E+00 * ( x1p3 - tx2 )
  98     end do
  99
 100     return
 101   end if
 102
 103   m1 = 0
 104   do m = 1, lj, jump
 105     m1 = m1 + 1
 106     dsum(m1) = x(m,1) - x(m,n)
 107     x(m,1) = x(m,1) + x(m,n)
 108   end do
 109
 110   do k = 2, ns2
 111     m1 = 0
 112     do m = 1, lj, jump
 113       m1 = m1 + 1
 114       kc = np1 - k
 115       t1 = x(m,k) + x(m,kc)
 116       t2 = x(m,k) - x(m,kc)
 117       dsum(m1) = dsum(m1) + wsave(kc) * t2
 118       t2 = wsave(k) * t2
 119       x(m,k) = t1 - t2
 120       x(m,kc) = t1 + t2
 121     end do
 122   end do
 123
 124   modn = mod ( n, 2 )
 125
 126   if ( modn /= 0 ) then
 127     do m = 1, lj, jump
 128       x(m,ns2+1) = x(m,ns2+1) + x(m,ns2+1)
 129     end do
 130   end if
 131
 132   lenx = ( lot - 1 ) * jump + inc * ( nm1 - 1 )  + 1
 133   lnsv = nm1 + int ( log ( real ( nm1, kind = 4 ) ) ) + 4
 134   lnwk = lot * nm1
 135
 136   call rfftmf ( lot, jump, nm1, inc, x, lenx, wsave(n+1), lnsv, work, &
 137     lnwk, ier1 )
 138
 139   if ( ier1 /= 0 ) then
 140     ier = 20
 141     call xerfft ( 'mcstf1', -5 )
 142     return
 143   end if
 144
 145   snm1 = 1.0E+00 / real ( nm1, kind = 4 )
 146   do m = 1, lot
 147     dsum(m) = snm1 * dsum(m)
 148   end do
 149
 150   if ( mod ( nm1, 2 ) == 0 ) then
 151     do m = 1, lj, jump
 152       x(m,nm1) = x(m,nm1) + x(m,nm1)
 153     end do
 154   end if
 155
 156   do i = 3, n, 2
 157     m1 = 0
 158     do m = 1, lj, jump
 159       m1 = m1 + 1
 160       xi = 0.5E+00 * x(m,i)
 161       x(m,i) = 0.5E+00 * x(m,i-1)
 162       x(m,i-1) = dsum(m1)
 163       dsum(m1) = dsum(m1) + xi
 164     end do
 165   end do
 166
 167   if ( modn == 0 ) then
 168
 169     m1 = 0
 170     do m = 1, lj, jump
 171       m1 = m1 + 1
 172       x(m,n) = dsum(m1)
 173     end do
 174
 175   end if
 176
 177   do m = 1, lj, jump
 178     x(m,1) = 0.5E+00 * x(m,1)
 179     x(m,n) = 0.5E+00 * x(m,n)
 180   end do
 181
 182   return
 183 end