1 /* temperton/gpfa.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static integer c__2 = 2;
21 static integer c__3 = 3;
22
23 /* SUBROUTINE 'GPFA' */
24 /* SELF-SORTING IN-PLACE GENERALIZED PRIME FACTOR (COMPLEX) FFT */
25
26 /* *** THIS IS THE ALL-FORTRAN VERSION *** */
27 /* ------------------------------- */
28
29 /* CALL GPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN,NPQR) */
30
31 /* A IS FIRST REAL INPUT/OUTPUT VECTOR */
32 /* B IS FIRST IMAGINARY INPUT/OUTPUT VECTOR */
33 /* TRIGS IS A TABLE OF TWIDDLE FACTORS, PRECALCULATED */
34 /* BY CALLING SUBROUTINE 'SETGPFA' */
35 /* INC IS THE INCREMENT WITHIN EACH DATA VECTOR */
36 /* JUMP IS THE INCREMENT BETWEEN DATA VECTORS */
37 /* N IS THE LENGTH OF THE TRANSFORMS: */
38 /* ----------------------------------- */
39 /* N = (2**IP) * (3**IQ) * (5**IR) */
40 /* ----------------------------------- */
41 /* LOT IS THE NUMBER OF TRANSFORMS */
42 /* ISIGN = +1 FOR FORWARD TRANSFORM */
43 /* = -1 FOR INVERSE TRANSFORM */
44 /* NPQR = NPQR OBTAINED FROM SETGPFA */
45
46 /* WRITTEN BY CLIVE TEMPERTON */
47 /* RECHERCHE EN PREVISION NUMERIQUE */
48 /* ATMOSPHERIC ENVIRONMENT SERVICE, CANADA */
49
50 /* MODIFIED FOR VXL PROJECT TO ADD NPQR ARGUMENT */
51
52 /* ---------------------------------------------------------------------- */
53
54 /* DEFINITION OF TRANSFORM */
55 /* ----------------------- */
56
57 /* X(J) = SUM(K=0,...,N-1)(C(K)*EXP(ISIGN*2*I*J*K*PI/N)) */
58
59 /* --------------------------------------------------------------------- */
60
61 /* FOR A MATHEMATICAL DEVELOPMENT OF THE ALGORITHM USED, */
62 /* SEE: */
63
64 /* C TEMPERTON : "A GENERALIZED PRIME FACTOR FFT ALGORITHM */
65 /* FOR ANY N = (2**P)(3**Q)(5**R)", */
66 /* SIAM J. SCI. STAT. COMP., MAY 1992. */
67
68 /* ---------------------------------------------------------------------- */
69
70 /*< SUBROUTINE GPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN,NPQR) >*/
gpfa_(real * a,real * b,real * trigs,integer * inc,integer * jump,integer * n,integer * lot,integer * isign,integer * npqr)71 /* Subroutine */ int gpfa_(real *a, real *b, real *trigs, integer *inc,
72 integer *jump, integer *n, integer *lot, integer *isign, integer *
73 npqr)
74 {
75 /* Builtin functions */
76 integer pow_ii(integer *, integer *);
77
78 /* Local variables */
79 integer i__, ip, iq, ir;
80 extern /* Subroutine */ int gpfa2f_(real *, real *, real *, integer *,
81 integer *, integer *, integer *, integer *, integer *), gpfa3f_(
82 real *, real *, real *, integer *, integer *, integer *, integer *
83 , integer *, integer *), gpfa5f_(real *, real *, real *, integer *
84 , integer *, integer *, integer *, integer *, integer *);
85
86
87 /*< REAL A(*), B(*), TRIGS(*) >*/
88 /*< INTEGER INC, JUMP, N, LOT, ISIGN, NPQR(3) >*/
89
90 /*< IP = NPQR(1) >*/
91 /* Parameter adjustments */
92 --npqr;
93 --trigs;
94 --b;
95 --a;
96
97 /* Function Body */
98 ip = npqr[1];
99 /*< IQ = NPQR(2) >*/
100 iq = npqr[2];
101 /*< IR = NPQR(3) >*/
102 ir = npqr[3];
103
104 /* COMPUTE THE TRANSFORM */
105 /* --------------------- */
106 /*< I = 1 >*/
107 i__ = 1;
108 /*< IF (IP.GT.0) THEN >*/
109 if (ip > 0) {
110 /*< CALL GPFA2F(A,B,TRIGS,INC,JUMP,N,IP,LOT,ISIGN) >*/
111 gpfa2f_(&a[1], &b[1], &trigs[1], inc, jump, n, &ip, lot, isign);
112 /*< I = I + 2 * ( 2**IP) >*/
113 i__ += pow_ii(&c__2, &ip) << 1;
114 /*< ENDIF >*/
115 }
116 /*< IF (IQ.GT.0) THEN >*/
117 if (iq > 0) {
118 /*< CALL GPFA3F(A,B,TRIGS(I),INC,JUMP,N,IQ,LOT,ISIGN) >*/
119 gpfa3f_(&a[1], &b[1], &trigs[i__], inc, jump, n, &iq, lot, isign);
120 /*< I = I + 2 * (3**IQ) >*/
121 i__ += pow_ii(&c__3, &iq) << 1;
122 /*< ENDIF >*/
123 }
124 /*< IF (IR.GT.0) THEN >*/
125 if (ir > 0) {
126 /*< CALL GPFA5F(A,B,TRIGS(I),INC,JUMP,N,IR,LOT,ISIGN) >*/
127 gpfa5f_(&a[1], &b[1], &trigs[i__], inc, jump, n, &ir, lot, isign);
128 /*< ENDIF >*/
129 }
130
131 /*< RETURN >*/
132 return 0;
133 /*< END >*/
134 } /* gpfa_ */
135
136 #ifdef __cplusplus
137 }
138 #endif
139