external/fftpack/fftpack5/msntb1.F

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