cosqb1.f 1.31 KB
Newer Older
rtoy's avatar
rtoy committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   FFTPACK 5.0 
C
C   Authors:  Paul N. Swarztrauber and Richard A. Valent
C
C   $Id$
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE COSQB1 (N,INC,X,WSAVE,WORK,IER)
      DIMENSION       X(INC,*)     ,WSAVE(*)     ,WORK(*)
      IER = 0
      NS2 = (N+1)/2
      NP2 = N+2
      DO 101 I=3,N,2
         XIM1 = X(1,I-1)+X(1,I)
         X(1,I) = .5*(X(1,I-1)-X(1,I))
         X(1,I-1) = .5*XIM1
  101 CONTINUE
      X(1,1) = .5*X(1,1)
      MODN = MOD(N,2)
      IF (MODN .NE. 0) GO TO 302
      X(1,N) = .5*X(1,N)
  302 LENX = INC*(N-1)  + 1
Raymond Toy's avatar
Raymond Toy committed
26
      LNSV = N + INT(LOG(REAL(N))/LOG(2.)) + 4
rtoy's avatar
rtoy committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
      LNWK = N
C
      CALL RFFT1B(N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1)
      IF (IER1 .NE. 0) THEN
        IER = 20
        CALL XERFFT ('COSQB1',-5)
        GO TO 400
      ENDIF
C
      DO 102 K=2,NS2
         KC = NP2-K
         WORK(K) = WSAVE(K-1)*X(1,KC)+WSAVE(KC-1)*X(1,K)
         WORK(KC) = WSAVE(K-1)*X(1,K)-WSAVE(KC-1)*X(1,KC)
  102 CONTINUE
      IF (MODN .NE. 0) GO TO 305
      X(1,NS2+1) = WSAVE(NS2)*(X(1,NS2+1)+X(1,NS2+1))
  305 DO 103 K=2,NS2
         KC = NP2-K
         X(1,K) = WORK(K)+WORK(KC)
         X(1,KC) = WORK(K)-WORK(KC)
  103 CONTINUE
      X(1,1) = X(1,1)+X(1,1)
  400 RETURN
      END