SUBROUTINE VCOST(M,N,X,XT,MDIMX,WSAVE)
C***BEGIN PROLOGUE VCOST
C***DATE WRITTEN 860701 (YYMMDD)
C***REVISION DATE 900509 (YYMMDD)
C***CATEGORY NO. J1A3
C***KEYWORDS FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE
C SEQUENCES
C***AUTHOR BOISVERT, R. F. (NIST)
C***PURPOSE Cosine transform of one or more real, even sequences.
C***DESCRIPTION
C
C Subroutine VCOST computes the discrete Fourier cosine transform
C of M even sequences X(J,I), J=1,...,M. The transform is defined
C below at output parameter X.
C
C The array WSAVE which is used by subroutine VCOST must be
C initialized by calling subroutine VCOSTI(N,WSAVE).
C
C Input Parameters
C
C M the number of sequences to be transformed.
C
C N the length of the sequence to be transformed. N must be
C greater than 1. The method is most efficient when N-1 is
C is a product of small primes.
C
C X an array of size at least X(MDIMX,N) which contains the
C the sequences to be transformed. The sequences are stored
C in the ROWS of X. Thus, the Jth sequence is stored in
C X(J,I), I=1,..,N.
C
C XT a work array of size at least XT(MDIMX,N-1).
C
C MDIMX the first dimension of the array X exactly as it appears in
C the calling program.
C
C WSAVE a work array which must be dimensioned at least 3*N+15
C in the program that calls VCOST. The WSAVE array must be
C initialized by calling subroutine VCOSTI(N,WSAVE), and a
C different WSAVE array must be used for each different
C value of N. This initialization does not have to be
C repeated so long as N remains unchanged. Thus subsequent
C transforms can be obtained faster than the first.
C
C Output Parameters
C
C X For I=1,...,N and J=1,...,M
C
C X(J,I) = ( X(J,1)+(-1)**(I-1)*X(J,N)
C
C + the sum from K=2 to K=N-1
C
C 2*X(J,K)*COS((K-1)*(I-1)*PI/(N-1)) )/SQRT(2*(N-1))
C
C WSAVE contains initialization calculations which must not be
C destroyed between calls of VCOST.
C
C -----------------------------------------------------------------
C
C NOTE - A call of VCOST followed immediately by another call
C of VCOST will return the original sequences X. Thus,
C VCOST is the correctly normalized inverse of itself.
C
C -----------------------------------------------------------------
C
C VCOST is a straightforward extension of the subprogram COST to
C handle M simultaneous sequences. The scaling of the sequences
C computed by VCOST is different than that of COST. COST was
C originally developed by P. N. Swarztrauber of NCAR.
C
C***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
C Computations, (G. Rodrigue, ed.), Academic Press, 1982,
C pp. 51-83.
C***ROUTINES CALLED VRFFTF
C***END PROLOGUE VCOST
DIMENSION X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
C***FIRST EXECUTABLE STATEMENT VCOST
IF (M .LE. 0) GO TO 900
IF (N .LE. 1) GO TO 900
IF (N .GT. 3) GO TO 400
IF (N .EQ. 3) GO TO 300
C
C CASE N = 2
C
SCALE = SQRT(0.50E0)
DO 210 J=1,M
X1H = SCALE*(X(J,1)+X(J,2))
X(J,2) = SCALE*(X(J,1)-X(J,2))
X(J,1) = X1H
210 CONTINUE
GO TO 900
C
C CASE N = 3
C
300 CONTINUE
SCALE = 0.50E0
DO 310 J=1,M
X1P3 = X(J,1)+X(J,3)
TX2 = X(J,2)+X(J,2)
X(J,2) = SCALE*(X(J,1)-X(J,3))
X(J,1) = SCALE*(X1P3+TX2)
X(J,3) = SCALE*(X1P3-TX2)
310 CONTINUE
GO TO 900
C
C CASE N .GT. 3
C
C ... PREPROCESSING
C
400 CONTINUE
NM1 = N-1
NP1 = N+1
NS2 = N/2
DO 410 J=1,M
XT(J,1) = X(J,1)-X(J,N)
X(J,1) = X(J,1)+X(J,N)
410 CONTINUE
DO 420 K=2,NS2
KC = NP1-K
DO 420 J=1,M
T1 = X(J,K)+X(J,KC)
T2 = X(J,K)-X(J,KC)
XT(J,1) = XT(J,1)+WSAVE(KC)*T2
T2 = WSAVE(K)*T2
X(J,K) = T1-T2
X(J,KC) = T1+T2
420 CONTINUE
MODN = MOD(N,2)
IF (MODN .NE. 0) THEN
DO 430 J=1,M
X(J,NS2+1) = X(J,NS2+1)+X(J,NS2+1)
430 CONTINUE
ENDIF
DO 435 J=1,M
X(J,N) = XT(J,1)
435 CONTINUE
C
C ... REAL PERIODIC TRANSFORM
C
CALL VRFFTF (M,NM1,X,XT,MDIMX,WSAVE(NP1))
C
C ... POSTPROCESSING
C
FACTOR = 1.0/SQRT(REAL(NM1))
DO 440 J=1,M
XT(J,1) = X(J,2)
X(J,2) = FACTOR*X(J,N)
440 CONTINUE
DO 450 I=4,N,2
DO 450 J=1,M
XI = X(J,I)
X(J,I) = X(J,I-2)-X(J,I-1)
X(J,I-1) = XT(J,1)
XT(J,1) = XI
450 CONTINUE
IF (MODN .NE. 0) THEN
DO 460 J=1,M
X(J,N) = XT(J,1)
460 CONTINUE
ENDIF
C
C ... NORMALIZATION
C
SCALE = SQRT(0.5)
DO 490 I=1,N
DO 490 J=1,M
X(J,I) = SCALE*X(J,I)
490 CONTINUE
C
C EXIT
C
900 CONTINUE
RETURN
END
.
