Main Page Class Hierarchy Compound List File List Compound Members File Members

# fft99.c

Go to the documentation of this file.

    00001 /*
    00002  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
    00003  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
    00004  *
    00005  *    $Id: fft99_8c-source.html,v 1.1 2006/05/19 16:55:01 paulf Exp $
    00006  *
    00007  *    Revision history:
    00008  *     $Log: fft99_8c-source.html,v $
    00008  *     Revision 1.1  2006/05/19 16:55:01  paulf
    00008  *     first inclusion
    00008  *
    00009  *     Revision 1.1  2001/03/31 00:46:48  lombard
    00010  *     Initial revision
    00011  *
    00012  *
    00013  *
    00014  */
    00015
    00016 /* fft99.f -- translated by f2c (version 19991025).*/
    00017 /*
    00018  * Hand-editted after f2c to make it look a little more like C.
    00019  * But there is no hiding the fact that this was once Fortran.
    00020  * Read the comments with caution! Pete Lombard, January 2001
    00021  */
    00022 #include <math.h>
    00023 #include "fft99.h"
    00024
    00025 /* Internal function prototypes */
    00026 static void fact(long n, long *ifax);
    00027 static void cftrig(long n, double *trigs);
    00028 static void vpassm(double *a, double *b, double *c, double *d,
    00029                    double *trigs, long inc1, long inc2, long inc3, long inc4,
    00030                    long lot, long n, long ifac, long la);
    00031 static void fft99a(double *a, double *work, double *trigs, long inc, long jump,
    00032                    long n, long lot);
    00033 static void fft99b(double *work, double *a, double *trigs, long inc, long jump,
    00034                    long n, long lot);
    00035 static void fax(long *ifax, long n);
    00036 static void fftrig(double *trigs, long n);
    00037
    00038
    00039 /* Function definitions */
    00040 void cfft99(double *a, double *work, double *trigs, long *ifax, long inc,
    00041             long jump, long n, long lot, long isign)
    00042 {
    00043   /* Local variables */
    00044   static long nfax, i, j, k, l, m, ibase, jbase;
    00045   static double himag, hreal;
    00046   static long ilast, i1, i2, la, nh, nn;
    00047   static long igo, ink, jnk, jum, jst;
    00048
    00049
    00050   /* PURPOSE      PERFORMS MULTIPLE FAST FOURIER TRANSFORMS.  THIS PACKAGE
    00051  *              WILL PERFORM A NUMBER OF SIMULTANEOUS COMPLEX PERIODIC
    00052  *              FOURIER TRANSFORMS OR CORRESPONDING INVERSE TRANSFORMS.
    00053  *              THAT IS, GIVEN A SET OF COMPLEX GRIDPOINT VECTORS, THE
    00054  *              PACKAGE RETURNS A SET OF COMPLEX FOURIER
    00055  *              COEFFICIENT VECTORS, OR VICE VERSA.  THE LENGTH OF THE
    00056  *              TRANSFORMS MUST BE A NUMBER GREATER THAN 1 THAT HAS
    00057  *              NO PRIME FACTORS OTHER THAN 2, 3, AND 5.
    00058  *
    00059  *              THE PACKAGE CFFT99 CONTAINS SEVERAL USER-LEVEL ROUTINES:
    00060  *
    00061  *            CFTFAX
    00062  *                AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE
    00063  *                BEFORE A SEQUENCE OF CALLS TO CFFT99
    00064  *                (PROVIDED THAT N IS NOT CHANGED).
    00065  *
    00066  *            CFFT99
    00067  *                THE ACTUAL TRANSFORM ROUTINE ROUTINE, CABABLE OF
    00068  *                PERFORMING BOTH THE TRANSFORM AND ITS INVERSE.
    00069  *                HOWEVER, AS THE TRANSFORMS ARE NOT NORMALIZED,
    00070  *                THE APPLICATION OF A TRANSFORM FOLLOWED BY ITS
    00071  *                INVERSE WILL YIELD THE ORIGINAL VALUES MULTIPLIED
    00072  *                BY N.
    00073  *
    00074  *
    00075  *
    00076  * USAGE        LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 0,
    00077  *              Q .GE. 0, AND R .GE. 0.  THEN A TYPICAL SEQUENCE OF
    00078  *              CALLS TO TRANSFORM A GIVEN SET OF COMPLEX VECTORS OF
    00079  *              LENGTH N TO A SET OF (UNSCALED) COMPLEX FOURIER
    00080  *              COEFFICIENT VECTORS OF LENGTH N IS
    00081  *
    00082  *                   long ifax[13];
    00083  *                   double trigs[2*n];
    00084  *                   double a[...], work[...]
    00085  *
    00086  *                   cftfax (n, ifax, trigs);
    00087  *                   cfft99 (a,work,trigs,ifax,inc,jump,n,lot,isign);
    00088  *
    00089  *              THE OUTPUT VECTORS OVERWRITE THE INPUT VECTORS, AND
    00090  *              THESE ARE STORED IN A.  WITH APPROPRIATE CHOICES FOR
    00091  *              THE OTHER ARGUMENTS, THESE VECTORS MAY BE CONSIDERED
    00092  *              EITHER THE ROWS OR THE COLUMNS OF THE ARRAY A.
    00093  *              SEE THE INDIVIDUAL WRITE-UPS FOR CFTFAX AND
    00094  *              CFFT99 BELOW, FOR A DETAILED DESCRIPTION OF THE
    00095  *              ARGUMENTS.
    00096  *
    00097  * HISTORY      THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN
    00098  *              NOVEMBER, 1978.  IT WAS MODIFIED, DOCUMENTED, AND TESTED
    00099  *              FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.  IT WAS
    00100  *              FURTHER MODIFIED FOR THE FULLY COMPLEX CASE BY DAVE
    00101  *              FULKER IN NOVEMBER, 1980.
    00102  *
    00103  * -----------------------------------------------------------------------
    00104  *
    00105  *  cftfax (n,ifax,trigs)
    00106  *
    00107  * PURPOSE      A SET-UP ROUTINE FOR CFFT99.  IT NEED ONLY BE
    00108  *              CALLED ONCE BEFORE A SEQUENCE OF CALLS TO CFFT99,
    00109  *              PROVIDED THAT N IS NOT CHANGED.
    00110  *
    00111  * ARGUMENT     long ifax[13]; double trigs[2*n];
    00112  * DIMENSIONS
    00113  *
    00114  * ARGUMENTS
    00115  *
    00116  * ON INPUT     N
    00117  *               AN EVEN NUMBER GREATER THAN 1 THAT HAS NO PRIME FACTOR
    00118  *               GREATER THAN 5.  N IS THE LENGTH OF THE TRANSFORMS (SEE
    00119  *               THE DOCUMENTATION FOR CFFT99 FOR THE DEFINITION OF
    00120  *               THE TRANSFORMS).
    00121  *
    00122  *              IFAX
    00123  *               AN INTEGER ARRAY.  THE NUMBER OF ELEMENTS ACTUALLY USED
    00124  *               WILL DEPEND ON THE FACTORIZATION OF N.  DIMENSIONING
    00125  *               IFAX FOR 13 SUFFICES FOR ALL N LESS THAN 1 MILLION.
    00126  *
    00127  *              TRIGS
    00128  *               AN ARRAY OF DIMENSION 2*N
    00129  *
    00130  * ON OUTPUT    IFAX
    00131  *               CONTAINS THE FACTORIZATION OF N.  IFAX[0] IS THE
    00132  *               NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED
    00133  *               IN IFAX[1],IFAX[2],...  IF N HAS ANY PRIME FACTORS
    00134  *               GREATER THAN 5, IFAX[0] IS SET TO -99.
    00135  *
    00136  *              TRIGS
    00137  *               AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY
    00138  *               USED BY THE CFT ROUTINES.
    00139  *
    00140  * -----------------------------------------------------------------------
    00141  *
    00142  * cfft99 (a,work,trigs,ifax,inc,jump,n,lot,isign)
    00143  *
    00144  * PURPOSE      PERFORM A NUMBER OF SIMULTANEOUS (UNNORMALIZED) COMPLEX
    00145  *              PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
    00146  *              TRANSFORMS.  GIVEN A SET OF COMPLEX GRIDPOINT
    00147  *              VECTORS, THE PACKAGE RETURNS A SET OF
    00148  *              COMPLEX FOURIER COEFFICIENT VECTORS, OR VICE
    00149  *              VERSA.  THE LENGTH OF THE TRANSFORMS MUST BE A
    00150  *              NUMBER HAVING NO PRIME FACTORS OTHER THAN
    00151  *              2, 3, AND 5.  THIS ROUTINE IS
    00152  *              OPTIMIZED FOR USE ON THE CRAY-1.
    00153  *
    00154  * ARGUMENT     double a[2*n*inc+(lot-1)*jump], work[2*n*lot]
    00155  * DIMENSIONS   double trigs[2*n], long ifax[13]
    00156  *
    00157  * ARGUMENTS
    00158  *
    00159  * ON INPUT     A
    00160  *               A COMPLEX ARRAY OF LENGTH N*INC+(LOT-1)*JUMP CONTAINING
    00161  *               THE INPUT GRIDPOINT OR COEFFICIENT VECTORS.  THIS ARRAY
    00162  *               OVERWRITTEN BY THE RESULTS.
    00163  *
    00164  *               N.B. ALTHOUGH THE ARRAY A IS USUALLY CONSIDERED TO BE OF
    00165  *               TYPE COMPLEX IN THE CALLING PROGRAM, IT IS TREATED AS
    00166  *               REAL WITHIN THE TRANSFORM PACKAGE.  THIS REQUIRES THAT
    00167  *               SUCH TYPE CONFLICTS ARE PERMITTED IN THE USER^S
    00168  *               ENVIRONMENT, AND THAT THE STORAGE OF COMPLEX NUMBERS
    00169  *               MATCHES THE ASSUMPTIONS OF THIS ROUTINE.  THIS ROUTINE
    00170  *               ASSUMES THAT THE REAL AND IMAGINARY PORTIONS OF A
    00171  *               COMPLEX NUMBER OCCUPY ADJACENT ELEMENTS OF MEMORY.  IF
    00172  *               THESE CONDITIONS ARE NOT MET, THE USER MUST TREAT THE
    00173  *               ARRAY A AS REAL (AND OF TWICE THE ABOVE LENGTH), AND
    00174  *               WRITE THE CALLING PROGRAM TO TREAT THE REAL AND
    00175  *               IMAGINARY PORTIONS EXPLICITLY.
    00176  *
    00177  *              WORK
    00178  *               A COMPLEX WORK ARRAY OF LENGTH N*LOT OR A REAL ARRAY
    00179  *               OF LENGTH 2*N*LOT.  SEE N.B. ABOVE.
    00180  *
    00181  *              TRIGS
    00182  *               AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST.
    00183  *
    00184  *              IFAX
    00185  *               AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST.
    00186  *
    00187  *
    00188  *               N.B. IN THE FOLLOWING ARGUMENTS, INCREMENTS ARE MEASURED
    00189  *               IN WORD PAIRS, BECAUSE EACH COMPLEX ELEMENT IS ASSUMED
    00190  *               TO OCCUPY AN ADJACENT PAIR OF WORDS IN MEMORY.
    00191  *
    00192  *              INC
    00193  *               THE INCREMENT (IN WORD PAIRS) BETWEEN SUCCESSIVE ELEMENT
    00194  *               OF EACH (COMPLEX) GRIDPOINT OR COEFFICIENT VECTOR
    00195  *               (E.G.  INC=1 FOR CONSECUTIVELY STORED DATA).
    00196  *
    00197  *              JUMP
    00198  *               THE INCREMENT (IN WORD PAIRS) BETWEEN THE FIRST ELEMENTS
    00199  *               OF SUCCESSIVE DATA OR COEFFICIENT VECTORS.  ON THE CRAY-
    00200  *               TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8
    00201  *               (TO AVOID MEMORY BANK CONFLICTS).  FOR CLARIFICATION OF
    00202  *               INC AND JUMP, SEE THE EXAMPLES BELOW.
    00203  *
    00204  *              N
    00205  *               THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF
    00206  *               TRANSFORMS, BELOW).
    00207  *
    00208  *              LOT
    00209  *               THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY.
    00210  *
    00211  *              ISIGN
    00212  *               = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER
    00213  *                    COEFFICIENTS.
    00214  *               = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO
    00215  *                    GRIDPOINT VALUES.
    00216  *
    00217  * ON OUTPUT    A
    00218  *               IF ISIGN = -1, AND LOT GRIDPOINT VECTORS ARE SUPPLIED,
    00219  *               EACH CONTAINING THE COMPLEX SEQUENCE:
    00220  *
    00221  *               G(0),G(1), ... ,G(N-1)  (N COMPLEX VALUES)
    00222  *
    00223  *               THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH
    00224  *               CONTAINING THE CORRESPONDING N COEFFICIENT VALUES:
    00225  *
    00226  *               C(0),C(1), ... ,C(N-1)  (N COMPLEX VALUES)
    00227  *
    00228  *               DEFINED BY:
    00229  *                 C(K) = SUM(J=0,...,N-1)( G(J)*EXP(-2*I*J*K*PI/N) )
    00230  *                 WHERE I = SQRT(-1)
    00231  *
    00232  *
    00233  *               IF ISIGN = +1, AND LOT COEFFICIENT VECTORS ARE SUPPLIED,
    00234  *               EACH CONTAINING THE COMPLEX SEQUENCE:
    00235  *
    00236  *               C(0),C(1), ... ,C(N-1)  (N COMPLEX VALUES)
    00237  *
    00238  *               THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH
    00239  *               CONTAINING THE CORRESPONDING N GRIDPOINT VALUES:
    00240  *
    00241  *               G(0),G(1), ... ,G(N-1)  (N COMPLEX VALUES)
    00242  *
    00243  *               DEFINED BY:
    00244  *                 G(J) = SUM(K=0,...,N-1)( G(K)*EXP(+2*I*J*K*PI/N) )
    00245  *                 WHERE I = SQRT(-1)
    00246  *
    00247  *
    00248  *               A CALL WITH ISIGN=-1 FOLLOWED BY A CALL WITH ISIGN=+1
    00249  *               (OR VICE VERSA) RETURNS THE ORIGINAL DATA, MULTIPLIED
    00250  *               BY THE FACTOR N.
    00251  *
    00252  *
    00253  * EXAMPLE       GIVEN A 64 BY 9 GRID OF COMPLEX VALUES, STORED IN
    00254  *               A 66 BY 9 COMPLEX ARRAY, A, COMPUTE THE TWO DIMENSIONAL
    00255  *               FOURIER TRANSFORM OF THE GRID.  FROM TRANSFORM THEORY,
    00256  *               IT IS KNOWN THAT A TWO DIMENSIONAL TRANSFORM CAN BE
    00257  *               OBTAINED BY FIRST TRANSFORMING THE GRID ALONG ONE
    00258  *               DIRECTION, THEN TRANSFORMING THESE RESULTS ALONG THE
    00259  *               ORTHOGONAL DIRECTION.
    00260  *
    00261  *               COMPLEX A(66,9), WORK(64,9)
    00262  *               REAL TRIGS1(128), TRIGS2(18)
    00263  *               INTEGER IFAX1(13), IFAX2(13)
    00264  *
    00265  *               SET UP THE IFAX AND TRIGS ARRAYS FOR EACH DIRECTION:
    00266  *
    00267  *               CALL CFTFAX(64, IFAX1, TRIGS1)
    00268  *               CALL CFTFAX( 9, IFAX2, TRIGS2)
    00269  *
    00270  *               IN THIS CASE, THE COMPLEX VALUES OF THE GRID ARE
    00271  *               STORED IN MEMORY AS FOLLOWS (USING U AND V TO
    00272  *               DENOTE THE REAL AND IMAGINARY COMPONENTS, AND
    00273  *               ASSUMING CONVENTIONAL FORTRAN STORAGE):
    00274  *
    00275  *   U(1,1), V(1,1), U(2,1), V(2,1),  ...  U(64,1), V(64,1), 4 NULLS,
    00276  *
    00277  *   U(1,2), V(1,2), U(2,2), V(2,2),  ...  U(64,2), V(64,2), 4 NULLS,
    00278  *
    00279  *   .       .       .       .         .   .        .        .
    00280  *   .       .       .       .         .   .        .        .
    00281  *   .       .       .       .         .   .        .        .
    00282  *
    00283  *   U(1,9), V(1,9), U(2,9), V(2,9),  ...  U(64,9), V(64,9), 4 NULLS.
    00284  *
    00285  *               WE CHOOSE (ARBITRARILY) TO TRANSORM FIRST ALONG THE
    00286  *               DIRECTION OF THE FIRST SUBSCRIPT.  THUS WE DEFINE
    00287  *               THE LENGTH OF THE TRANSFORMS, N, TO BE 64, THE
    00288  *               NUMBER OF TRANSFORMS, LOT, TO BE 9, THE INCREMENT
    00289  *               BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 1,
    00290  *               AND THE INCREMENT BETWEEN THE STARTING POINTS
    00291  *               FOR EACH TRANSFORM, JUMP, TO BE 66 (THE FIRST
    00292  *               DIMENSION OF A).
    00293  *
    00294  *               CALL CFFT99( A, WORK, TRIGS1, IFAX1, 1, 66, 64, 9, -1)
    00295  *
    00296  *               TO TRANSFORM ALONG THE DIRECTION OF THE SECOND SUBSCRIPT
    00297  *               THE ROLES OF THE INCREMENTS ARE REVERSED.  THUS WE DEFIN
    00298  *               THE LENGTH OF THE TRANSFORMS, N, TO BE 9, THE
    00299  *               NUMBER OF TRANSFORMS, LOT, TO BE 64, THE INCREMENT
    00300  *               BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 66,
    00301  *               AND THE INCREMENT BETWEEN THE STARTING POINTS
    00302  *               FOR EACH TRANSFORM, JUMP, TO BE 1
    00303  *
    00304  *               CALL CFFT99( A, WORK, TRIGS2, IFAX2, 66, 1, 9, 64, -1)
    00305  *
    00306  *               THESE TWO SEQUENTIAL STEPS RESULTS IN THE TWO-DIMENSIONA
    00307  *               FOURIER COEFFICIENT ARRAY OVERWRITING THE INPUT
    00308  *               GRIDPOINT ARRAY, A.  THE SAME TWO STEPS APPLIED AGAIN
    00309  *               WITH ISIGN = +1 WOULD RESULT IN THE RECONSTRUCTION OF
    00310  *               THE GRIDPOINT ARRAY (MULTIPLIED BY A FACTOR OF 64*9).
    00311  *
    00312  *
    00313  * -----------------------------------------------------------------------
    00314  */
    00315
    00316 /*     function ^cfft99^ - MULTIPLE FAST COMPLEX FOURIER TRANSFORM
    00317
    00318  *     A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
    00319  *     WORK IS AN AREA OF SIZE N*LOT
    00320  *     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
    00321  *     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
    00322  *     INC IS THE INCREMENT WITHIN EACH DATA #VECTOR#
    00323  *         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
    00324  *     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
    00325  *     N IS THE LENGTH OF THE DATA VECTORS
    00326  *     LOT IS THE NUMBER OF DATA VECTORS
    00327  *     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
    00328  *           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
    00329  */
    00330
    00331 /*  VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN  PARALLEL. */
    00332
    00333     /* Parameter adjustments */
    00334   --work;
    00335   --a;
    00336
    00337     /* Function Body */
    00338   nn = n + n;
    00339   ink = inc + inc;
    00340   jum = jump + jump;
    00341   nfax = ifax[0];
    00342   jnk = 2;
    00343   jst = 2;
    00344   if (isign >= 0) {
    00345     goto L30;
    00346   }
    00347
    00348   /*     THE INNERMOST TEMPERTON ROUTINES HAVE NO FACILITY FOR THE */
    00349   /*     FORWARD (ISIGN = -1) TRANSFORM.  THEREFORE, THE INPUT MUST BE */
    00350   /*     REARRANGED AS FOLLOWS: */
    00351
    00352 /*     THE ORDER OF EACH INPUT VECTOR, */
    00353
    00354 /*     G(0), G(1), G(2), ... , G(N-2), G(N-1) */
    00355
    00356 /*     IS REVERSED (EXCLUDING G(0)) TO YIELD */
    00357
    00358 /*     G(0), G(N-1), G(N-2), ... , G(2), G(1). */
    00359
    00360 /*     WITHIN THE TRANSFORM, THE CORRESPONDING EXPONENTIAL MULTIPLIER */
    00361 /*     IS THEN PRECISELY THE CONJUGATE OF THAT FOR THE NORMAL */
    00362 /*     ORDERING.  THUS THE FORWARD (ISIGN = -1) TRANSFORM IS */
    00363 /*     ACCOMPLISHED */
    00364
    00365 /*     FOR NFAX ODD, THE INPUT MUST BE TRANSFERRED TO THE WORK ARRAY, */
    00366 /*     AND THE REARRANGEMENT CAN BE DONE DURING THE MOVE. */
    00367
    00368   jnk = -2;
    00369   jst = nn - 2;
    00370   if (nfax % 2 == 1) {
    00371     goto L40;
    00372   }
    00373
    00374   /*     FOR NFAX EVEN, THE REARRANGEMENT MUST BE APPLIED DIRECTLY TO */
    00375   /*     THE INPUT ARRAY.  THIS CAN BE DONE BY SWAPPING ELEMENTS. */
    00376
    00377   ibase = 1;
    00378   ilast = (n - 1) * ink;
    00379   nh = n / 2;
    00380   for (l = 1; l <= lot; ++l) {
    00381     i1 = ibase + ink;
    00382     i2 = ibase + ilast;
    00383     /* DIR$ IVDEP */
    00384     for (m = 1; m <= nh; ++m) {
    00385       /*     SWAP REAL AND IMAGINARY PORTIONS */
    00386       hreal = a[i1];
    00387       himag = a[i1 + 1];
    00388       a[i1] = a[i2];
    00389       a[i1 + 1] = a[i2 + 1];
    00390       a[i2] = hreal;
    00391       a[i2 + 1] = himag;
    00392       i1 += ink;
    00393       i2 -= ink;
    00394       /* L10: */
    00395     }
    00396     ibase += jum;
    00397     /* L20: */
    00398   }
    00399   goto L100;
    00400
    00401  L30:
    00402   if (nfax % 2 == 0) {
    00403     goto L100;
    00404   }
    00405
    00406  L40:
    00407
    00408 /*     DURING THE TRANSFORM PROCESS, NFAX STEPS ARE TAKEN, AND THE */
    00409 /*     RESULTS ARE STORED ALTERNATELY IN WORK AND IN A.  IF NFAX IS */
    00410 /*     ODD, THE INPUT DATA ARE FIRST MOVED TO WORK SO THAT THE FINAL */
    00411 /*     RESULT (AFTER NFAX STEPS) IS STORED IN ARRAY A. */
    00412
    00413   ibase = 1;
    00414   jbase = 1;
    00415   for (l = 1; l <= lot; ++l) {
    00416     /*     MOVE REAL AND IMAGINARY PORTIONS OF ELEMENT ZERO */
    00417     work[jbase] = a[ibase];
    00418     work[jbase + 1] = a[ibase + 1];
    00419     i = ibase + ink;
    00420     j = jbase + jst;
    00421     /* DIR$ IVDEP */
    00422     for (m = 2; m <= n; ++m) {
    00423       /*     MOVE REAL AND IMAGINARY PORTIONS OF OTHER ELEMENTS (POSSIBLY IN */
    00424       /*     REVERSE ORDER, DEPENDING ON JST AND JNK) */
    00425       work[j] = a[i];
    00426       work[j + 1] = a[i + 1];
    00427       i += ink;
    00428       j += jnk;
    00429       /* L50: */
    00430     }
    00431     ibase += jum;
    00432     jbase += nn;
    00433     /* L60: */
    00434   }
    00435
    00436  L100:
    00437
    00438 /*     PERFORM THE TRANSFORM PASSES, ONE PASS FOR EACH FACTOR.  DURING */
    00439 /*     EACH PASS THE DATA ARE MOVED FROM A TO WORK OR FROM WORK TO A. */
    00440
    00441 /*     FOR NFAX EVEN, THE FIRST PASS MOVES FROM A TO WORK */
    00442   igo = 110;
    00443   /*     FOR NFAX ODD, THE FIRST PASS MOVES FROM WORK TO A */
    00444   if (nfax % 2 == 1) {
    00445     igo = 120;
    00446   }
    00447   la = 1;
    00448   for (k = 1; k <= nfax; ++k) {
    00449     if (igo == 120) {
    00450       goto L120;
    00451     }
    00452     /* L110: */
    00453     vpassm(&a[1], &a[2], &work[1], &work[2], trigs, ink, 2, jum, nn, lot, n,
    00454            ifax[k], la);
    00455     igo = 120;
    00456     goto L130;
    00457   L120:
    00458     vpassm(&work[1], &work[2], &a[1], &a[2], trigs, 2, ink, nn,
    00459            jum, lot, n, ifax[k], la);
    00460     igo = 110;
    00461   L130:
    00462     la *= ifax[k];
    00463     /* L140: */
    00464   }
    00465
    00466   /*     AT THIS POINT THE FINAL TRANSFORM RESULT IS STORED IN A. */
    00467
    00468   return;
    00469 } /* cfft99 */
    00470
    00471 void cftfax(long n, long *ifax, double *trigs)
    00472 {
    00473   /* Local variables */
    00474   static long k;
    00475
    00476 /*     THIS ROUTINE WAS MODIFIED FROM TEMPERTON^S ORIGINAL */
    00477 /*     BY DAVE FULKER.  IT NO LONGER PRODUCES FACTORS IN ASCENDING */
    00478 /*     ORDER, AND THERE ARE NONE OF THE ORIGINAL #MODE# OPTIONS. */
    00479
    00480 /* ON INPUT     N */
    00481 /*               THE LENGTH OF EACH COMPLEX TRANSFORM TO BE PERFORMED */
    00482
    00483 /*               N MUST BE GREATER THAN 1 AND CONTAIN NO PRIME */
    00484 /*               FACTORS GREATER THAN 5. */
    00485
    00486 /* ON OUTPUT    IFAX */
    00487 /*               IFAX(1) */
    00488 /*                 THE NUMBER OF FACTORS CHOSEN OR -99 IN CASE OF ERROR */
    00489 /*               IFAX(2) THRU IFAX( IFAX(1)+1 ) */
    00490 /*                 THE FACTORS OF N IN THE FOLLOWIN ORDER:  APPEARING */
    00491 /*                 FIRST ARE AS MANY FACTORS OF 4 AS CAN BE OBTAINED. */
    00492 /*                 SUBSEQUENT FACTORS ARE PRIMES, AND APPEAR IN */
    00493 /*                 ASCENDING ORDER, EXCEPT FOR MULTIPLE FACTORS. */
    00494
    00495 /*              TRIGS */
    00496 /*               2N SIN AND COS VALUES FOR USE BY THE TRANSFORM ROUTINE */
    00497
    00498     /* Function Body */
    00499   fact(n, ifax);
    00500   k = ifax[0];
    00501   if (k < 1 || ifax[k] > 5) {
    00502     ifax[0] = -99;
    00503   }
    00504   if (ifax[0] <= 0)
    00505     return;
    00506
    00507   cftrig(n, trigs);
    00508   return;
    00509 } /* cftfax */
    00510
    00511 static void fact(long n, long *ifax)
    00512 {
    00513   /* Local variables */
    00514   static long k, l, nn, inc, max;
    00515
    00516 /*     FACTORIZATION ROUTINE THAT FIRST EXTRACTS ALL FACTORS OF 4 */
    00517     /* Function Body */
    00518   if (n > 1) {
    00519     goto L10;
    00520   }
    00521   ifax[0] = 0;
    00522   if (n < 1) {
    00523     ifax[0] = -99;
    00524   }
    00525   return;
    00526  L10:
    00527   nn = n;
    00528   k = 0;
    00529   /*     TEST FOR FACTORS OF 4 */
    00530  L20:
    00531   if (nn % 4 != 0) {
    00532     goto L30;
    00533   }
    00534   ++k;
    00535   ifax[k] = 4;
    00536   nn /= 4;
    00537   if (nn == 1) {
    00538     goto L80;
    00539   }
    00540   goto L20;
    00541   /*     TEST FOR EXTRA FACTOR OF 2 */
    00542  L30:
    00543   if (nn % 2 != 0) {
    00544     goto L40;
    00545   }
    00546   ++k;
    00547   ifax[k] = 2;
    00548   nn /= 2;
    00549   if (nn == 1) {
    00550     goto L80;
    00551   }
    00552   /*     TEST FOR FACTORS OF 3 */
    00553  L40:
    00554   if (nn % 3 != 0) {
    00555     goto L50;
    00556   }
    00557   ++k;
    00558   ifax[k] = 3;
    00559   nn /= 3;
    00560   if (nn == 1) {
    00561     goto L80;
    00562   }
    00563   goto L40;
    00564   /*     NOW FIND REMAINING FACTORS */
    00565  L50:
    00566   l = 5;
    00567   max = sqrt((double) nn);
    00568   inc = 2;
    00569   /*     INC ALTERNATELY TAKES ON VALUES 2 AND 4 */
    00570  L60:
    00571   if (nn % l != 0) {
    00572     goto L70;
    00573   }
    00574   ++k;
    00575   ifax[k] = l;
    00576   nn /= l;
    00577   if (nn == 1) {
    00578     goto L80;
    00579   }
    00580   goto L60;
    00581  L70:
    00582   if (l > max) {
    00583     goto L75;
    00584   }
    00585   l += inc;
    00586   inc = 6 - inc;
    00587   goto L60;
    00588  L75:
    00589   ++k;
    00590   ifax[k] = nn;
    00591  L80:
    00592   ifax[0] = k;
    00593   /*     IFAX(0) NOW CONTAINS NUMBER OF FACTORS */
    00594   return;
    00595 } /* fact */
    00596
    00597 static void cftrig(long n, double *trigs)
    00598 {
    00599   /* Local variables */
    00600   static long i, n2;
    00601   static double angle, del;
    00602
    00603   /* Function Body */
    00604   del = (PI + PI) / (double) (n);
    00605   n2 = n + n;
    00606   for (i = 0; i < n2; i += 2) {
    00607     angle = (double) i * 0.5 * del;
    00608     trigs[i] = cos(angle);
    00609     trigs[i + 1] = sin(angle);
    00610   }
    00611   return;
    00612 } /* cftrig */
    00613
    00614 static void vpassm(double *a, double *b, double *c, double *d,
    00615                    double *trigs, long inc1, long inc2, long inc3, long inc4,
    00616                    long lot, long n, long ifac, long la)
    00617 {
    00618   /* Initialized data */
    00619
    00620   static double sin36 = 0.587785252292473;
    00621   static double cos36 = 0.809016994374947;
    00622   static double sin72 = 0.951056516295154;
    00623   static double cos72 = 0.309016994374947;
    00624   static double sin60 = 0.866025403784437;
    00625
    00626     /* Local variables */
    00627   static long iink, jink, jump, i, j, k, l, m, ibase, jbase;
    00628   static double c1, c2, c3, c4, s1, s2, s3, s4;
    00629   static long ia, ja, ib, jb, kb, ic, jc, kc, id, jd, kd, ie, je, ke,
    00630     la1, ijk, igo;
    00631
    00632
    00633 /*     SUBROUTINE ^VPASSM^ - MULTIPLE VERSION OF ^VPASSA^ */
    00634 /*     PERFORMS ONE PASS THROUGH DATA */
    00635 /*     AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE */
    00636 /*     A IS FIRST REAL INPUT VECTOR */
    00637 /*     B IS FIRST IMAGINARY INPUT VECTOR */
    00638 /*     C IS FIRST REAL OUTPUT VECTOR */
    00639 /*     D IS FIRST IMAGINARY OUTPUT VECTOR */
    00640 /*     TRIGS IS PRECALCULATED TABLE OF SINES \ COSINES */
    00641 /*     INC1 IS ADDRESSING INCREMENT FOR A AND B */
    00642 /*     INC2 IS ADDRESSING INCREMENT FOR C AND D */
    00643 /*     INC3 IS ADDRESSING INCREMENT BETWEEN A^S \ B^S */
    00644 /*     INC4 IS ADDRESSING INCREMENT BETWEEN C^S \ D^S */
    00645 /*     LOT IS THE NUMBER OF VECTORS */
    00646 /*     N IS LENGTH OF VECTORS */
    00647 /*     IFAC IS CURRENT FACTOR OF N */
    00648 /*     LA IS PRODUCT OF PREVIOUS FACTORS */
    00649
    00650     /* Parameter adjustments */
    00651   --d;
    00652   --c;
    00653   --b;
    00654   --a;
    00655
    00656     /* Function Body */
    00657
    00658   m = n / ifac;
    00659   iink = m * inc1;
    00660   jink = la * inc2;
    00661   jump = (ifac - 1) * jink;
    00662   ibase = 0;
    00663   jbase = 0;
    00664   igo = ifac - 1;
    00665   if (igo > 4) {
    00666     return;
    00667   }
    00668   switch ((int)igo) {
    00669   case 1:  goto L10;
    00670   case 2:  goto L50;
    00671   case 3:  goto L90;
    00672   case 4:  goto L130;
    00673   }
    00674
    00675   /*     CODING FOR FACTOR 2 */
    00676
    00677  L10:
    00678   ia = 1;
    00679   ja = 1;
    00680   ib = ia + iink;
    00681   jb = ja + jink;
    00682   for (l = 1; l <= la; ++l) {
    00683     i = ibase;
    00684     j = jbase;
    00685     /* DIR$ IVDEP */
    00686     for (ijk = 1; ijk <= lot; ++ijk) {
    00687       c[ja + j] = a[ia + i] + a[ib + i];
    00688       d[ja + j] = b[ia + i] + b[ib + i];
    00689       c[jb + j] = a[ia + i] - a[ib + i];
    00690       d[jb + j] = b[ia + i] - b[ib + i];
    00691       i += inc3;
    00692       j += inc4;
    00693       /* L15: */
    00694     }
    00695     ibase += inc1;
    00696     jbase += inc2;
    00697     /* L20: */
    00698   }
    00699   if (la == m) {
    00700     return;
    00701   }
    00702   la1 = la + 1;
    00703   jbase += jump;
    00704   for (k = la1; la < 0 ? k >= m : k <= m; k += la) {
    00705     kb = k + k - 2;
    00706     c1 = trigs[kb];
    00707     s1 = trigs[kb + 1];
    00708     for (l = 1; l <= la; ++l) {
    00709       i = ibase;
    00710       j = jbase;
    00711       /* DIR$ IVDEP */
    00712       for (ijk = 1; ijk <= lot; ++ijk) {
    00713         c[ja + j] = a[ia + i] + a[ib + i];
    00714         d[ja + j] = b[ia + i] + b[ib + i];
    00715         c[jb + j] = c1 * (a[ia + i] - a[ib + i]) - s1 *
    00716           (b[ia +i] - b[ib + i]);
    00717         d[jb + j] = s1 * (a[ia + i] - a[ib + i]) + c1 *
    00718           (b[ia +i] - b[ib + i]);
    00719         i += inc3;
    00720         j += inc4;
    00721         /* L25: */
    00722       }
    00723       ibase += inc1;
    00724       jbase += inc2;
    00725       /* L30: */
    00726     }
    00727     jbase += jump;
    00728     /* L40: */
    00729   }
    00730   return;
    00731
    00732 /*     CODING FOR FACTOR 3 */
    00733
    00734  L50:
    00735   ia = 1;
    00736   ja = 1;
    00737   ib = ia + iink;
    00738   jb = ja + jink;
    00739   ic = ib + iink;
    00740   jc = jb + jink;
    00741   for (l = 1; l <= la; ++l) {
    00742     i = ibase;
    00743     j = jbase;
    00744     /* DIR$ IVDEP */
    00745     for (ijk = 1; ijk <= lot; ++ijk) {
    00746       c[ja + j] = a[ia + i] + (a[ib + i] + a[ic + i]);
    00747       d[ja + j] = b[ia + i] + (b[ib + i] + b[ic + i]);
    00748       c[jb + j] = a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 -
    00749         sin60 * (b[ib + i] - b[ic + i]);
    00750       c[jc + j] = a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 +
    00751         sin60 * (b[ib + i] - b[ic + i]);
    00752       d[jb + j] = b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 +
    00753         sin60 * (a[ib + i] - a[ic + i]);
    00754       d[jc + j] = b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 -
    00755         sin60 * (a[ib + i] - a[ic + i]);
    00756       i += inc3;
    00757       j += inc4;
    00758       /* L55: */
    00759     }
    00760     ibase += inc1;
    00761     jbase += inc2;
    00762     /* L60: */
    00763   }
    00764   if (la == m) {
    00765     return;
    00766   }
    00767   la1 = la + 1;
    00768   jbase += jump;
    00769   for (k = la1; la < 0 ? k >= m : k <= m; k += la) {
    00770     kb = k + k - 2;
    00771     kc = kb + kb;
    00772     c1 = trigs[kb];
    00773     s1 = trigs[kb + 1];
    00774     c2 = trigs[kc];
    00775     s2 = trigs[kc + 1];
    00776     for (l = 1; l <= la; ++l) {
    00777       i = ibase;
    00778       j = jbase;
    00779       /* DIR$ IVDEP */
    00780       for (ijk = 1; ijk <= lot; ++ijk) {
    00781         c[ja + j] = a[ia + i] + (a[ib + i] + a[ic + i]);
    00782         d[ja + j] = b[ia + i] + (b[ib + i] + b[ic + i]);
    00783         c[jb + j] = c1 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 -
    00784                           sin60 * (b[ib + i] - b[ic + i])) -
    00785           s1 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 +
    00786                 sin60 * (a[ib + i] - a[ic + i]));
    00787         d[jb + j] = s1 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 -
    00788                           sin60 * (b[ib + i] - b[ic + i])) +
    00789           c1 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 +
    00790                 sin60 * (a[ib + i] - a[ic + i]));
    00791         c[jc + j] = c2 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 +
    00792                           sin60 * (b[ib + i] - b[ic + i])) -
    00793           s2 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 -
    00794                 sin60 * (a[ib + i] - a[ic + i]));
    00795         d[jc + j] = s2 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 +
    00796                           sin60 * (b[ib + i] - b[ic + i])) +
    00797           c2 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 -
    00798                 sin60 * (a[ib + i] - a[ic + i]));
    00799         i += inc3;
    00800         j += inc4;
    00801         /* L65: */
    00802       }
    00803       ibase += inc1;
    00804       jbase += inc2;
    00805       /* L70: */
    00806     }
    00807     jbase += jump;
    00808     /* L80: */
    00809   }
    00810   return;
    00811
    00812 /*     CODING FOR FACTOR 4 */
    00813
    00814  L90:
    00815   ia = 1;
    00816   ja = 1;
    00817   ib = ia + iink;
    00818   jb = ja + jink;
    00819   ic = ib + iink;
    00820   jc = jb + jink;
    00821   id = ic + iink;
    00822   jd = jc + jink;
    00823   for (l = 1; l <= la; ++l) {
    00824     i = ibase;
    00825     j = jbase;
    00826     /* DIR$ IVDEP */
    00827     for (ijk = 1; ijk <= lot; ++ijk) {
    00828       c[ja + j] = a[ia + i] + a[ic + i] + (a[ib + i] + a[id + i]);
    00829       c[jc + j] = a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i]);
    00830       d[ja + j] = b[ia + i] + b[ic + i] + (b[ib + i] + b[id + i]);
    00831       d[jc + j] = b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i]);
    00832       c[jb + j] = a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i]);
    00833       c[jd + j] = a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i]);
    00834       d[jb + j] = b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i]);
    00835       d[jd + j] = b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i]);
    00836       i += inc3;
    00837       j += inc4;
    00838       /* L95: */
    00839     }
    00840     ibase += inc1;
    00841     jbase += inc2;
    00842     /* L100: */
    00843   }
    00844   if (la == m) {
    00845     return;
    00846   }
    00847   la1 = la + 1;
    00848   jbase += jump;
    00849   for (k = la1; la < 0 ? k >= m : k <= m; k += la) {
    00850     kb = k + k - 2;
    00851     kc = kb + kb;
    00852     kd = kc + kb;
    00853     c1 = trigs[kb];
    00854     s1 = trigs[kb + 1];
    00855     c2 = trigs[kc];
    00856     s2 = trigs[kc + 1];
    00857     c3 = trigs[kd];
    00858     s3 = trigs[kd + 1];
    00859     for (l = 1; l <= la; ++l) {
    00860       i = ibase;
    00861       j = jbase;
    00862       /* DIR$ IVDEP */
    00863       for (ijk = 1; ijk <= lot; ++ijk) {
    00864         c[ja + j] = a[ia + i] + a[ic + i] + (a[ib + i] + a[id + i]);
    00865         d[ja + j] = b[ia + i] + b[ic + i] + (b[ib + i] + b[id + i]);
    00866         c[jc + j] = c2 * (a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i])) -
    00867           s2 * (b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i]));
    00868         d[jc + j] = s2 * (a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i])) +
    00869           c2 * (b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i]));
    00870         c[jb + j] = c1 * (a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i])) -
    00871           s1 * (b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i]));
    00872         d[jb + j] = s1 * (a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i])) +
    00873           c1 * (b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i]));
    00874         c[jd + j] = c3 * (a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i])) -
    00875           s3 * (b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i]));
    00876         d[jd + j] = s3 * (a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i])) +
    00877           c3 * (b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i]));
    00878         i += inc3;
    00879         j += inc4;
    00880         /* L105: */
    00881       }
    00882       ibase += inc1;
    00883       jbase += inc2;
    00884       /* L110: */
    00885     }
    00886     jbase += jump;
    00887     /* L120: */
    00888   }
    00889   return;
    00890
    00891 /*     CODING FOR FACTOR 5 */
    00892
    00893  L130:
    00894   ia = 1;
    00895   ja = 1;
    00896   ib = ia + iink;
    00897   jb = ja + jink;
    00898   ic = ib + iink;
    00899   jc = jb + jink;
    00900   id = ic + iink;
    00901   jd = jc + jink;
    00902   ie = id + iink;
    00903   je = jd + jink;
    00904   for (l = 1; l <= la; ++l) {
    00905     i = ibase;
    00906     j = jbase;
    00907     /* DIR$ IVDEP */
    00908     for (ijk = 1; ijk <= lot; ++ijk) {
    00909       c[ja + j] = a[ia + i] + (a[ib + i] + a[ie + i]) +
    00910         (a[ic + i] + a[id + i]);
    00911       d[ja + j] = b[ia + i] + (b[ib + i] + b[ie + i]) +
    00912         (b[ic + i] + b[id + i]);
    00913       c[jb + j] = a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    00914         cos36 * (a[ic + i] + a[id + i]) - (sin72 * (b[ib + i] - b[ie + i]) +
    00915                                            sin36 * (b[ic + i] - b[id + i]));
    00916       c[je + j] = a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    00917         cos36 * (a[ic + i] + a[id + i]) + (sin72 * (b[ib + i] - b[ie + i]) +
    00918                                            sin36 * (b[ic + i] - b[id + i]));
    00919       d[jb + j] = b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    00920         cos36 * (b[ic + i] + b[id + i]) + (sin72 * (a[ib + i] - a[ie + i]) +
    00921                                            sin36 * (a[ic + i] - a[id + i]));
    00922       d[je + j] = b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    00923         cos36 * (b[ic + i] + b[id + i]) - (sin72 * (a[ib + i] - a[ie + i]) +
    00924                                            sin36 * (a[ic + i] - a[id + i]));
    00925       c[jc + j] = a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    00926         cos72 * (a[ic + i] + a[id + i]) - (sin36 * (b[ib + i] - b[ie + i]) -
    00927                                            sin72 * (b[ic + i] - b[id + i]));
    00928       c[jd + j] = a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    00929         cos72 * (a[ic + i] + a[id + i]) + (sin36 * (b[ib + i] - b[ie + i]) -
    00930                                            sin72 * (b[ic + i] - b[id + i]));
    00931       d[jc + j] = b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    00932         cos72 * (b[ic + i] + b[id + i]) + (sin36 * (a[ib + i] - a[ie + i]) -
    00933                                            sin72 * (a[ic + i] - a[id + i]));
    00934       d[jd + j] = b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    00935         cos72 * (b[ic + i] + b[id + i]) - (sin36 * (a[ib + i] - a[ie + i]) -
    00936                                            sin72 * (a[ic + i] - a[id + i]));
    00937       i += inc3;
    00938       j += inc4;
    00939       /* L135: */
    00940     }
    00941     ibase += inc1;
    00942     jbase += inc2;
    00943     /* L140: */
    00944   }
    00945   if (la == m) {
    00946     return;
    00947   }
    00948   la1 = la + 1;
    00949   jbase += jump;
    00950   for (k = la1; la < 0 ? k >= m : k <= m; k += la) {
    00951     kb = k + k - 2;
    00952     kc = kb + kb;
    00953     kd = kc + kb;
    00954     ke = kd + kb;
    00955     c1 = trigs[kb];
    00956     s1 = trigs[kb + 1];
    00957     c2 = trigs[kc];
    00958     s2 = trigs[kc + 1];
    00959     c3 = trigs[kd];
    00960     s3 = trigs[kd + 1];
    00961     c4 = trigs[ke];
    00962     s4 = trigs[ke + 1];
    00963     for (l = 1; l <= la; ++l) {
    00964       i = ibase;
    00965       j = jbase;
    00966       /* DIR$ IVDEP */
    00967       for (ijk = 1; ijk <= lot; ++ijk) {
    00968         c[ja + j] = a[ia + i] + (a[ib + i] + a[ie + i]) +
    00969           (a[ic + i] + a[id + i]);
    00970         d[ja + j] = b[ia + i] + (b[ib + i] + b[ie + i]) +
    00971           (b[ic + i] + b[id + i]);
    00972         c[jb + j] =
    00973           c1 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    00974                 cos36 * (a[ic + i] + a[id + i]) -
    00975                 (sin72 * (b[ib + i] - b[ie + i]) +
    00976                  sin36 * (b[ic + i] - b[id + i]))) -
    00977           s1 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    00978                 cos36 * (b[ic + i] + b[id + i]) +
    00979                 (sin72 * (a[ib + i] - a[ie + i]) +
    00980                  sin36 * (a[ic + i] - a[id + i])));
    00981         d[jb + j] =
    00982           s1 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    00983                 cos36 * (a[ic + i] + a[id + i]) -
    00984                 (sin72 * (b[ib + i] - b[ie + i]) +
    00985                  sin36 * (b[ic + i] - b[id + i]))) +
    00986           c1 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    00987                 cos36 * (b[ic + i] + b[id + i]) +
    00988                 (sin72 * (a[ib + i] - a[ie + i]) +
    00989                  sin36 * (a[ic + i] - a[id + i])));
    00990         c[je + j] =
    00991           c4 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    00992                 cos36 * (a[ic + i] + a[id + i]) +
    00993                 (sin72 * (b[ib + i] - b[ie + i]) +
    00994                  sin36 * (b[ic + i] - b[id + i]))) -
    00995           s4 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    00996                 cos36 * (b[ic + i] + b[id + i]) -
    00997                 (sin72 * (a[ib + i] - a[ie + i]) +
    00998                  sin36 * (a[ic + i] - a[id + i])));
    00999         d[je + j] =
    01000           s4 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) -
    01001                 cos36 * (a[ic + i] + a[id + i]) +
    01002                 (sin72 * (b[ib + i] - b[ie + i]) +
    01003                  sin36 * (b[ic + i] - b[id + i]))) +
    01004           c4 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) -
    01005                 cos36 * (b[ic + i] + b[id + i]) -
    01006                 (sin72 * (a[ib + i] - a[ie + i]) +
    01007                  sin36 * (a[ic + i] - a[id + i])));
    01008         c[jc + j] =
    01009           c2 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    01010                 cos72 * (a[ic + i] + a[id + i]) -
    01011                 (sin36 * (b[ib + i] - b[ie + i]) -
    01012                  sin72 * (b[ic + i] - b[id + i]))) -
    01013           s2 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    01014                 cos72 * (b[ic + i] + b[id + i]) +
    01015                 (sin36 * (a[ib + i] - a[ie + i]) -
    01016                  sin72 * (a[ic + i] - a[id + i])));
    01017         d[jc + j] =
    01018           s2 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    01019                 cos72 * (a[ic + i] + a[id + i]) -
    01020                 (sin36 * (b[ib + i] - b[ie + i]) -
    01021                  sin72 * (b[ic + i] - b[id + i]))) +
    01022           c2 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    01023                 cos72 * (b[ic + i] + b[id + i]) +
    01024                 (sin36 * (a[ib + i] - a[ie + i]) -
    01025                  sin72 * (a[ic + i] - a[id + i])));
    01026         c[jd + j] =
    01027           c3 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    01028                 cos72 * (a[ic + i] + a[id + i]) +
    01029                 (sin36 * (b[ib + i] - b[ie + i]) -
    01030                  sin72 * (b[ic +i] - b[id + i]))) -
    01031           s3 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    01032                 cos72 * (b[ic + i] + b[id + i]) -
    01033                 (sin36 * (a[ib + i] - a[ie + i]) -
    01034                  sin72 * (a[ic + i] - a[id + i])));
    01035         d[jd + j] =
    01036           s3 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) +
    01037                 cos72 * (a[ic + i] + a[id + i]) +
    01038                 (sin36 * (b[ib + i] - b[ie + i]) -
    01039                  sin72 * (b[ic + i] - b[id + i]))) +
    01040           c3 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) +
    01041                 cos72 * (b[ic + i] + b[id + i]) -
    01042                 (sin36 * (a[ib + i] - a[ie + i]) -
    01043                  sin72 * (a[ic + i] - a[id + i])));
    01044         i += inc3;
    01045         j += inc4;
    01046         /* L145: */
    01047       }
    01048       ibase += inc1;
    01049       jbase += inc2;
    01050       /* L150: */
    01051     }
    01052     jbase += jump;
    01053     /* L160: */
    01054   }
    01055   return;
    01056 } /* vpassm */
    01057
    01058 void fft99(double *a, double *work, double *trigs, long *ifax, long inc,
    01059            long jump, long n, long lot, long isign)
    01060 {
    01061   /* Local variables */
    01062   static long nfax, i, j, k, l, m, ibase, jbase;
    01063   static long ia, ib, la, nh, nx;
    01064   static long igo, ink;
    01065
    01066
    01067 /* PURPOSE      PERFORMS MULTIPLE FAST FOURIER TRANSFORMS.  THIS PACKAGE
    01068  *              WILL PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
    01069  *              PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
    01070  *              TRANSFORMS, I.E.  GIVEN A SET OF REAL DATA VECTORS, THE
    01071  *              PACKAGE RETURNS A SET OF #HALF-COMPLEX# FOURIER
    01072  *              COEFFICIENT VECTORS, OR VICE VERSA.  THE LENGTH OF THE
    01073  *              TRANSFORMS MUST BE AN EVEN NUMBER GREATER THAN 4 THAT HAS
    01074  *              NO OTHER FACTORS EXCEPT POSSIBLY POWERS OF 2, 3, AND 5.
    01075  *
    01076  *              THE PACKAGE FFT99F CONTAINS SEVERAL USER-LEVEL ROUTINES:
    01077  *
    01078  *             fftfax
    01079  *                AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE
    01080  *                BEFORE A SEQUENCE OF CALLS TO THE FFT ROUTINES
    01081  *                (PROVIDED THAT N IS NOT CHANGED).
    01082  *
    01083  *             fft99 AND fft991
    01084  *                TWO FFT ROUTINES THAT RETURN SLIGHTLY DIFFERENT
    01085  *                ARRANGEMENTS OF THE DATA IN GRIDPOINT SPACE.
    01086  *
    01087  * USAGE        LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 1,
    01088  *              Q .GE. 0, AND R .GE. 0.  THEN A TYPICAL SEQUENCE OF
    01089  *              CALLS TO TRANSFORM A GIVEN SET OF REAL VECTORS OF LENGTH
    01090  *              N TO A SET OF #HALF-COMPLEX# FOURIER COEFFICIENT VECTORS
    01091  *              OF LENGTH N IS
    01092  *
    01093  *                   long ifax[13];
    01094  *                   double trigs[3*n/2+1],a[m*(n+2)], work[m*(n+1)];
    01095  *
    01096  *                    fftfax (n, ifax, trigs);
    01097  *                    fft99 (a,work,trigs,ifax,inc,jump,n,m,isign);
    01098  *
    01099  *              SEE THE INDIVIDUAL WRITE-UPS FOR FFTFAX, FFT99, AND
    01100  *              FFT991 BELOW, FOR A DETAILED DESCRIPTION OF THE
    01101  *              ARGUMENTS.
    01102  *
    01103  * HISTORY      THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN
    01104  *              NOVEMBER, 1978.  IT WAS MODIFIED, DOCUMENTED, AND TESTED
    01105  *              FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.
    01106  *
    01107  * -----------------------------------------------------------------------
    01108  *
    01109  *  fftfax (n,ifax,trigs)
    01110  *
    01111  * PURPOSE      A SET-UP ROUTINE FOR FFT99 AND FFT991.  IT NEED ONLY BE
    01112  *              CALLED ONCE BEFORE A SEQUENCE OF CALLS TO THE FFT
    01113  *              ROUTINES (PROVIDED THAT N IS NOT CHANGED).
    01114  *
    01115  * ARGUMENT     IFAX(13),TRIGS(3*N/2+1)
    01116  * DIMENSIONS
    01117  *
    01118  * ARGUMENTS
    01119  *
    01120  * ON INPUT     N
    01121  *               AN EVEN NUMBER GREATER THAN 4 THAT HAS NO PRIME FACTOR
    01122  *               GREATER THAN 5.  N IS THE LENGTH OF THE TRANSFORMS (SEE
    01123  *               THE DOCUMENTATION FOR FFT99 AND FFT991 FOR THE
    01124  *               DEFINITIONS OF THE TRANSFORMS).
    01125  *
    01126  *              IFAX
    01127  *               AN INTEGER ARRAY.  THE NUMBER OF ELEMENTS ACTUALLY USED
    01128  *               WILL DEPEND ON THE FACTORIZATION OF N.  DIMENSIONING
    01129  *               IFAX FOR 13 SUFFICES FOR ALL N LESS THAN A MILLION.
    01130  *
    01131  *              TRIGS
    01132  *               A FLOATING POINT ARRAY OF DIMENSION 3*N/2 IF N/2 IS
    01133  *               EVEN, OR 3*N/2+1 IF N/2 IS ODD.
    01134  *
    01135  * ON OUTPUT    IFAX
    01136  *               CONTAINS THE FACTORIZATION OF N/2.  IFAX(1) IS THE
    01137  *               NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED
    01138  *               IN IFAX(2),IFAX(3),...  IF FFTFAX IS CALLED WITH N ODD,
    01139  *               OR IF N HAS ANY PRIME FACTORS GREATER THAN 5, IFAX(1)
    01140  *               IS SET TO -99.
    01141  *
    01142  *              TRIGS
    01143  *               AN ARRAY OF TRIGNOMENTRIC FUNCTION VALUES SUBSEQUENTLY
    01144  *               USED BY THE FFT ROUTINES.
    01145  *
    01146  * -----------------------------------------------------------------------
    01147  *
    01148  *  fft991 (a,work,trigs,ifax,inc,jump,n,m,isign)
    01149  *                       AND
    01150  *  fft99 (a,work,trigs,ifax,inc,jump,n,m,isign)
    01151  *
    01152  * PURPOSE      PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
    01153  *              PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
    01154  *              TRANSFORMS, USING ORDINARY SPATIAL ORDER OF GRIDPOINT
    01155  *              VALUES (FFT991) OR EXPLICIT CYCLIC CONTINUITY IN THE
    01156  *              GRIDPOINT VALUES (FFT99).  GIVEN A SET
    01157  *              OF REAL DATA VECTORS, THE PACKAGE RETURNS A SET OF
    01158  *              #HALF-COMPLEX# FOURIER COEFFICIENT VECTORS, OR VICE
    01159  *              VERSA.  THE LENGTH OF THE TRANSFORMS MUST BE AN EVEN
    01160  *              NUMBER THAT HAS NO OTHER FACTORS EXCEPT POSSIBLY POWERS
    01161  *              OF 2, 3, AND 5.  THESE VERSION OF FFT991 AND FFT99 ARE
    01162  *              OPTIMIZED FOR USE ON THE CRAY-1.
    01163  *
    01164  * ARGUMENT     A(M*(N+2)), WORK(M*(N+1)), TRIGS(3*N/2+1), IFAX(13)
    01165  * DIMENSIONS
    01166  *
    01167  * ARGUMENTS
    01168  *
    01169  * ON INPUT     A
    01170  *               AN ARRAY OF LENGTH M*(N+2) CONTAINING THE INPUT DATA
    01171  *               OR COEFFICIENT VECTORS.  THIS ARRAY IS OVERWRITTEN BY
    01172  *               THE RESULTS.
    01173  *
    01174  *              WORK
    01175  *               A WORK ARRAY OF DIMENSION M*(N+1)
    01176  *
    01177  *              TRIGS
    01178  *               AN ARRAY SET UP BY FFTFAX, WHICH MUST BE CALLED FIRST.
    01179  *
    01180  *              IFAX
    01181  *               AN ARRAY SET UP BY FFTFAX, WHICH MUST BE CALLED FIRST.
    01182  *
    01183  *              INC
    01184  *               THE INCREMENT (IN WORDS) BETWEEN SUCCESSIVE ELEMENTS OF
    01185  *               EACH DATA OR COEFFICIENT VECTOR (E.G.  INC=1 FOR
    01186  *               CONSECUTIVELY STORED DATA).
    01187  *
    01188  *              JUMP
    01189  *               THE INCREMENT (IN WORDS) BETWEEN THE FIRST ELEMENTS OF
    01190  *               SUCCESSIVE DATA OR COEFFICIENT VECTORS.  ON THE CRAY-1,
    01191  *               TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8
    01192  *               (TO AVOID MEMORY BANK CONFLICTS).  FOR CLARIFICATION OF
    01193  *               INC AND JUMP, SEE THE EXAMPLES BELOW.
    01194  *
    01195  *              N
    01196  *               THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF
    01197  *               TRANSFORMS, BELOW).
    01198  *
    01199  *              M
    01200  *               THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY.
    01201  *
    01202  *              ISIGN
    01203  *               = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO
    01204  *                    GRIDPOINT VALUES.
    01205  *               = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER
    01206  *                    COEFFICIENTS.
    01207  *
    01208  * ON OUTPUT    A
    01209  *               IF ISIGN = +1, AND M COEFFICIENT VECTORS ARE SUPPLIED
    01210  *               EACH CONTAINING THE SEQUENCE:
    01211  *
    01212  *               A(0),B(0),A(1),B(1),...,A(N/2),B(N/2)  (N+2 VALUES)
    01213  *
    01214  *               THEN THE RESULT CONSISTS OF M DATA VECTORS EACH
    01215  *               CONTAINING THE CORRESPONDING N+2 GRIDPOINT VALUES:
    01216  *
    01217  *               FOR FFT991, X(0), X(1), X(2),...,X(N-1),0,0.
    01218  *               FOR FFT99, X(N-1),X(0),X(1),X(2),...,X(N-1),X(0).
    01219  *                   (EXPLICIT CYCLIC CONTINUITY)
    01220  *
    01221  *               WHEN ISIGN = +1, THE TRANSFORM IS DEFINED BY:
    01222  *                 X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
    01223  *                 WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
    01224  *                 AND I=SQRT (-1)
    01225  *
    01226  *               IF ISIGN = -1, AND M DATA VECTORS ARE SUPPLIED EACH
    01227  *               CONTAINING A SEQUENCE OF GRIDPOINT VALUES X(J) AS
    01228  *               DEFINED ABOVE, THEN THE RESULT CONSISTS OF M VECTORS
    01229  *               EACH CONTAINING THE CORRESPONDING FOURIER COFFICIENTS
    01230  *               A(K), B(K), 0 .LE. K .LE N/2.
    01231  *
    01232  *               WHEN ISIGN = -1, THE INVERSE TRANSFORM IS DEFINED BY:
    01233  *                 C(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*EXP(-2*I*J*K*PI/N))
    01234  *                 WHERE C(K)=A(K)+I*B(K) AND I=SQRT(-1)
    01235  *
    01236  *               A CALL WITH ISIGN=+1 FOLLOWED BY A CALL WITH ISIGN=-1
    01237  *               (OR VICE VERSA) RETURNS THE ORIGINAL DATA.
    01238  *
    01239  *               NOTE: THE FACT THAT THE GRIDPOINT VALUES X(J) ARE REAL
    01240  *               IMPLIES THAT B(0)=B(N/2)=0.  FOR A CALL WITH ISIGN=+1,
    01241  *               IT IS NOT ACTUALLY NECESSARY TO SUPPLY THESE ZEROS.
    01242  *
    01243  * EXAMPLES      GIVEN 19 DATA VECTORS EACH OF LENGTH 64 (+2 FOR EXPLICIT
    01244  *               CYCLIC CONTINUITY), COMPUTE THE CORRESPONDING VECTORS OF
    01245  *               FOURIER COEFFICIENTS.  THE DATA MAY, FOR EXAMPLE, BE
    01246  *               ARRANGED LIKE THIS:
    01247  *
    01248  * FIRST DATA   A(1)=    . . .                A(66)=             A(70)
    01249  * VECTOR       X(63) X(0) X(1) X(2) ... X(63) X(0)  (4 EMPTY LOCATIONS)
    01250  *
    01251  * SECOND DATA  A(71)=   . . .                                  A(140)
    01252  * VECTOR       X(63) X(0) X(1) X(2) ... X(63) X(0)  (4 EMPTY LOCATIONS)
    01253  *
    01254  *               AND SO ON.  HERE INC=1, JUMP=70, N=64, M=19, ISIGN=-1,
    01255  *               AND FFT99 SHOULD BE USED (BECAUSE OF THE EXPLICIT CYCLIC
    01256  *               CONTINUITY).
    01257  *
    01258  *               ALTERNATIVELY THE DATA MAY BE ARRANGED LIKE THIS:
    01259  *
    01260  *                FIRST         SECOND                          LAST
    01261  *                DATA          DATA                            DATA
    01262  *                VECTOR        VECTOR                          VECTOR
    01263  *
    01264  *                 A(1)=         A(2)=                           A(19)=
    01265  *
    01266  *                 X(63)         X(63)       . . .               X(63)
    01267  *        A(20)=   X(0)          X(0)        . . .               X(0)
    01268  *        A(39)=   X(1)          X(1)        . . .               X(1)
    01269  *                  .             .                               .
    01270  *                  .             .                               .
    01271  *                  .             .                               .
    01272  *
    01273  *               IN WHICH CASE WE HAVE INC=19, JUMP=1, AND THE REMAINING
    01274  *               PARAMETERS ARE THE SAME AS BEFORE.  IN EITHER CASE, EACH
    01275  *               COEFFICIENT VECTOR OVERWRITES THE CORRESPONDING INPUT
    01276  *               DATA VECTOR.
    01277  *
    01278  * -----------------------------------------------------------------------
    01279  *
    01280  *     SUBROUTINE ^FFT99^ - MULTIPLE FAST REAL PERIODIC TRANSFORM
    01281  *     CORRESPONDING TO OLD SCALAR ROUTINE FFT9
    01282  *     PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM
    01283  *     IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12
    01284  *     (1970), 315-337)
    01285  *
    01286  *     A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
    01287  *     WORK IS AN AREA OF SIZE (N+1)*LOT
    01288  *     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
    01289  *     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
    01290  *     INC IS THE INCREMENT WITHIN EACH DATA #VECTOR#
    01291  *         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
    01292  *     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
    01293  *     N IS THE LENGTH OF THE DATA VECTORS
    01294  *     LOT IS THE NUMBER OF DATA VECTORS
    01295  *     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
    01296  *           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
    01297  *
    01298  *     ORDERING OF COEFFICIENTS:
    01299  *         A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
    01300  *         WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
    01301  *
    01302  *     ORDERING OF DATA:
    01303  *         X(N-1),X(0),X(1),X(2),...,X(N),X(0)
    01304  *         I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED
    01305  *
    01306  *     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
    01307  *     PARALLEL
    01308  *
    01309  *     *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
    01310  *
    01311  *     DEFINITION OF TRANSFORMS:
    01312  *     -------------------------
    01313  *
    01314  *     ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
    01315  *         WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
    01316  *
    01317  *     ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
    01318  *               B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
    01319  */
    01320
    01321     /* Parameter adjustments */
    01322   --work;
    01323   --a;
    01324
    01325     /* Function Body */
    01326   nfax = ifax[0];
    01327   nx = n + 1;
    01328   nh = n / 2;
    01329   ink = inc + inc;
    01330   if (isign == 1) {
    01331     goto L30;
    01332   }
    01333
    01334   /*     IF NECESSARY, TRANSFER DATA TO WORK AREA */
    01335   igo = 50;
    01336   if (nfax % 2 == 1) {
    01337     goto L40;
    01338   }
    01339   ibase = inc + 1;
    01340   jbase = 1;
    01341   for (l = 1; l <= lot; ++l) {
    01342     i = ibase;
    01343     j = jbase;
    01344     /* DIR$ IVDEP */
    01345     for (m = 1; m <= n; ++m) {
    01346       work[j] = a[i];
    01347       i += inc;
    01348       ++j;
    01349       /* L10: */
    01350     }
    01351     ibase += jump;
    01352     jbase += nx;
    01353     /* L20: */
    01354   }
    01355
    01356   igo = 60;
    01357   goto L40;
    01358
    01359   /*     PREPROCESSING (ISIGN=+1) */
    01360   /*     ------------------------ */
    01361
    01362  L30:
    01363   fft99a(&a[1], &work[1], trigs, inc, jump, n, lot);
    01364   igo = 60;
    01365
    01366 /*     COMPLEX TRANSFORM */
    01367 /*     ----------------- */
    01368
    01369  L40:
    01370   ia = inc + 1;
    01371   la = 1;
    01372   for (k = 1; k <= nfax; ++k) {
    01373     if (igo == 60) {
    01374       goto L60;
    01375     }
    01376     /* L50: */
    01377     vpassm(&a[ia], &a[ia + inc], &work[1], &work[2], trigs, ink, 2, jump, nx,
    01378            lot, nh, ifax[k], la);
    01379     igo = 60;
    01380     goto L70;
    01381   L60:
    01382     vpassm(&work[1], &work[2], &a[ia], &a[ia + inc], trigs, 2, ink, nx, jump,
    01383            lot, nh, ifax[k], la);
    01384     igo = 50;
    01385   L70:
    01386     la *= ifax[k];
    01387     /* L80: */
    01388   }
    01389
    01390   if (isign == -1) {
    01391     goto L130;
    01392   }
    01393
    01394   /*     IF NECESSARY, TRANSFER DATA FROM WORK AREA */
    01395   if (nfax % 2 == 1) {
    01396     goto L110;
    01397   }
    01398   ibase = 1;
    01399   jbase = ia;
    01400   for (l = 1; l <= lot; ++l) {
    01401     i = ibase;
    01402     j = jbase;
    01403     /* DIR$ IVDEP */
    01404     for (m = 1; m <= n; ++m) {
    01405       a[j] = work[i];
    01406       ++i;
    01407       j += inc;
    01408       /* L90: */
    01409     }
    01410     ibase += nx;
    01411     jbase += jump;
    01412     /* L100: */
    01413   }
    01414
    01415   /*     FILL IN CYCLIC BOUNDARY POINTS */
    01416  L110:
    01417   ia = 1;
    01418   ib = n * inc + 1;
    01419   /* DIR$ IVDEP */
    01420   for (l = 1; l <= lot; ++l) {
    01421     a[ia] = a[ib];
    01422     a[ib + inc] = a[ia + inc];
    01423     ia += jump;
    01424     ib += jump;
    01425     /* L120: */
    01426   }
    01427   goto L140;
    01428
    01429 /*     POSTPROCESSING (ISIGN=-1): */
    01430 /*     -------------------------- */
    01431
    01432  L130:
    01433   fft99b(&work[1], &a[1], trigs, inc, jump, n, lot);
    01434
    01435  L140:
    01436   return;
    01437 } /* fft99 */
    01438
    01439 static void fft99a(double *a, double *work, double *trigs, long inc, long jump,
    01440                    long n, long lot)
    01441 {
    01442   /* Local variables */
    01443   static double c;
    01444   static long k, l;
    01445   static double s;
    01446   static long ia, ib, ja, jb, iabase, ibbase, nh, jabase, jbbase, nx,
    01447     ink;
    01448
    01449
    01450 /*     SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1 */
    01451 /*     (SPECTRAL TO GRIDPOINT TRANSFORM) */
    01452
    01453     /* Parameter adjustments */
    01454   --work;
    01455   --a;
    01456
    01457     /* Function Body */
    01458   nh = n / 2;
    01459   nx = n + 1;
    01460   ink = inc + inc;
    01461
    01462   /*     A(0) AND A(N/2) */
    01463   ia = 1;
    01464   ib = n * inc + 1;
    01465   ja = 1;
    01466   jb = 2;
    01467   /* DIR$ IVDEP */
    01468   for (l = 1; l <= lot; ++l) {
    01469     work[ja] = a[ia] + a[ib];
    01470     work[jb] = a[ia] - a[ib];
    01471     ia += jump;
    01472     ib += jump;
    01473     ja += nx;
    01474     jb += nx;
    01475     /* L10: */
    01476   }
    01477
    01478   /*     REMAINING WAVENUMBERS */
    01479   iabase = (inc << 1) + 1;
    01480   ibbase = (n - 2) * inc + 1;
    01481   jabase = 3;
    01482   jbbase = n - 1;
    01483
    01484   for (k = 3; k <= nh; k += 2) {
    01485     ia = iabase;
    01486     ib = ibbase;
    01487     ja = jabase;
    01488     jb = jbbase;
    01489     c = trigs[n + k - 1];
    01490     s = trigs[n + k];
    01491     /* DIR$ IVDEP */
    01492     for (l = 1; l <= lot; ++l) {
    01493       work[ja] = a[ia] + a[ib] - (s * (a[ia] - a[ib]) +
    01494                                   c * (a[ia + inc] + a[ib + inc]));
    01495       work[jb] = a[ia] + a[ib] + (s * (a[ia] - a[ib]) +
    01496                                   c * (a[ia + inc] + a[ib + inc]));
    01497       work[ja + 1] = c * (a[ia] - a[ib]) -
    01498         s * (a[ia + inc] + a[ib + inc]) + (a[ia + inc] - a[ib + inc]);
    01499       work[jb + 1] = c * (a[ia] - a[ib]) -
    01500         s * (a[ia + inc] + a[ib + inc]) - (a[ia + inc] - a[ib + inc]);
    01501       ia += jump;
    01502       ib += jump;
    01503       ja += nx;
    01504       jb += nx;
    01505       /* L20: */
    01506     }
    01507     iabase += ink;
    01508     ibbase -= ink;
    01509     jabase += 2;
    01510     jbbase += -2;
    01511     /* L30: */
    01512   }
    01513
    01514   if (iabase != ibbase) {
    01515     goto L50;
    01516   }
    01517   /*     WAVENUMBER N/4 (IF IT EXISTS) */
    01518   ia = iabase;
    01519   ja = jabase;
    01520   /* DIR$ IVDEP */
    01521   for (l = 1; l <= lot; ++l) {
    01522     work[ja] = a[ia] * 2.0;
    01523     work[ja + 1] = a[ia + inc] * -2.0;
    01524     ia += jump;
    01525     ja += nx;
    01526     /* L40: */
    01527   }
    01528
    01529  L50:
    01530   return;
    01531 } /* fft99a */
    01532
    01533 static void fft99b(double *work, double *a, double *trigs, long inc, long jump,
    01534                    long n, long lot)
    01535 {
    01536   /* Local variables */
    01537   static double c;
    01538   static long k, l;
    01539   static double s, scale;
    01540   static long ia, ib, ja, jb, iabase, ibbase, nh, jabase, jbbase, nx,
    01541     ink;
    01542
    01543
    01544 /*     SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1 */
    01545 /*     (GRIDPOINT TO SPECTRAL TRANSFORM) */
    01546
    01547     /* Parameter adjustments */
    01548   --a;
    01549   --work;
    01550
    01551     /* Function Body */
    01552   nh = n / 2;
    01553   nx = n + 1;
    01554   ink = inc + inc;
    01555
    01556   /*     A(0) AND A(N/2) */
    01557   scale = 1.0 / (double) n;
    01558   ia = 1;
    01559   ib = 2;
    01560   ja = 1;
    01561   jb = n * inc + 1;
    01562   /* DIR$ IVDEP */
    01563   for (l = 1; l <= lot; ++l) {
    01564     a[ja] = scale * (work[ia] + work[ib]);
    01565     a[jb] = scale * (work[ia] - work[ib]);
    01566     a[ja + inc] = 0.0;
    01567     a[jb + inc] = 0.0;
    01568     ia += nx;
    01569     ib += nx;
    01570     ja += jump;
    01571     jb += jump;
    01572     /* L10: */
    01573   }
    01574
    01575   /*     REMAINING WAVENUMBERS */
    01576   scale *= 0.5;
    01577   iabase = 3;
    01578   ibbase = n - 1;
    01579   jabase = (inc << 1) + 1;
    01580   jbbase = (n - 2) * inc + 1;
    01581
    01582   for (k = 3; k <= nh; k += 2) {
    01583     ia = iabase;
    01584     ib = ibbase;
    01585     ja = jabase;
    01586     jb = jbbase;
    01587     c = trigs[n + k - 1];
    01588     s = trigs[n + k];
    01589     /* DIR$ IVDEP */
    01590     for (l = 1; l <= lot; ++l) {
    01591       a[ja] = scale * (work[ia] + work[ib] +
    01592                        (c * (work[ia + 1] + work[ib + 1]) +
    01593                         s * (work[ia] - work[ib])));
    01594       a[jb] = scale * (work[ia] + work[ib] -
    01595                        (c * (work[ia + 1] + work[ib + 1]) +
    01596                         s * (work[ia] - work[ib])));
    01597       a[ja + inc] = scale * (c * (work[ia] - work[ib]) -
    01598                              s * (work[ia + 1] + work[ib + 1]) +
    01599                              (work[ib + 1] - work[ia + 1]));
    01600       a[jb + inc] = scale * (c * (work[ia] - work[ib]) -
    01601                              s * (work[ia + 1] + work[ib + 1]) -
    01602                              (work[ib + 1] - work[ia + 1]));
    01603       ia += nx;
    01604       ib += nx;
    01605       ja += jump;
    01606       jb += jump;
    01607       /* L20: */
    01608     }
    01609     iabase += 2;
    01610     ibbase += -2;
    01611     jabase += ink;
    01612     jbbase -= ink;
    01613     /* L30: */
    01614   }
    01615
    01616   if (iabase != ibbase) {
    01617     goto L50;
    01618   }
    01619   /*     WAVENUMBER N/4 (IF IT EXISTS) */
    01620   ia = iabase;
    01621   ja = jabase;
    01622   scale *= 2.0;
    01623   /* DIR$ IVDEP */
    01624   for (l = 1; l <= lot; ++l) {
    01625     a[ja] = scale * work[ia];
    01626     a[ja + inc] = -scale * work[ia + 1];
    01627     ia += nx;
    01628     ja += jump;
    01629     /* L40: */
    01630   }
    01631
    01632  L50:
    01633   return;
    01634 } /* fft99b */
    01635
    01636 void fft991(double *a, double *work, double *trigs, long *ifax, long inc,
    01637             long jump, long n, long lot, long isign)
    01638 {
    01639   /* Local variables */
    01640   static long nfax, i, j, k, l, m, ibase, jbase;
    01641   static long ia, ib, la, nh, nx;
    01642   static long igo, ink;
    01643
    01644
    01645 /*     SUBROUTINE ^FFT991^ - MULTIPLE REAL/HALF-COMPLEX PERIODIC */
    01646 /*     FAST FOURIER TRANSFORM */
    01647
    01648 /*     SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO */
    01649 /*     THAT IN MRFFT2 */
    01650
    01651 /*     PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM */
    01652 /*     IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 */
    01653 /*     (1970), 315-337) */
    01654
    01655 /*     A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA */
    01656 /*     WORK IS AN AREA OF SIZE (N+1)*LOT */
    01657 /*     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES */
    01658 /*     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 */
    01659 /*     INC IS THE INCREMENT WITHIN EACH DATA #VECTOR# */
    01660 /*         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) */
    01661 /*     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR */
    01662 /*     N IS THE LENGTH OF THE DATA VECTORS */
    01663 /*     LOT IS THE NUMBER OF DATA VECTORS */
    01664 /*     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT */
    01665 /*           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL */
    01666
    01667 /*     ORDERING OF COEFFICIENTS: */
    01668 /*         A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) */
    01669 /*         WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED */
    01670
    01671 /*     ORDERING OF DATA: */
    01672 /*         X(0),X(1),X(2),...,X(N-1) */
    01673
    01674 /*     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN */
    01675 /*     PARALLEL */
    01676
    01677 /*     *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER */
    01678
    01679 /*     DEFINITION OF TRANSFORMS: */
    01680 /*     ------------------------- */
    01681
    01682 /*     ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) */
    01683 /*         WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) */
    01684
    01685 /*     ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) */
    01686 /*               B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) */
    01687
    01688     /* Parameter adjustments */
    01689   --work;
    01690   --a;
    01691
    01692     /* Function Body */
    01693   nfax = ifax[0];
    01694   nx = n + 1;
    01695   nh = n / 2;
    01696   ink = inc + inc;
    01697   if (isign == 1) {
    01698     goto L30;
    01699   }
    01700
    01701   /*     IF NECESSARY, TRANSFER DATA TO WORK AREA */
    01702   igo = 50;
    01703   if (nfax % 2 == 1) {
    01704     goto L40;
    01705   }
    01706   ibase = 1;
    01707   jbase = 1;
    01708   for (l = 1; l <= lot; ++l) {
    01709     i = ibase;
    01710     j = jbase;
    01711     /* DIR$ IVDEP */
    01712     for (m = 1; m <= n; ++m) {
    01713       work[j] = a[i];
    01714       i += inc;
    01715       ++j;
    01716       /* L10: */
    01717     }
    01718     ibase += jump;
    01719     jbase += nx;
    01720     /* L20: */
    01721   }
    01722
    01723   igo = 60;
    01724   goto L40;
    01725
    01726   /*     PREPROCESSING (ISIGN=+1) */
    01727   /*     ------------------------ */
    01728
    01729  L30:
    01730   fft99a(&a[1], &work[1], trigs, inc, jump, n, lot);
    01731   igo = 60;
    01732
    01733 /*     COMPLEX TRANSFORM */
    01734 /*     ----------------- */
    01735
    01736  L40:
    01737   ia = 1;
    01738   la = 1;
    01739   for (k = 1; k <= nfax; ++k) {
    01740     if (igo == 60) {
    01741       goto L60;
    01742     }
    01743     /* L50: */
    01744     vpassm(&a[ia], &a[ia + inc], &work[1], &work[2], trigs, ink, 2, jump, nx,
    01745            lot, nh, ifax[k], la);
    01746     igo = 60;
    01747     goto L70;
    01748   L60:
    01749     vpassm(&work[1], &work[2], &a[ia], &a[ia + inc], trigs, 2, ink, nx, jump,
    01750            lot, nh, ifax[k], la);
    01751     igo = 50;
    01752   L70:
    01753     la *= ifax[k];
    01754     /* L80: */
    01755   }
    01756
    01757   if (isign == -1) {
    01758     goto L130;
    01759   }
    01760
    01761   /*     IF NECESSARY, TRANSFER DATA FROM WORK AREA */
    01762   if (nfax % 2 == 1) {
    01763     goto L110;
    01764   }
    01765   ibase = 1;
    01766   jbase = 1;
    01767   for (l = 1; l <= lot; ++l) {
    01768     i = ibase;
    01769     j = jbase;
    01770     /* DIR$ IVDEP */
    01771     for (m = 1; m <= n; ++m) {
    01772       a[j] = work[i];
    01773       ++i;
    01774       j += inc;
    01775       /* L90: */
    01776     }
    01777     ibase += nx;
    01778     jbase += jump;
    01779     /* L100: */
    01780   }
    01781
    01782   /*     FILL IN ZEROS AT END */
    01783  L110:
    01784   ib = n * inc + 1;
    01785   /* DIR$ IVDEP */
    01786   for (l = 1; l <= lot; ++l) {
    01787     a[ib] = 0.0;
    01788     a[ib + inc] = 0.0;
    01789     ib += jump;
    01790     /* L120: */
    01791   }
    01792   goto L140;
    01793
    01794 /*     POSTPROCESSING (ISIGN=-1): */
    01795 /*     -------------------------- */
    01796
    01797  L130:
    01798   fft99b(&work[1], &a[1], trigs, inc, jump, n, lot);
    01799
    01800  L140:
    01801   return;
    01802 } /* fft991 */
    01803
    01804 void fftfax(long n, long *ifax, double *trigs)
    01805 {
    01806     /* Local variables */
    01807   static long i;
    01808
    01809   /* Function Body */
    01810   fax(ifax, n);
    01811   i = ifax[0];
    01812   if (ifax[i] > 5 || n <= 4) {
    01813     ifax[0] = -99;
    01814   }
    01815   if (ifax[0] <= 0)
    01816     return;
    01817
    01818   fftrig(trigs, n);
    01819   return;
    01820 } /* fftfax */
    01821
    01822 static void fax(long *ifax, long n)
    01823 {
    01824   /* Local variables */
    01825   static long nfax, item, i, k, l, istop, ii, nn, inc;
    01826
    01827     /* Function Body */
    01828   nn = n / 2;
    01829   if (nn + nn == n) {
    01830     goto L10;
    01831   }
    01832   ifax[0] = -99;
    01833   return;
    01834  L10:
    01835   k = 0;
    01836   /*     TEST FOR FACTORS OF 4 */
    01837  L20:
    01838   if (nn % 4 != 0) {
    01839     goto L30;
    01840   }
    01841   ++k;
    01842   ifax[k] = 4;
    01843   nn /= 4;
    01844   if (nn == 1) {
    01845     goto L80;
    01846   }
    01847   goto L20;
    01848   /*     TEST FOR EXTRA FACTOR OF 2 */
    01849  L30:
    01850   if (nn % 2 != 0) {
    01851     goto L40;
    01852   }
    01853   ++k;
    01854   ifax[k] = 2;
    01855   nn /= 2;
    01856   if (nn == 1) {
    01857     goto L80;
    01858   }
    01859   /*     TEST FOR FACTORS OF 3 */
    01860  L40:
    01861   if (nn % 3 != 0) {
    01862     goto L50;
    01863   }
    01864   ++k;
    01865   ifax[k] = 3;
    01866   nn /= 3;
    01867   if (nn == 1) {
    01868     goto L80;
    01869   }
    01870   goto L40;
    01871   /*     NOW FIND REMAINING FACTORS */
    01872  L50:
    01873   l = 5;
    01874   inc = 2;
    01875   /*     INC ALTERNATELY TAKES ON VALUES 2 AND 4 */
    01876  L60:
    01877   if (nn % l != 0) {
    01878     goto L70;
    01879   }
    01880   ++k;
    01881   ifax[k] = l;
    01882   nn /= l;
    01883   if (nn == 1) {
    01884     goto L80;
    01885   }
    01886   goto L60;
    01887  L70:
    01888   l += inc;
    01889   inc = 6 - inc;
    01890   goto L60;
    01891  L80:
    01892   ifax[0] = k;
    01893   /*     IFAX(0) CONTAINS NUMBER OF FACTORS */
    01894   nfax = ifax[0];
    01895   /*     SORT FACTORS INTO ASCENDING ORDER */
    01896   if (nfax == 1)
    01897     return;
    01898
    01899   for (ii = 2; ii <= nfax; ++ii) {
    01900     istop = nfax + 2 - ii;
    01901     for (i = 2; i <= istop; ++i) {
    01902       if (ifax[i] < ifax[i - 1]) {
    01903         item = ifax[i - 1];
    01904         ifax[i - 1] = ifax[i];
    01905         ifax[i] = item;
    01906       }
    01907     }
    01908   }
    01909   return;
    01910 } /* fax */
    01911
    01912 static void fftrig(double *trigs, long n)
    01913 {
    01914   /* Local variables */
    01915   static long i, n2;
    01916   static double angle;
    01917   static long la, nh;
    01918   static long nn;
    01919   static double del;
    01920
    01921     /* Function Body */
    01922   nn = n / 2;
    01923   del = (PI + PI) / (double) nn;
    01924   n2 = nn + nn;
    01925   for (i = 0; i < n2; i += 2) {
    01926     angle = (double) i * 0.5 * del;
    01927     trigs[i] = cos(angle);
    01928     trigs[i + 1] = sin(angle);
    01929   }
    01930   del *= 0.5;
    01931   nh = (nn + 1) / 2;
    01932   n2 = nh + nh;
    01933   la = nn + nn;
    01934   for (i = 0; i < n2; i += 2) {
    01935     angle = (double) i * 0.5 * del;
    01936     trigs[la + i] = cos(angle);
    01937     trigs[la + i + 1] = sin(angle);
    01938   }
    01939
    01940   return;
    01941 } /* fftrig */
    01942