website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/rfftb1.f

   1 *DECK RFFTB1
   2       SUBROUTINE RFFTB1 (N, C, CH, WA, IFAC)
   3 C***BEGIN PROLOGUE  RFFTB1
   4 C***PURPOSE  Compute the backward fast Fourier transform of a real
   5 C            coefficient array.
   6 C***LIBRARY   SLATEC (FFTPACK)
   7 C***CATEGORY  J1A1
   8 C***TYPE      SINGLE PRECISION (RFFTB1-S, CFFTB1-C)
   9 C***KEYWORDS  FFTPACK, FOURIER TRANSFORM
  10 C***AUTHOR  Swarztrauber, P. N., (NCAR)
  11 C***DESCRIPTION
  12 C
  13 C   Subroutine RFFTB1 computes the real periodic sequence from its
  14 C   Fourier coefficients (Fourier synthesis).  The transform is defined
  15 C   below at output parameter C.
  16 C
  17 C   The arrays WA and IFAC which are used by subroutine RFFTB1 must be
  18 C   initialized by calling subroutine RFFTI1.
  19 C
  20 C   Input Arguments
  21 C
  22 C   N       the length of the array R to be transformed.  The method
  23 C           is most efficient when N is a product of small primes.
  24 C           N may change so long as different work arrays are provided.
  25 C
  26 C   C       a real array of length N which contains the sequence
  27 C           to be transformed.
  28 C
  29 C   CH      a real work array of length at least N.
  30 C
  31 C   WA      a real work array which must be dimensioned at least N.
  32 C
  33 C   IFAC    an integer work array which must be dimensioned at least 15.
  34 C
  35 C           The WA and IFAC arrays must be initialized by calling
  36 C           subroutine RFFTI1, and different WA and IFAC arrays must be
  37 C           used for each different value of N.  This initialization
  38 C           does not have to be repeated so long as N remains unchanged.
  39 C           Thus subsequent transforms can be obtained faster than the
  40 C           first.  The same WA and IFAC arrays can be used by RFFTF1
  41 C           and RFFTB1.
  42 C
  43 C   Output Argument
  44 C
  45 C   C       For N even and for I = 1,...,N
  46 C
  47 C                C(I) = C(1)+(-1)**(I-1)*C(N)
  48 C
  49 C                     plus the sum from K=2 to K=N/2 of
  50 C
  51 C                      2.*C(2*K-2)*COS((K-1)*(I-1)*2*PI/N)
  52 C
  53 C                     -2.*C(2*K-1)*SIN((K-1)*(I-1)*2*PI/N)
  54 C
  55 C           For N odd and for I = 1,...,N
  56 C
  57 C                C(I) = C(1) plus the sum from K=2 to K=(N+1)/2 of
  58 C
  59 C                     2.*C(2*K-2)*COS((K-1)*(I-1)*2*PI/N)
  60 C
  61 C                    -2.*C(2*K-1)*SIN((K-1)*(I-1)*2*PI/N)
  62 C
  63 C   Notes:  This transform is unnormalized since a call of RFFTF1
  64 C           followed by a call of RFFTB1 will multiply the input
  65 C           sequence by N.
  66 C
  67 C           WA and IFAC contain initialization calculations which must
  68 C           not be destroyed between calls of subroutine RFFTF1 or
  69 C           RFFTB1.
  70 C
  71 C***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
  72 C                 Computations (G. Rodrigue, ed.), Academic Press,
  73 C                 1982, pp. 51-83.
  74 C***ROUTINES CALLED  RADB2, RADB3, RADB4, RADB5, RADBG
  75 C***REVISION HISTORY  (YYMMDD)
  76 C   790601  DATE WRITTEN
  77 C   830401  Modified to use SLATEC library source file format.
  78 C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
  79 C           changing dummy array size declarations (1) to (*).
  80 C   881128  Modified by Dick Valent to meet prologue standards.
  81 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  82 C   900131  Routine changed from subsidiary to user-callable.  (WRB)
  83 C   920501  Reformatted the REFERENCES section.  (WRB)
  84 C***END PROLOGUE  RFFTB1
  85       DIMENSION CH(*), C(*), WA(*), IFAC(*)
  86 C***FIRST EXECUTABLE STATEMENT  RFFTB1
  87       NF = IFAC(2)
  88       NA = 0
  89       L1 = 1
  90       IW = 1
  91       DO 116 K1=1,NF
  92          IP = IFAC(K1+2)
  93          L2 = IP*L1
  94          IDO = N/L2
  95          IDL1 = IDO*L1
  96          IF (IP .NE. 4) GO TO 103
  97          IX2 = IW+IDO
  98          IX3 = IX2+IDO
  99          IF (NA .NE. 0) GO TO 101
 100          CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
 101          GO TO 102
 102   101    CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
 103   102    NA = 1-NA
 104          GO TO 115
 105   103    IF (IP .NE. 2) GO TO 106
 106          IF (NA .NE. 0) GO TO 104
 107          CALL RADB2 (IDO,L1,C,CH,WA(IW))
 108          GO TO 105
 109   104    CALL RADB2 (IDO,L1,CH,C,WA(IW))
 110   105    NA = 1-NA
 111          GO TO 115
 112   106    IF (IP .NE. 3) GO TO 109
 113          IX2 = IW+IDO
 114          IF (NA .NE. 0) GO TO 107
 115          CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
 116          GO TO 108
 117   107    CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
 118   108    NA = 1-NA
 119          GO TO 115
 120   109    IF (IP .NE. 5) GO TO 112
 121          IX2 = IW+IDO
 122          IX3 = IX2+IDO
 123          IX4 = IX3+IDO
 124          IF (NA .NE. 0) GO TO 110
 125          CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
 126          GO TO 111
 127   110    CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
 128   111    NA = 1-NA
 129          GO TO 115
 130   112    IF (NA .NE. 0) GO TO 113
 131          CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
 132          GO TO 114
 133   113    CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
 134   114    IF (IDO .EQ. 1) NA = 1-NA
 135   115    L1 = L2
 136          IW = IW+(IP-1)*IDO
 137   116 CONTINUE
 138       IF (NA .EQ. 0) RETURN
 139       DO 117 I=1,N
 140          C(I) = CH(I)
 141   117 CONTINUE
 142       RETURN
 143       END