share/fftpack5/fortran/mcstb1.f

   1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   2 C
   3 C   FFTPACK 5.0
   4 C
   5 C   Authors:  Paul N. Swarztrauber and Richard A. Valent
   6 C
   7 C   $Id$
   8 C
   9 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  10
  11       SUBROUTINE MCSTB1(LOT,JUMP,N,INC,X,WSAVE,DSUM,WORK,IER)
  12       REAL       X(INC,*)       ,WSAVE(*)
  13       DOUBLE PRECISION           DSUM(*)
  14       IER = 0
  15       NM1 = N-1
  16       NP1 = N+1
  17       NS2 = N/2
  18       LJ = (LOT-1)*JUMP+1
  19       IF (N-2) 106,101,102
  20   101 DO 111 M=1,LJ,JUMP
  21          X1H = X(M,1)+X(M,2)
  22          X(M,2) = X(M,1)-X(M,2)
  23          X(M,1) = X1H
  24   111 CONTINUE
  25       RETURN
  26   102 IF (N .GT. 3) GO TO 103
  27       DO 112 M=1,LJ,JUMP
  28          X1P3 = X(M,1)+X(M,3)
  29          X2 = X(M,2)
  30          X(M,2) = X(M,1)-X(M,3)
  31          X(M,1) = X1P3+X2
  32          X(M,3) = X1P3-X2
  33   112 CONTINUE
  34       RETURN
  35  103  DO 118 M=1,LJ,JUMP
  36       X(M,1) = X(M,1)+X(M,1)
  37       X(M,N) = X(M,N)+X(M,N)
  38  118  CONTINUE
  39       M1 = 0
  40       DO 113 M=1,LJ,JUMP
  41          M1 = M1+1
  42          DSUM(M1) = X(M,1)-X(M,N)
  43          X(M,1) = X(M,1)+X(M,N)
  44   113 CONTINUE
  45       DO 104 K=2,NS2
  46          M1 = 0
  47          DO 114 M=1,LJ,JUMP
  48          M1 = M1+1
  49          KC = NP1-K
  50          T1 = X(M,K)+X(M,KC)
  51          T2 = X(M,K)-X(M,KC)
  52          DSUM(M1) = DSUM(M1)+WSAVE(KC)*T2
  53          T2 = WSAVE(K)*T2
  54          X(M,K) = T1-T2
  55          X(M,KC) = T1+T2
  56   114    CONTINUE
  57   104 CONTINUE
  58       MODN = MOD(N,2)
  59       IF (MODN .EQ. 0) GO TO 124
  60          DO 123 M=1,LJ,JUMP
  61          X(M,NS2+1) = X(M,NS2+1)+X(M,NS2+1)
  62   123    CONTINUE
  63  124  CONTINUE
  64       LENX = (LOT-1)*JUMP + INC*(NM1-1)  + 1
  65       LNSV = NM1 + INT(LOG(REAL(NM1))/LOG(2.)) + 4
  66       LNWK = LOT*NM1
  67 C
  68       CALL RFFTMF(LOT,JUMP,NM1,INC,X,LENX,WSAVE(N+1),LNSV,WORK,
  69      1            LNWK,IER1)
  70       IF (IER1 .NE. 0) THEN
  71         IER = 20
  72         CALL XERFFT ('MCSTB1',-5)
  73         GO TO 106
  74       ENDIF
  75 C
  76       FNM1S2 = FLOAT(NM1)/2.
  77       M1 = 0
  78       DO 10 M=1,LJ,JUMP
  79       M1 = M1+1
  80       DSUM(M1) = .5*DSUM(M1)
  81       X(M,1) = FNM1S2*X(M,1)
  82    10 CONTINUE
  83       IF(MOD(NM1,2) .NE. 0) GO TO 30
  84       DO 20 M=1,LJ,JUMP
  85       X(M,NM1) = X(M,NM1)+X(M,NM1)
  86    20 CONTINUE
  87  30   FNM1S4 = FLOAT(NM1)/4.
  88       DO 105 I=3,N,2
  89          M1 = 0
  90          DO 115 M=1,LJ,JUMP
  91             M1 = M1+1
  92             XI = FNM1S4*X(M,I)
  93             X(M,I) = FNM1S4*X(M,I-1)
  94             X(M,I-1) = DSUM(M1)
  95             DSUM(M1) = DSUM(M1)+XI
  96   115 CONTINUE
  97   105 CONTINUE
  98       IF (MODN .NE. 0) RETURN
  99       M1 = 0
 100       DO 116 M=1,LJ,JUMP
 101          M1 = M1+1
 102          X(M,N) = DSUM(M1)
 103   116 CONTINUE
 104   106 CONTINUE
 105       RETURN
 106       END