# BEGIN PROLOGUE VCOST DATE WRITTEN RE

## Found at: ftp.icm.edu.pl:70/packages/netlib/vfftpack/vcost.f

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

.