external/fftpack/fftpack5/mradb5.F

   1 subroutine mradb5 ( m, ido, l1, cc, im1, in1, ch, im2, in2, &
   2   wa1, wa2, wa3, wa4 )
   3
   4 !*****************************************************************************80
   5 !
   6 !! MRADB5 is an FFTPACK5 auxiliary routine.
   7 !
   8 !
   9 !    Copyright (C) 1995-2004, Scientific Computing Division,
  10 !    University Corporation for Atmospheric Research
  11 !
  12 !  Modified:
  13 !
  14 !    27 March 2009
  15 !
  16 !  Author:
  17 !
  18 !    Paul Swarztrauber
  19 !    Richard Valent
  20 !
  21 !  Reference:
  22 !
  23 !    Paul Swarztrauber,
  24 !    Vectorizing the Fast Fourier Transforms,
  25 !    in Parallel Computations,
  26 !    edited by G. Rodrigue,
  27 !    Academic Press, 1982.
  28 !
  29 !    Paul Swarztrauber,
  30 !    Fast Fourier Transform Algorithms for Vector Computers,
  31 !    Parallel Computing, pages 45-63, 1984.
  32 !
  33 !  Parameters:
  34 !
  35   implicit none
  36
  37   integer ( kind = 4 ) ido
  38   integer ( kind = 4 ) in1
  39   integer ( kind = 4 ) in2
  40   integer ( kind = 4 ) l1
  41
  42   real ( kind = 4 ) arg
  43   real ( kind = 4 ) cc(in1,ido,5,l1)
  44   real ( kind = 4 ) ch(in2,ido,l1,5)
  45   integer ( kind = 4 ) i
  46   integer ( kind = 4 ) ic
  47   integer ( kind = 4 ) idp2
  48   integer ( kind = 4 ) im1
  49   integer ( kind = 4 ) im2
  50   integer ( kind = 4 ) k
  51   integer ( kind = 4 ) m
  52   integer ( kind = 4 ) m1
  53   integer ( kind = 4 ) m1d
  54   integer ( kind = 4 ) m2
  55   integer ( kind = 4 ) m2s
  56   real ( kind = 4 ) ti11
  57   real ( kind = 4 ) ti12
  58   real ( kind = 4 ) tr11
  59   real ( kind = 4 ) tr12
  60   real ( kind = 4 ) wa1(ido)
  61   real ( kind = 4 ) wa2(ido)
  62   real ( kind = 4 ) wa3(ido)
  63   real ( kind = 4 ) wa4(ido)
  64
  65   m1d = ( m - 1 ) * im1 + 1
  66   m2s = 1 - im2
  67   arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 5.E+00
  68   tr11 = cos ( arg )
  69   ti11 = sin ( arg )
  70   tr12 = cos ( 2.0E+00 * arg )
  71   ti12 = sin ( 2.0E+00 * arg )
  72
  73   do k = 1, l1
  74     m2 = m2s
  75     do m1 = 1, m1d, im1
  76       m2 = m2 + im2
  77       ch(m2,1,k,1) = cc(m1,1,1,k) + 2.0E+00 * cc(m1,ido,2,k) &
  78         + 2.0E+00 * cc(m1,ido,4,k)
  79       ch(m2,1,k,2) = ( cc(m1,1,1,k) + tr11 * 2.0E+00 * cc(m1,ido,2,k) &
  80         + tr12 * 2.0E+00 * cc(m1,ido,4,k) ) - ( ti11 * 2.0E+00 * cc(m1,1,3,k) &
  81         + ti12 * 2.0E+00 * cc(m1,1,5,k) )
  82       ch(m2,1,k,3) = ( cc(m1,1,1,k) + tr12 * 2.0E+00 * cc(m1,ido,2,k) &
  83         + tr11 * 2.0E+00 * cc(m1,ido,4,k) ) - ( ti12 * 2.0E+00 * cc(m1,1,3,k) &
  84         - ti11 * 2.0E+00 * cc(m1,1,5,k) )
  85       ch(m2,1,k,4) = ( cc(m1,1,1,k) + tr12 * 2.0E+00 * cc(m1,ido,2,k) &
  86         + tr11 * 2.0E+00 * cc(m1,ido,4,k) ) + ( ti12 * 2.0E+00 * cc(m1,1,3,k) &
  87         - ti11 * 2.0E+00 * cc(m1,1,5,k) )
  88       ch(m2,1,k,5) = ( cc(m1,1,1,k) + tr11 * 2.0E+00 * cc(m1,ido,2,k) &
  89         + tr12 * 2.0E+00 * cc(m1,ido,4,k) ) + ( ti11 * 2.0E+00 * cc(m1,1,3,k) &
  90         + ti12 * 2.0E+00 * cc(m1,1,5,k) )
  91     end do
  92   end do
  93
  94   if ( ido == 1 ) then
  95     return
  96   end if
  97
  98   idp2 = ido + 2
  99   do k = 1, l1
 100     do i = 3, ido, 2
 101       ic = idp2 - i
 102       m2 = m2s
 103       do m1 = 1, m1d, im1
 104         m2 = m2 + im2
 105         ch(m2,i-1,k,1) = cc(m1,i-1,1,k)+(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 106           +(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k))
 107         ch(m2,i,k,1) = cc(m1,i,1,k)+(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
 108           +(cc(m1,i,5,k)-cc(m1,ic,4,k))
 109         ch(m2,i-1,k,2) = wa1(i-2)*((cc(m1,i-1,1,k)+tr11* &
 110           (cc(m1,i-1,3,k)+cc(m1,ic-1,2,k))+tr12 &
 111           *(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti11*(cc(m1,i,3,k) &
 112           +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
 113           -wa1(i-1)*((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
 114           +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))+(ti11*(cc(m1,i-1,3,k) &
 115           -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
 116         ch(m2,i,k,2) = wa1(i-2)*((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k) &
 117           -cc(m1,ic,2,k))+tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
 118           +(ti11*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))+ti12 &
 119           *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) + wa1(i-1) &
 120           *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k) &
 121           +cc(m1,ic-1,2,k))+tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k))) &
 122           -(ti11*(cc(m1,i,3,k)+cc(m1,ic,2,k))+ti12 &
 123           *(cc(m1,i,5,k)+cc(m1,ic,4,k))))
 124         ch(m2,i-1,k,3) = wa2(i-2) &
 125           *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 126           +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti12*(cc(m1,i,3,k) &
 127           +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
 128           -wa2(i-1) &
 129           *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
 130         cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
 131           +(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
 132           *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
 133         ch(m2,i,k,3) = wa2(i-2) &
 134           *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
 135           cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
 136           +(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
 137           *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
 138           + wa2(i-1) &
 139           *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 140           +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti12*(cc(m1,i,3,k) &
 141           +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
 142         ch(m2,i-1,k,4) = wa3(i-2) &
 143           *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 144           +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti12*(cc(m1,i,3,k) &
 145           +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
 146           -wa3(i-1) &
 147           *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
 148           cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
 149           -(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
 150           *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
 151         ch(m2,i,k,4) = wa3(i-2) &
 152           *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
 153           cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
 154           -(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
 155           *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
 156           + wa3(i-1) &
 157           *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 158           +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti12*(cc(m1,i,3,k) &
 159           +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
 160         ch(m2,i-1,k,5) = wa4(i-2) &
 161           *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 162           +tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti11*(cc(m1,i,3,k) &
 163           +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
 164           -wa4(i-1) &
 165           *((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
 166           +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))-(ti11*(cc(m1,i-1,3,k) &
 167           -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
 168         ch(m2,i,k,5) = wa4(i-2) &
 169           *((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
 170           +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))-(ti11*(cc(m1,i-1,3,k) &
 171           -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
 172           + wa4(i-1) &
 173           *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
 174           +tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti11*(cc(m1,i,3,k) &
 175           +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
 176       end do
 177     end do
 178   end do
 179
 180   return
 181 end