share/fftpack5/fortran/mcsqb1.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 MCSQB1 (LOT,JUMP,N,INC,X,WSAVE,WORK,IER)
  12       DIMENSION       X(INC,*)     ,WSAVE(*)     ,WORK(LOT,*)
  13       IER = 0
  14       LJ = (LOT-1)*JUMP+1
  15       NS2 = (N+1)/2
  16       NP2 = N+2
  17       DO 101 I=3,N,2
  18          DO 201 M=1,LJ,JUMP
  19          XIM1 = X(M,I-1)+X(M,I)
  20          X(M,I) = .5*(X(M,I-1)-X(M,I))
  21          X(M,I-1) = .5*XIM1
  22  201     CONTINUE
  23   101 CONTINUE
  24       DO 301 M=1,LJ,JUMP
  25       X(M,1) = .5*X(M,1)
  26  301  CONTINUE
  27       MODN = MOD(N,2)
  28       IF (MODN .NE. 0) GO TO 302
  29       DO 303 M=1,LJ,JUMP
  30       X(M,N) = .5*X(M,N)
  31  303  CONTINUE
  32  302  CONTINUE
  33       LENX = (LOT-1)*JUMP + INC*(N-1)  + 1
  34       LNSV = N + INT(LOG(REAL(N))/LOG(2.)) + 4
  35       LNWK = LOT*N
  36 C
  37       CALL RFFTMB(LOT,JUMP,N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1)
  38       IF (IER1 .NE. 0) THEN
  39         IER = 20
  40         CALL XERFFT ('MCSQB1',-5)
  41         GO TO 400
  42       ENDIF
  43 C
  44       DO 102 K=2,NS2
  45          KC = NP2-K
  46          M1 = 0
  47          DO 202 M=1,LJ,JUMP
  48          M1 = M1 + 1
  49          WORK(M1,K) = WSAVE(K-1)*X(M,KC)+WSAVE(KC-1)*X(M,K)
  50          WORK(M1,KC) = WSAVE(K-1)*X(M,K)-WSAVE(KC-1)*X(M,KC)
  51  202     CONTINUE
  52   102 CONTINUE
  53       IF (MODN .NE. 0) GO TO 305
  54       DO 304 M=1,LJ,JUMP
  55          X(M,NS2+1) = WSAVE(NS2)*(X(M,NS2+1)+X(M,NS2+1))
  56  304     CONTINUE
  57  305  DO 103 K=2,NS2
  58          KC = NP2-K
  59          M1 = 0
  60          DO 203 M=1,LJ,JUMP
  61             M1 = M1 + 1
  62             X(M,K) = WORK(M1,K)+WORK(M1,KC)
  63             X(M,KC) = WORK(M1,K)-WORK(M1,KC)
  64  203     CONTINUE
  65   103 CONTINUE
  66       DO 104 M=1,LJ,JUMP
  67       X(M,1) = X(M,1)+X(M,1)
  68  104  CONTINUE
  69   400 CONTINUE
  70       RETURN
  71       END