share/fftpack5/fortran/mcstf1.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 MCSTF1(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) 200,101,102
  20   101 DO 111 M=1,LJ,JUMP
  21          X1H = X(M,1)+X(M,2)
  22          X(M,2) = .5*(X(M,1)-X(M,2))
  23          X(M,1) = .5*X1H
  24   111 CONTINUE
  25       GO TO 200
  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          TX2 = X(M,2)+X(M,2)
  30          X(M,2) = .5*(X(M,1)-X(M,3))
  31          X(M,1) = .25*(X1P3+TX2)
  32          X(M,3) = .25*(X1P3-TX2)
  33   112 CONTINUE
  34       GO TO 200
  35   103 M1 = 0
  36       DO 113 M=1,LJ,JUMP
  37          M1 = M1+1
  38          DSUM(M1) = X(M,1)-X(M,N)
  39          X(M,1) = X(M,1)+X(M,N)
  40   113 CONTINUE
  41       DO 104 K=2,NS2
  42          M1 = 0
  43          DO 114 M=1,LJ,JUMP
  44          M1 = M1+1
  45          KC = NP1-K
  46          T1 = X(M,K)+X(M,KC)
  47          T2 = X(M,K)-X(M,KC)
  48          DSUM(M1) = DSUM(M1)+WSAVE(KC)*T2
  49          T2 = WSAVE(K)*T2
  50          X(M,K) = T1-T2
  51          X(M,KC) = T1+T2
  52   114    CONTINUE
  53   104 CONTINUE
  54       MODN = MOD(N,2)
  55       IF (MODN .EQ. 0) GO TO 124
  56          DO 123 M=1,LJ,JUMP
  57          X(M,NS2+1) = X(M,NS2+1)+X(M,NS2+1)
  58   123    CONTINUE
  59  124  CONTINUE
  60       LENX = (LOT-1)*JUMP + INC*(NM1-1)  + 1
  61       LNSV = NM1 + INT(LOG(REAL(NM1))/LOG(2.)) + 4
  62       LNWK = LOT*NM1
  63 C
  64       CALL RFFTMF(LOT,JUMP,NM1,INC,X,LENX,WSAVE(N+1),LNSV,WORK,
  65      1            LNWK,IER1)
  66       IF (IER1 .NE. 0) THEN
  67         IER = 20
  68         CALL XERFFT ('MCSTF1',-5)
  69         GO TO 200
  70       ENDIF
  71 C
  72       SNM1 = 1./FLOAT(NM1)
  73       DO 10 M=1,LOT
  74       DSUM(M) = SNM1*DSUM(M)
  75    10 CONTINUE
  76       IF(MOD(NM1,2) .NE. 0) GO TO 30
  77       DO 20 M=1,LJ,JUMP
  78       X(M,NM1) = X(M,NM1)+X(M,NM1)
  79    20 CONTINUE
  80  30   DO 105 I=3,N,2
  81          M1 = 0
  82          DO 115 M=1,LJ,JUMP
  83             M1 = M1+1
  84             XI = .5*X(M,I)
  85             X(M,I) = .5*X(M,I-1)
  86             X(M,I-1) = DSUM(M1)
  87             DSUM(M1) = DSUM(M1)+XI
  88   115 CONTINUE
  89   105 CONTINUE
  90       IF (MODN .NE. 0) GO TO 117
  91       M1 = 0
  92       DO 116 M=1,LJ,JUMP
  93          M1 = M1+1
  94          X(M,N) = DSUM(M1)
  95   116 CONTINUE
  96  117  DO 118 M=1,LJ,JUMP
  97       X(M,1) = .5*X(M,1)
  98       X(M,N) = .5*X(M,N)
  99  118  CONTINUE
 100 C
 101  200  CONTINUE
 102       RETURN
 103       END