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 2161 2006-05-19 16:55:03Z paulf $
00006  *
00007  *    Revision history:
00008  *     $Log$
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 

Generated on Tue May 6 09:16:01 2003 for Earthworm Libs by doxygen1.3-rc3