external/fftpack/fftpack5/mradbg.F

   1 subroutine mradbg ( m, ido, ip, l1, idl1, cc, c1, c2, im1, in1, &
   2   ch, ch2, im2, in2, wa )
   3
   4 !*****************************************************************************80
   5 !
   6 !! MRADBG 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 ) idl1
  38   integer ( kind = 4 ) ido
  39   integer ( kind = 4 ) in1
  40   integer ( kind = 4 ) in2
  41   integer ( kind = 4 ) ip
  42   integer ( kind = 4 ) l1
  43
  44   real ( kind = 4 ) ai1
  45   real ( kind = 4 ) ai2
  46   real ( kind = 4 ) ar1
  47   real ( kind = 4 ) ar1h
  48   real ( kind = 4 ) ar2
  49   real ( kind = 4 ) ar2h
  50   real ( kind = 4 ) arg
  51   real ( kind = 4 ) c1(in1,ido,l1,ip)
  52   real ( kind = 4 ) c2(in1,idl1,ip)
  53   real ( kind = 4 ) cc(in1,ido,ip,l1)
  54   real ( kind = 4 ) ch(in2,ido,l1,ip)
  55   real ( kind = 4 ) ch2(in2,idl1,ip)
  56   real ( kind = 4 ) dc2
  57   real ( kind = 4 ) dcp
  58   real ( kind = 4 ) ds2
  59   real ( kind = 4 ) dsp
  60   integer ( kind = 4 ) i
  61   integer ( kind = 4 ) ic
  62   integer ( kind = 4 ) idij
  63   integer ( kind = 4 ) idp2
  64   integer ( kind = 4 ) ik
  65   integer ( kind = 4 ) im1
  66   integer ( kind = 4 ) im2
  67   integer ( kind = 4 ) ipp2
  68   integer ( kind = 4 ) ipph
  69   integer ( kind = 4 ) is
  70   integer ( kind = 4 ) j
  71   integer ( kind = 4 ) j2
  72   integer ( kind = 4 ) jc
  73   integer ( kind = 4 ) k
  74   integer ( kind = 4 ) l
  75   integer ( kind = 4 ) lc
  76   integer ( kind = 4 ) m
  77   integer ( kind = 4 ) m1
  78   integer ( kind = 4 ) m1d
  79   integer ( kind = 4 ) m2
  80   integer ( kind = 4 ) m2s
  81   integer ( kind = 4 ) nbd
  82   real ( kind = 4 ) tpi
  83   real ( kind = 4 ) wa(ido)
  84
  85   m1d = ( m - 1 ) * im1 + 1
  86   m2s = 1 - im2
  87   tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 )
  88   arg = tpi / real ( ip, kind = 4 )
  89   dcp = cos ( arg )
  90   dsp = sin ( arg )
  91   idp2 = ido + 2
  92   nbd = ( ido - 1 ) / 2
  93   ipp2 = ip + 2
  94   ipph = ( ip + 1 ) / 2
  95
  96   if ( ido < l1 ) then
  97
  98     do i = 1, ido
  99       do k = 1, l1
 100         m2 = m2s
 101         do m1 = 1, m1d, im1
 102           m2 = m2 + im2
 103           ch(m2,i,k,1) = cc(m1,i,1,k)
 104         end do
 105       end do
 106     end do
 107
 108   else
 109
 110     do k = 1, l1
 111       do i = 1, ido
 112         m2 = m2s
 113         do m1 = 1, m1d, im1
 114           m2 = m2 + im2
 115           ch(m2,i,k,1) = cc(m1,i,1,k)
 116         end do
 117       end do
 118     end do
 119
 120   end if
 121
 122   do j = 2, ipph
 123     jc = ipp2 - j
 124     j2 = j + j
 125     do k = 1, l1
 126       m2 = m2s
 127       do m1 = 1, m1d, im1
 128         m2 = m2 + im2
 129         ch(m2,1,k,j) = cc(m1,ido,j2-2,k) + cc(m1,ido,j2-2,k)
 130         ch(m2,1,k,jc) = cc(m1,1,j2-1,k) + cc(m1,1,j2-1,k)
 131       end do
 132     end do
 133   end do
 134
 135   if ( ido == 1 ) then
 136
 137   else if ( nbd < l1 ) then
 138
 139     do j = 2, ipph
 140       jc = ipp2 - j
 141       do i = 3, ido, 2
 142         ic = idp2 - i
 143         do k = 1, l1
 144           m2 = m2s
 145           do m1 = 1, m1d, im1
 146             m2 = m2 + im2
 147             ch(m2,i-1,k,j)  = cc(m1,i-1,2*j-1,k) + cc(m1,ic-1,2*j-2,k)
 148             ch(m2,i-1,k,jc) = cc(m1,i-1,2*j-1,k) - cc(m1,ic-1,2*j-2,k)
 149             ch(m2,i,k,j)    = cc(m1,i,2*j-1,k)   - cc(m1,ic,2*j-2,k)
 150             ch(m2,i,k,jc)   = cc(m1,i,2*j-1,k)   + cc(m1,ic,2*j-2,k)
 151           end do
 152         end do
 153       end do
 154     end do
 155
 156   else
 157
 158     do j = 2, ipph
 159       jc = ipp2 - j
 160       do k = 1, l1
 161         do i = 3, ido, 2
 162           ic = idp2 - i
 163           m2 = m2s
 164           do m1 = 1, m1d, im1
 165             m2 = m2 + im2
 166             ch(m2,i-1,k,j)  = cc(m1,i-1,2*j-1,k) + cc(m1,ic-1,2*j-2,k)
 167             ch(m2,i-1,k,jc) = cc(m1,i-1,2*j-1,k) - cc(m1,ic-1,2*j-2,k)
 168             ch(m2,i,k,j)    = cc(m1,i,2*j-1,k)   - cc(m1,ic,2*j-2,k)
 169             ch(m2,i,k,jc)   = cc(m1,i,2*j-1,k)   + cc(m1,ic,2*j-2,k)
 170           end do
 171         end do
 172       end do
 173     end do
 174
 175   end if
 176
 177   ar1 = 1.0E+00
 178   ai1 = 0.0E+00
 179
 180   do l = 2, ipph
 181
 182     lc = ipp2 - l
 183     ar1h = dcp * ar1 - dsp * ai1
 184     ai1 =  dcp * ai1 + dsp * ar1
 185     ar1 = ar1h
 186
 187     do ik = 1, idl1
 188       m2 = m2s
 189       do m1 = 1, m1d, im1
 190         m2 = m2 + im2
 191         c2(m1,ik,l)  = ch2(m2,ik,1) + ar1 * ch2(m2,ik,2)
 192         c2(m1,ik,lc) =                ai1 * ch2(m2,ik,ip)
 193       end do
 194     end do
 195
 196     dc2 = ar1
 197     ds2 = ai1
 198     ar2 = ar1
 199     ai2 = ai1
 200
 201     do j = 3, ipph
 202       jc = ipp2 - j
 203       ar2h = dc2 * ar2 - ds2 * ai2
 204       ai2  = dc2 * ai2 + ds2 * ar2
 205       ar2 = ar2h
 206       do ik = 1, idl1
 207         m2 = m2s
 208         do m1 = 1, m1d, im1
 209           m2 = m2 + im2
 210           c2(m1,ik,l)  = c2(m1,ik,l)  + ar2 * ch2(m2,ik,j)
 211           c2(m1,ik,lc) = c2(m1,ik,lc) + ai2 * ch2(m2,ik,jc)
 212         end do
 213       end do
 214     end do
 215
 216   end do
 217
 218   do j = 2, ipph
 219     do ik = 1, idl1
 220       m2 = m2s
 221       do m1 = 1, m1d, im1
 222         m2 = m2 + im2
 223         ch2(m2,ik,1) = ch2(m2,ik,1) + ch2(m2,ik,j)
 224       end do
 225     end do
 226   end do
 227
 228   do j = 2, ipph
 229     jc = ipp2 - j
 230     do k = 1, l1
 231       m2 = m2s
 232       do m1 = 1, m1d, im1
 233         m2 = m2 + im2
 234         ch(m2,1,k,j)  = c1(m1,1,k,j) - c1(m1,1,k,jc)
 235         ch(m2,1,k,jc) = c1(m1,1,k,j) + c1(m1,1,k,jc)
 236       end do
 237     end do
 238   end do
 239
 240   if ( ido == 1 ) then
 241
 242   else if ( nbd < l1 ) then
 243
 244     do j = 2, ipph
 245       jc = ipp2 - j
 246       do i = 3, ido, 2
 247         do k = 1, l1
 248           m2 = m2s
 249           do m1 = 1, m1d, im1
 250             m2 = m2 + im2
 251             ch(m2,i-1,k,j)  = c1(m1,i-1,k,j) - c1(m1,i,k,jc)
 252             ch(m2,i-1,k,jc) = c1(m1,i-1,k,j) + c1(m1,i,k,jc)
 253             ch(m2,i,k,j)    = c1(m1,i,k,j)   + c1(m1,i-1,k,jc)
 254             ch(m2,i,k,jc)   = c1(m1,i,k,j)   - c1(m1,i-1,k,jc)
 255           end do
 256         end do
 257       end do
 258     end do
 259
 260   else
 261
 262     do j = 2, ipph
 263       jc = ipp2 - j
 264       do k = 1, l1
 265         do i = 3, ido, 2
 266           m2 = m2s
 267           do m1 = 1, m1d, im1
 268             m2 = m2 + im2
 269             ch(m2,i-1,k,j)  = c1(m1,i-1,k,j) - c1(m1,i,k,jc)
 270             ch(m2,i-1,k,jc) = c1(m1,i-1,k,j) + c1(m1,i,k,jc)
 271             ch(m2,i,k,j)    = c1(m1,i,k,j)   + c1(m1,i-1,k,jc)
 272             ch(m2,i,k,jc)   = c1(m1,i,k,j)   - c1(m1,i-1,k,jc)
 273           end do
 274         end do
 275       end do
 276     end do
 277
 278   end if
 279
 280   if ( ido == 1 ) then
 281     return
 282   end if
 283
 284   do ik = 1, idl1
 285     m2 = m2s
 286     do m1 = 1, m1d, im1
 287       m2 = m2 + im2
 288       c2(m1,ik,1) = ch2(m2,ik,1)
 289     end do
 290   end do
 291
 292   do j = 2, ip
 293     do k = 1, l1
 294       m2 = m2s
 295       do m1 = 1, m1d, im1
 296         m2 = m2 + im2
 297         c1(m1,1,k,j) = ch(m2,1,k,j)
 298       end do
 299     end do
 300   end do
 301
 302   if ( l1 < nbd ) then
 303
 304     is = -ido
 305
 306     do j = 2, ip
 307       is = is + ido
 308       do k = 1, l1
 309         idij = is
 310         do i = 3, ido, 2
 311           idij = idij + 2
 312           m2 = m2s
 313           do m1 = 1, m1d, im1
 314             m2 = m2 + im2
 315             c1(m1,i-1,k,j) = wa(idij-1) * ch(m2,i-1,k,j) &
 316                            - wa(idij)   * ch(m2,i,k,j)
 317             c1(m1,i,k,j) =   wa(idij-1) * ch(m2,i,k,j) &
 318                            + wa(idij)   * ch(m2,i-1,k,j)
 319           end do
 320         end do
 321       end do
 322     end do
 323
 324   else
 325
 326     is = -ido
 327
 328     do j = 2, ip
 329       is = is + ido
 330       idij = is
 331       do i = 3, ido, 2
 332         idij = idij + 2
 333         do k = 1, l1
 334           m2 = m2s
 335           do m1 = 1, m1d, im1
 336             m2 = m2 + im2
 337             c1(m1,i-1,k,j) = wa(idij-1) * ch(m2,i-1,k,j) &
 338                            - wa(idij)   * ch(m2,i,k,j)
 339             c1(m1,i,k,j) =   wa(idij-1) * ch(m2,i,k,j) &
 340                            + wa(idij)   * ch(m2,i-1,k,j)
 341           end do
 342         end do
 343       end do
 344     end do
 345
 346   end if
 347
 348   return
 349 end