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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
* SUBROUTINE 'DGPFA'
* SELF-SORTING IN-PLACE GENERALIZED PRIME FACTOR (COMPLEX) FFT
*
* *** THIS IS THE ALL-FORTRAN VERSION ***
* -------------------------------
*
* CALL DGPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN,NPQR)
*
* A IS FIRST REAL INPUT/OUTPUT VECTOR
* B IS FIRST IMAGINARY INPUT/OUTPUT VECTOR
* TRIGS IS A TABLE OF TWIDDLE FACTORS, PRECALCULATED
* BY CALLING SUBROUTINE 'SETDGPFA'
* INC IS THE INCREMENT WITHIN EACH DATA VECTOR
* JUMP IS THE INCREMENT BETWEEN DATA VECTORS
* N IS THE LENGTH OF THE TRANSFORMS:
* -----------------------------------
* N = (2**IP) * (3**IQ) * (5**IR)
* -----------------------------------
* LOT IS THE NUMBER OF TRANSFORMS
* ISIGN = +1 FOR FORWARD TRANSFORM
* = -1 FOR INVERSE TRANSFORM
* NPQR = NPQR OBTAINED FROM SETDGPFA
*
* WRITTEN BY CLIVE TEMPERTON
* RECHERCHE EN PREVISION NUMERIQUE
* ATMOSPHERIC ENVIRONMENT SERVICE, CANADA
*
* MODIFIED FOR VXL PROJECT TO ADD NPQR ARGUMENT
*
*----------------------------------------------------------------------
*
* DEFINITION OF TRANSFORM
* -----------------------
*
* X(J) = SUM(K=0,...,N-1)(C(K)*EXP(ISIGN*2*I*J*K*PI/N))
*
*---------------------------------------------------------------------
*
* FOR A MATHEMATICAL DEVELOPMENT OF THE ALGORITHM USED,
* SEE:
*
* C TEMPERTON : "A GENERALIZED PRIME FACTOR FFT ALGORITHM
* FOR ANY N = (2**P)(3**Q)(5**R)",
* SIAM J. SCI. STAT. COMP., MAY 1992.
*
*----------------------------------------------------------------------
*
SUBROUTINE DGPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN,NPQR)
*
DOUBLE PRECISION A(*), B(*), TRIGS(*)
INTEGER INC, JUMP, N, LOT, ISIGN, NPQR(3)
*
IP = NPQR(1)
IQ = NPQR(2)
IR = NPQR(3)
*
* COMPUTE THE TRANSFORM
* ---------------------
I = 1
IF (IP.GT.0) THEN
CALL DGPFA2F(A,B,TRIGS,INC,JUMP,N,IP,LOT,ISIGN)
I = I + 2 * ( 2**IP)
ENDIF
IF (IQ.GT.0) THEN
CALL DGPFA3F(A,B,TRIGS(I),INC,JUMP,N,IQ,LOT,ISIGN)
I = I + 2 * (3**IQ)
ENDIF
IF (IR.GT.0) THEN
CALL DGPFA5F(A,B,TRIGS(I),INC,JUMP,N,IR,LOT,ISIGN)
ENDIF
*
RETURN
END
|
