external/fftpack/fftpack5/mcsqb1.F

   1 subroutine mcsqb1 ( lot, jump, n, inc, x, wsave, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! MCSQB1 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   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 ) m
  50   integer ( kind = 4 ) m1
  51   integer ( kind = 4 ) modn
  52   integer ( kind = 4 ) n
  53   integer ( kind = 4 ) np2
  54   integer ( kind = 4 ) ns2
  55   real ( kind = 4 ) work(lot,*)
  56   real ( kind = 4 ) wsave(*)
  57   real ( kind = 4 ) x(inc,*)
  58   real ( kind = 4 ) xim1
  59
  60   ier = 0
  61   lj = ( lot - 1 ) * jump + 1
  62   ns2 = ( n + 1 ) / 2
  63   np2 = n + 2
  64
  65   do i = 3, n, 2
  66     do m = 1, lj, jump
  67       xim1 = x(m,i-1) + x(m,i)
  68       x(m,i) = 0.5E+00 * ( x(m,i-1) - x(m,i) )
  69       x(m,i-1) = 0.5E+00 * xim1
  70     end do
  71   end do
  72
  73   do m = 1, lj, jump
  74     x(m,1) = 0.5E+00 * x(m,1)
  75   end do
  76
  77   modn = mod ( n, 2 )
  78   if ( modn == 0 ) then
  79     do m = 1, lj, jump
  80       x(m,n) = 0.5E+00 * x(m,n)
  81     end do
  82   end if
  83
  84   lenx = ( lot - 1 ) * jump + inc * ( n - 1 )  + 1
  85   lnsv = n + int ( log ( real ( n, kind = 4 ) ) ) + 4
  86   lnwk = lot * n
  87
  88   call rfftmb ( lot, jump, n, inc, x, lenx, wsave(n+1), lnsv, &
  89     work, lnwk, ier1 )
  90
  91   if ( ier1 /= 0 ) then
  92     ier = 20
  93     call xerfft ( 'mcsqb1', -5 )
  94     return
  95   end if
  96
  97   do k = 2, ns2
  98     kc = np2 - k
  99     m1 = 0
 100     do m = 1, lj, jump
 101       m1 = m1 + 1
 102       work(m1,k) = wsave(k-1) * x(m,kc) + wsave(kc-1) * x(m,k)
 103       work(m1,kc) = wsave(k-1) * x(m,k) - wsave(kc-1) * x(m,kc)
 104     end do
 105   end do
 106
 107   if ( modn == 0 ) then
 108     do m = 1, lj, jump
 109       x(m,ns2+1) = wsave(ns2) * ( x(m,ns2+1) + x(m,ns2+1) )
 110     end do
 111   end if
 112
 113   do k = 2, ns2
 114     kc = np2 - k
 115     m1 = 0
 116     do m = 1, lj, jump
 117       m1 = m1 + 1
 118       x(m,k) = work(m1,k) + work(m1,kc)
 119       x(m,kc) = work(m1,k) - work(m1,kc)
 120     end do
 121   end do
 122
 123   do m = 1, lj, jump
 124     x(m,1) = x(m,1) + x(m,1)
 125   end do
 126
 127   return
 128 end