/* vfftpack.f -- translated by f2c (version 20031025).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/
#include <math.h>
#include "f2c.h"

/* Table of constant values */

static doublereal c_b26 = 1.;

/*                     * * * * * * * * * * * * * * */
/*                     *                         * */
/*                     *          VFFTPK         * */
/*                     *                         * */
/*                     * * * * * * * * * * * * * * */

/*                       (Version 2.1  May 1990) */

/*          A vectorized package of Fortran subprograms for the */
/*          fast Fourier transform of multiple real sequences */

/*                                  by */

/*                Roland A. Sweet and Linda L. Lindgren */

/*             Center for Computing and Applied Mathematics */
/*            National Institute of Standards and Technology */
/*                       Boulder, Colorado 80303 */

/*                         Ronald F. Boisvert */

/*             Center for Computing and Applied Mathematics */
/*            National Institute of Standards and Technology */
/*                     Gaithersburg, Maryland 20899 */



/*                     * * * * * * * * * * * * * * */
/*                     *                         * */
/*                     *       Background        * */
/*                     *                         * */
/*                     * * * * * * * * * * * * * * */


/*  NOTATION */
/*  -------- */

/*  To simulate the mathematical symbol 'sigma' that represents */
/*  summation, we define the notation */

/*       SUM(l=is,if)[ y(l) ] = y(is) + y(is+1) + . . . + y(if) . */

/*   DISCRETE FOURIER TRANSFORM */
/*  -------------------------- */

/*  Given a complex periodic sequence x(j), j=0,1,...,n-1, the */
/*  discrete Fourier transform is defined by */

/*      x(j) = sqrt(1/n)*SUM(k=0,n-1)[ f(k)*exp(2j*k*i*pi/n) ] , */

/*  for j = 0, 1, . . . , n-1 , where i = sqrt(-1).  The discrete */
/*  Fourier coefficients, f(k), are given by */

/*      f(k) = sqrt(1/n)*SUM(j=0,n-1)[ x(j)*exp(-2j*k*i*pi/n) ] , */

/*  for k = 0, 1, . . . , n-1 . */
/* C */

/*  REAL PERIODIC TRANSFORM */
/*  ----------------------- */

/*  To facilitate the following discussion we define the integer m by */

/*         m = n/2 when n is even, and */

/*         m = (n+1)/2 when n is odd. */

/*  If the given sequence x(j) is real, the Fourier coefficients are */
/*  conjugate symmetric, i.e. f(n-k) is equal to the complex conjugate */
/*  of f(k).  This relation implies that the sequence x(j) can be */
/*  completely determined from the n real quantities */

/*      real(f(0)) =  sqrt(1/n)*SUM(j=0,n-1)[ x(j) ], */

/*      real(f(k)) =  sqrt(1/n)*SUM(j=0,n-1)[ x(j)*cos(2j*k*pi/n) ] */

/*      imag(f(k)) = -sqrt(1/n)*SUM(j=0,n-1)[ x(j)*sin(2j*k*pi/n) ] */

/*  for k = 1, 2, . . . , m-1, and, when n is even, */

/*      real(f(n/2)) = sqrt(1/n)*SUM(j=0,n-1)[ (-1)**j*x(j) ]. */


/*  TRIGONOMETRIC REPRESENTATION */
/*  ---------------------------- */

/*  When the sequence x(j) is real the discrete Fourier transform may */
/*  be re-written as the trigonometric series */

/*      x(j) = sqrt(2/n)* c(0)/2 + (-1)**j*c(n-1)/2 + SUM(k=1,m-1) */

/*             [c(2k-1)*cos(2j*k*pi/n) + c(2k)*sin(2j*k*pi/n)] */

/*  when n is an even integer (n = 2*m), or as */

/*      x(j) = sqrt(2/n)* c(0)/2 + SUM(k=1,m-1) */

/*             [c(2k-1)*cos(2j*k*pi/n) + c(2k)*sin(2j*k*pi/n)] */

/*  when n is an odd integer (n = 2*m-1), */

/*  for j = 0, 1, . . . , n-1 . */


/*  The coefficients c(k) are defined by the relations */

/*      c(0) = sqrt(2/n)*SUM(j=0,n-1)[ x(j) ] */

/*  for k = 1, 2, . . . , m-1 , */

/*      c(2*k-1) = sqrt(2/n)*SUM(j=0,n-1)[ x(j)*cos(2j*k*pi/n) ] */

/*      c(2*k)   = sqrt(2/n)*SUM(j=0,n-1)[ x(j)*sin(2j*k*pi/n) ] */

/*  and, when n is even */

/*      c(n-1) = sqrt(2/n)*SUM(j=0,n-1)[ (-1)**j*x(j) ] . */


/*  The relationship between the coefficients c(k) and the conjugate */
/*  symmetric Fourier coefficients f(k) is */

/*      c(0) = sqrt(2)*f(0) , */

/*      c(2*k-1) = sqrt(2)*real(f(k)) and c(2*k) = -sqrt(2)*imag(f(k)) */

/*  for k = 1, 2, . . . , m-1, and when n is even */

/*      c(n-1) = sqrt(2)*f(n/2) . */


/*  SINE TRANSFORM OF AN ODD SEQUENCE */
/*  --------------------------------- */

/*  Let x(k) be a real odd sequence, that is, it can be expanded in terms */
/*  of a trigonometric series that contains only sine terms.  The discrete */
/*  odd Fourier transform is given by */

/*     c(k) = SUM(j=1,n-1)[ 2x(j)*sin(j*k*pi/n) ]/sqrt(2n). */

/*  for k = 1, 2, . . ., n-1. */

/*  The inverse transform is identical, that is, */

/*     x(k) = SUM(j=1,n-1)[ 2c(j)*sin(j*k*pi/n) ]/sqrt(2n) */

/*  for k = 1, 2, . . ., n-1. */


/*  COSINE TRANSFORM OF AN EVEN SEQUENCE */
/*  ------------------------------------ */

/*  Let x(k) be a real even sequence, that is, it can be expanded in terms */
/*  of a trigonometric series that contains only cosine terms. The discrete */
/*  even Fourier transform is given by */

/*     c(k) = ( x(0) + (-1)**k*x(n) + */

/*            SUM(j=1,n-1)[ 2x(j)*cos(j*k*pi/n) ] )/sqrt(2n) */

/*  for k = 0, 1, 2, . . ., n. */

/*  The inverse transform is identical, that is, */

/*     x(k) = ( c(0) + (-1)**k*c(n) + */

/*            SUM(j=1,n-1)[ 2c(j)*cos(j*k*pi/n) ] )/sqrt(2n) */

/*  for k = 0, 1, 2, . . ., n. */


/*  QUARTER-WAVE COSINE TRANSFORM */
/*  ----------------------------- */

/*  Let x(k) be a real even quarter-wave sequence, that is, it can be */
/*  expanded in terms of a cosine series with only odd wave numbers. */
/*  The discrete cosine quarter-wave Fourier transform is given by */

/*     c(k) =  ( x(1) + */

/*             SUM(j=2,n)[ 2x(j)*cos((j-1)*(2k-1)*pi/2n) ] )/sqrt(4n) */

/*  for k = 1, 2, . . ., n. The inverse transform is given by */

/*     x(k) =  SUM(j=1,n)[ 4c(j)*cos((2j-1)*(k-1)*pi/2n) ]/sqrt(4n) */

/*  for k = 1, 2, . . ., n. */


/*  QUARTER-WAVE SINE TRANSFORM */
/*  --------------------------- */

/*  Let x(k) be a real odd quarter-wave sequence; that is, it can be */
/*  expanded in terms of a sine series with only odd wave numbers. */
/*  The discrete quarter-wave sine transform is given by */

/*     c(k) =  ( SUM(j=1,n-1)[ 2x(j)*sin(j*(2k-1)*pi/2n) ] */

/*               + (-1)**(k-1)*x(n) )/sqrt(4n) */

/*  for k = 1, 2, . . ., n.  The inverse transform is given by */

/*     x(k) =  ( SUM(j=1,n) [ 4c(j)*sin((2j-1)*k*pi/2n) ] )/sqrt(4n) */

/*  for k = 1, 2, . . ., n. */




/*              * * * * * * * * * * * * * * * * * * * * * */
/*              *                                       * */
/*              *         Description of Package        * */
/*              *                                       * */
/*              * * * * * * * * * * * * * * * * * * * * * */


/*  This package consists of subroutines that compute the transforms */
/*  between sequences x(j) and the Fourier coefficients f(k) */
/*  described above for multiple real sequences.  The real periodic */
/*  transform is a variant of the Stockham autosort algorithm. */
/*  There is no restriction on the length of the sequence.  The */
/*  other specialized transforms are implemented by preprocessing */
/*  and postprocessing. */

/*  In the description below we use the following terms: */

/*    forward transform  = Fourier analysis. Computation of Fourier */
/*                         coefficients f(k) from the sequence x(k). */
/*    backward transform = Fourier synthesis. Computation of the */
/*                         sequence x(k) from its Fourier coefficients */
/*                         f(k). */

/*  The user-callable subroutines are : */

/*     1.  VRFFTI     Initialization routine for VRFFTF and VRFFTB */
/*     2.  VRFFTF     Forward transform of multiple real periodic */
/*                    sequences */
/*     3.  VRFFTB     Backward transform for  multiple  real periodic */
/*                    sequences */
/*     4.  VSINTI     Initialization routine for VSINT */
/*     5.  VSINT      Forward or backward transform for multiple odd */
/*                    sequences */
/*     6.  VCOSTI     Initialization routine for VCOST */
/*     7.  VCOST      Forward or backward transform for multiple even */
/*                    sequences */
/*     8.  VSINQI     Initialization routine for VSINQF and VSINQB */
/*     9.  VSINQF     Forward transform of multiple odd sequences with */
/*                    only odd wave numbers (quarter-wave sine transform) */
/*    10.  VSINQB     Backward transform for multiple odd sequences with */
/*                    only odd wave numbers (quarter-wave sine transform) */
/*    11   VCOSQI     Initialization routine for VCOSQF and VCOSQB */
/*    12   VCOSQF     Forward transform of multiple even sequences with */
/*                    only odd wave numbers (quarter-wave cosine */
/*                    transform) */
/*    13.  VCOSQB     Backward transform of multiple even sequences with */
/*                    only odd wave numbers (quarter-wave cosine */
/*                    transform) */

/* This package is a straightforward vectorization of the scalar package */
/* FFTPACK (Version 3, June 1979) written by Paul N. Swarztrauber of the */
/* National Center for Atmospheric Research, Boulder, Colorado, 80307. */
/* All vector lengths are equal to the number of sequences being */
/* transformed. One slight change has been made: symmetric scaling has */
/* been incorporated into the forward and backward transforms. */


/* Revision history: */

/*   Version 1.0 (Aug 1985) : VRFFTF, VRFFTB, VRFFTI developed by Sweet */
/*                            and Lindgren */
/*   Version 2.0 (Mar 1987) : VSINT, VSINTI, VCOST, VCOSTI, VSINQF, */
/*                            VSINQB, VSINQI, VCOSQF, VCOSQB, VCOSQI */
/*                            added by Boisvert */
/*   Version 2.1 (May 1990) : documentation revised */



/* To get a particular transform, use a request like */
/*    send vrffti vrfftf from vfftpk. */
/* To get the entire library, say */
/*    send all from vfftpk. */

doublereal pimach_(doublereal *dum)
{
    /* System generated locals */
    doublereal ret_val;

/* ***BEGIN PROLOGUE  PIMACH */

/*     This subprogram supplies the value of the constant PI correct to */
/*     machine precision where */

/*     PI=3.1415926535897932384626433832795028841971693993751058209749446 */
/* ***ROUTINES CALLED  (NONE) */
/* ***END PROLOGUE  PIMACH */

/* ***FIRST EXECUTABLE STATEMENT  PIMACH */
    ret_val = 3.14159265358979;
    return ret_val;
} /* pimach_ */

/* Subroutine */ int vcosqb_(integer *m, integer *n, doublereal *x, 
	doublereal *xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k;
    static doublereal x1;
    static integer kc, np2, ns2;
    static doublereal xim1;
    static integer modn;
    static doublereal scale;
    extern /* Subroutine */ int vrfftb_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VCOSQB */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Normalized inverse of VCOSQF. */
/* ***DESCRIPTION */

/*  Subroutine VCOSQB computes the backward fast Fourier cosine transform */
/*  of M quarter wave sequences.  That is, cosine series representations */
/*  with only odd wave numbers.  The transform is defined below at output */
/*  parameter X. */

/*  The array WSAVE which is used by subroutine VCOSQB must be */
/*  initialized by calling subroutine VCOSQI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequences to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  X       an array of size at least X(MDIMX,N) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N. */

/*  XT      a work array of size at least XT(MDIMX,N). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array which must be dimensioned at least 2*N+15 */
/*          in the program that calls VCOSQB.  The WSAVE array must be */
/*          initialized by calling subroutine VCOSQI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged. */

/*  Output Parameters */

/*  X       For I=1,...,N and J=1,...,M */

/*               X(I)= the sum from K=1 to K=N of */

/*                 4*X(K)*COS((2*K-1)*(I-1)*PI/(2*N)) /SQRT(4*N) */

/*  WSAVE   contains initialization calculations which must not */
/*          be destroyed between calls of VCOSQF or VCOSQB. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VCOSQF followed immediately by a call of */
/*           of VCOSQB will return the original sequences X.  Thus, */
/*           VCOSQB is the correctly normalized inverse VCOSQF. */

/*  ----------------------------------------------------------------- */

/*  VCOSQB is a straightforward extension of the subprogram COSQB to */
/*  handle M simultaneous sequences.  COSQB was originally developed */
/*  by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  (NONE) */
/* ***END PROLOGUE  VCOSQB */
/* ***FIRST EXECUTABLE STATEMENT  VCOSQB */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n > 2) {
	goto L300;
    }
    if (*n == 2) {
	goto L200;
    }
    goto L900;

/*  CASE  N = 2 */

L200:
    scale = sqrt(.5) * 2.;
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	x1 = scale * (x[j + x_dim1] + x[j + (x_dim1 << 1)]);
	x[j + (x_dim1 << 1)] = x[j + x_dim1] - x[j + (x_dim1 << 1)];
	x[j + x_dim1] = x1;
/* L210: */
    }
    goto L900;

/*  CASE N .GT. 2 */

/*     ... PREPROCESSING */

L300:
    ns2 = (*n + 1) / 2;
    np2 = *n + 2;
    i__1 = *n;
    for (i__ = 3; i__ <= i__1; i__ += 2) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xim1 = x[j + (i__ - 1) * x_dim1] + x[j + i__ * x_dim1];
	    x[j + i__ * x_dim1] -= x[j + (i__ - 1) * x_dim1];
	    x[j + (i__ - 1) * x_dim1] = xim1;
/* L310: */
	}
    }
    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	x[j + x_dim1] += x[j + x_dim1];
/* L320: */
    }
    modn = *n % 2;
    if (modn == 0) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + *n * x_dim1] += x[j + *n * x_dim1];
/* L330: */
	}
    }

/*     ... REAL, PERIODIC TRANSFORM */

    vrfftb_(m, n, &x[x_offset], &xt[xt_offset], mdimx, &wsave[*n + 1]);

/*     ... POSTPROCESSING */

    i__2 = ns2;
    for (k = 2; k <= i__2; ++k) {
	kc = np2 - k;
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    xt[j + k * xt_dim1] = wsave[k - 1] * x[j + kc * x_dim1] + wsave[
		    kc - 1] * x[j + k * x_dim1];
	    xt[j + kc * xt_dim1] = wsave[k - 1] * x[j + k * x_dim1] - wsave[
		    kc - 1] * x[j + kc * x_dim1];
/* L340: */
	}
    }
    if (modn == 0) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + (ns2 + 1) * x_dim1] = wsave[ns2] * (x[j + (ns2 + 1) * 
		    x_dim1] + x[j + (ns2 + 1) * x_dim1]);
/* L350: */
	}
    }
    i__1 = ns2;
    for (k = 2; k <= i__1; ++k) {
	kc = np2 - k;
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + k * x_dim1] = xt[j + k * xt_dim1] + xt[j + kc * xt_dim1];
	    x[j + kc * x_dim1] = xt[j + k * xt_dim1] - xt[j + kc * xt_dim1];
/* L360: */
	}
    }
    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	x[j + x_dim1] += x[j + x_dim1];
/* L370: */
    }

/*     ... NORMALIZATION */

    scale = .5;
    i__2 = *n;
    for (i__ = 1; i__ <= i__2; ++i__) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + i__ * x_dim1] = scale * x[j + i__ * x_dim1];
/* L380: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vcosqb_ */

/* Subroutine */ int vcosqf_(integer *m, integer *n, doublereal *x, 
	doublereal *xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k, kc, np2, ns2;
    static doublereal xim1;
    static integer modn;
    static doublereal tsqx, sqrt2, scale;
    extern /* Subroutine */ int vrfftf_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VCOSQF */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Forward cosine transform, odd wave numbers, M sequences. */
/* ***DESCRIPTION */

/*  Subroutine VCOSQF computes the forward fast Fourier cosine transform */
/*  of M quarter wave sequences.  That is, cosine series representations */
/*  with only odd wave numbers.  The transform is defined below at output */
/*  parameter X. */

/*  The array WSAVE which is used by subroutine VCOSQF must be */
/*  initialized by calling subroutine VCOSQI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequences to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  X       an array of size at least X(MDIMX,N) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N. */

/*  XT      a work array of size at least XT(MDIMX,N). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array which must be dimensioned at least 2*N+15 */
/*          in the program that calls VCOSQF.  The WSAVE array must be */
/*          initialized by calling subroutine VCOSQI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged. */

/*  Output Parameters */

/*  X       For I=1,...,N and J=1,...,M */

/*               X(I) = ( X(1) + the sum from K=2 to K=N of */

/*                  2*X(K)*COS((2*I-1)*(K-1)*PI/(2*N)) )/SQRT(4*N) */

/*  WSAVE   contains initialization calculations which must not */
/*          be destroyed between calls of VCOSQF or VCOSQB. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VCOSQF followed immediately by a call of */
/*           of VCOSQB will return the original sequences X.  Thus, */
/*           VCOSQB is the correctly normalized inverse VCOSQF. */

/*  ----------------------------------------------------------------- */

/*  VCOSQF is a straightforward extension of the subprogram COSQF to */
/*  handle M simultaneous sequences.  COSQF was originally developed */
/*  by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  (NONE) */
/* ***END PROLOGUE  VCOSQF */
/* ***FIRST EXECUTABLE STATEMENT  VCOSQF */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n > 2) {
	goto L300;
    }
    if (*n < 2) {
	goto L900;
    }

/*  CASE  N = 2 */

    sqrt2 = sqrt(2.);
    scale = .5 / sqrt2;
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	tsqx = sqrt2 * x[j + (x_dim1 << 1)];
	x[j + (x_dim1 << 1)] = scale * (x[j + x_dim1] - tsqx);
	x[j + x_dim1] = scale * (x[j + x_dim1] + tsqx);
/* L210: */
    }
    goto L900;

/*  CASE N .GT. 2 */

L300:

/*     ... PREPROCESSING */

    ns2 = (*n + 1) / 2;
    np2 = *n + 2;
    i__1 = ns2;
    for (k = 2; k <= i__1; ++k) {
	kc = np2 - k;
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xt[j + k * xt_dim1] = x[j + k * x_dim1] + x[j + kc * x_dim1];
	    xt[j + kc * xt_dim1] = x[j + k * x_dim1] - x[j + kc * x_dim1];
/* L310: */
	}
    }
    modn = *n % 2;
    if (modn == 0) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xt[j + (ns2 + 1) * xt_dim1] = x[j + (ns2 + 1) * x_dim1] + x[j + (
		    ns2 + 1) * x_dim1];
/* L320: */
	}
    }
    i__2 = ns2;
    for (k = 2; k <= i__2; ++k) {
	kc = np2 - k;
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + k * x_dim1] = wsave[k - 1] * xt[j + kc * xt_dim1] + wsave[
		    kc - 1] * xt[j + k * xt_dim1];
	    x[j + kc * x_dim1] = wsave[k - 1] * xt[j + k * xt_dim1] - wsave[
		    kc - 1] * xt[j + kc * xt_dim1];
/* L330: */
	}
    }
    if (modn == 0) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + (ns2 + 1) * x_dim1] = wsave[ns2] * xt[j + (ns2 + 1) * 
		    xt_dim1];
/* L340: */
	}
    }

/*     ... REAL, PERIODIC TRANSFORM */

    vrfftf_(m, n, &x[x_offset], &xt[xt_offset], mdimx, &wsave[*n + 1]);

/*     ... POSTPROCESSING */

    i__1 = *n;
    for (i__ = 3; i__ <= i__1; i__ += 2) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xim1 = x[j + (i__ - 1) * x_dim1] - x[j + i__ * x_dim1];
	    x[j + i__ * x_dim1] = x[j + (i__ - 1) * x_dim1] + x[j + i__ * 
		    x_dim1];
	    x[j + (i__ - 1) * x_dim1] = xim1;
/* L350: */
	}
    }

/*     ... NORMALIZATION */

    scale = .5;
    i__2 = *n;
    for (i__ = 1; i__ <= i__2; ++i__) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + i__ * x_dim1] = scale * x[j + i__ * x_dim1];
/* L360: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vcosqf_ */

/* Subroutine */ int vcosqi_(integer *n, doublereal *wsave)
{
    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    double cos(doublereal);

    /* Local variables */
    static integer k;
    static doublereal fk, dt, pih;
    extern doublereal pimach_(doublereal *);
    extern /* Subroutine */ int vrffti_(integer *, doublereal *);

/* ***BEGIN PROLOGUE  VCOSQI */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Initialize for VCOSQF and VCOSQB. */
/* ***DESCRIPTION */

/*  Subroutine VCOSQI initializes the array WSAVE which is used in */
/*  both VCOSQF and VCOSQB.  The prime factorization of N together with */
/*  a tabulation of the trigonometric functions are computed and */
/*  stored in WSAVE. */

/*  Input Parameter */

/*  N       the length of the array to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  Output Parameter */

/*  WSAVE   a work array which must be dimensioned at least 2*N+15. */
/*          The same work array can be used for both VCOSQF and VCOSQB */
/*          as long as N remains unchanged.  Different WSAVE arrays */
/*          are required for different values of N.  The contents of */
/*          WSAVE must not be changed between calls of VCOSQF or VCOSQB. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFFTI */
/* ***END PROLOGUE  VCOSQI */
/* ***FIRST EXECUTABLE STATEMENT  VCOSQI */
    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    pih = pimach_(&c_b26) * .5;
    dt = pih / (doublereal) (*n);
    fk = 0.;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	fk += 1.;
	wsave[k] = cos(fk * dt);
/* L101: */
    }
    vrffti_(n, &wsave[*n + 1]);
    return 0;
} /* vcosqi_ */

/* Subroutine */ int vcost_(integer *m, integer *n, doublereal *x, doublereal 
	*xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k;
    static doublereal t1, t2;
    static integer kc;
    static doublereal xi;
    static integer nm1, np1;
    static doublereal x1h;
    static integer ns2;
    static doublereal tx2, x1p3;
    static integer modn;
    static doublereal scale, factor;
    extern /* Subroutine */ int vrfftf_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VCOST */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE */
/*             SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Cosine transform of one or more real, even sequences. */
/* ***DESCRIPTION */

/*  Subroutine VCOST computes the discrete Fourier cosine transform */
/*  of M even sequences X(J,I), J=1,...,M.  The transform is defined */
/*  below at output parameter X. */

/*  The array WSAVE which is used by subroutine VCOST must be */
/*  initialized by calling subroutine VCOSTI(N,WSAVE). */

/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequence to be transformed.  N must be */
/*          greater than 1.  The method is most efficient when N-1 is */
/*          is a product of small primes. */

/*  X       an array of size at least X(MDIMX,N) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N. */

/*  XT      a work array of size at least XT(MDIMX,N-1). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array which must be dimensioned at least 3*N+15 */
/*          in the program that calls VCOST.  The WSAVE array must be */
/*          initialized by calling subroutine VCOSTI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged.  Thus subsequent */
/*          transforms can be obtained faster than the first. */

/*  Output Parameters */

/*  X       For I=1,...,N and J=1,...,M */

/*             X(J,I) = ( X(J,1)+(-1)**(I-1)*X(J,N) */

/*               + the sum from K=2 to K=N-1 */

/*                 2*X(J,K)*COS((K-1)*(I-1)*PI/(N-1)) )/SQRT(2*(N-1)) */

/*  WSAVE   contains initialization calculations which must not be */
/*          destroyed between calls of VCOST. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VCOST followed immediately by another call */
/*           of VCOST will return the original sequences X.  Thus, */
/*           VCOST is the correctly normalized inverse of itself. */

/*  ----------------------------------------------------------------- */

/*  VCOST is a straightforward extension of the subprogram COST to */
/*  handle M simultaneous sequences.  The scaling of the sequences */
/*  computed by VCOST is different than that of COST.  COST was */
/*  originally developed by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFFTF */
/* ***END PROLOGUE  VCOST */
/* ***FIRST EXECUTABLE STATEMENT  VCOST */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n <= 1) {
	goto L900;
    }
    if (*n > 3) {
	goto L400;
    }
    if (*n == 3) {
	goto L300;
    }

/*  CASE  N = 2 */

    scale = sqrt(.5);
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	x1h = scale * (x[j + x_dim1] + x[j + (x_dim1 << 1)]);
	x[j + (x_dim1 << 1)] = scale * (x[j + x_dim1] - x[j + (x_dim1 << 1)]);
	x[j + x_dim1] = x1h;
/* L210: */
    }
    goto L900;

/*  CASE  N = 3 */

L300:
    scale = .5;
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	x1p3 = x[j + x_dim1] + x[j + x_dim1 * 3];
	tx2 = x[j + (x_dim1 << 1)] + x[j + (x_dim1 << 1)];
	x[j + (x_dim1 << 1)] = scale * (x[j + x_dim1] - x[j + x_dim1 * 3]);
	x[j + x_dim1] = scale * (x1p3 + tx2);
	x[j + x_dim1 * 3] = scale * (x1p3 - tx2);
/* L310: */
    }
    goto L900;

/*  CASE  N .GT. 3 */

/*     ... PREPROCESSING */

L400:
    nm1 = *n - 1;
    np1 = *n + 1;
    ns2 = *n / 2;
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	xt[j + xt_dim1] = x[j + x_dim1] - x[j + *n * x_dim1];
	x[j + x_dim1] += x[j + *n * x_dim1];
/* L410: */
    }
    i__1 = ns2;
    for (k = 2; k <= i__1; ++k) {
	kc = np1 - k;
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    t1 = x[j + k * x_dim1] + x[j + kc * x_dim1];
	    t2 = x[j + k * x_dim1] - x[j + kc * x_dim1];
	    xt[j + xt_dim1] += wsave[kc] * t2;
	    t2 = wsave[k] * t2;
	    x[j + k * x_dim1] = t1 - t2;
	    x[j + kc * x_dim1] = t1 + t2;
/* L420: */
	}
    }
    modn = *n % 2;
    if (modn != 0) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + (ns2 + 1) * x_dim1] += x[j + (ns2 + 1) * x_dim1];
/* L430: */
	}
    }
    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	x[j + *n * x_dim1] = xt[j + xt_dim1];
/* L435: */
    }

/*     ... REAL PERIODIC TRANSFORM */

    vrfftf_(m, &nm1, &x[x_offset], &xt[xt_offset], mdimx, &wsave[np1]);

/*     ... POSTPROCESSING */

    factor = 1. / sqrt((doublereal) nm1);
    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	xt[j + xt_dim1] = x[j + (x_dim1 << 1)];
	x[j + (x_dim1 << 1)] = factor * x[j + *n * x_dim1];
/* L440: */
    }
    i__2 = *n;
    for (i__ = 4; i__ <= i__2; i__ += 2) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    xi = x[j + i__ * x_dim1];
	    x[j + i__ * x_dim1] = x[j + (i__ - 2) * x_dim1] - x[j + (i__ - 1) 
		    * x_dim1];
	    x[j + (i__ - 1) * x_dim1] = xt[j + xt_dim1];
	    xt[j + xt_dim1] = xi;
/* L450: */
	}
    }
    if (modn != 0) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + *n * x_dim1] = xt[j + xt_dim1];
/* L460: */
	}
    }

/*     ... NORMALIZATION */

    scale = sqrt(.5);
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + i__ * x_dim1] = scale * x[j + i__ * x_dim1];
/* L490: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vcost_ */

/* Subroutine */ int vcosti_(integer *n, doublereal *wsave)
{
    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    double sin(doublereal), cos(doublereal);

    /* Local variables */
    static integer k, kc;
    static doublereal fk, dt, pi;
    static integer nm1, np1, ns2;
    extern doublereal pimach_(doublereal *);
    extern /* Subroutine */ int vrffti_(integer *, doublereal *);

/* ***BEGIN PROLOGUE  VCOSTI */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE */
/*             SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Initialize for VCOST. */
/* ***DESCRIPTION */

/*  Subroutine VCOSTI initializes the array WSAVE which is used in */
/*  subroutine VCOST.  The prime factorization of N together with */
/*  a tabulation of the trigonometric functions are computed and */
/*  stored in WSAVE. */

/*  Input Parameter */

/*  N       the length of the sequence to be transformed.  The method */
/*          is most efficient when N-1 is a product of small primes. */

/*  Output Parameter */

/*  WSAVE   a work array which must be dimensioned at least 3*N+15. */
/*          Different WSAVE arrays are required for different values */
/*          of N.  The contents of WSAVE must not be changed between */
/*          calls of VCOST. */

/*  ----------------------------------------------------------------- */

/*  VCOSTI is a straightforward extension of the subprogram COSTI to */
/*  handle M simultaneous sequences.  COSTI was originally developed */
/*  by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFFTI */
/* ***END PROLOGUE  VCOSTI */
/* ***FIRST EXECUTABLE STATEMENT  VCOSTI */
    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    pi = pimach_(&c_b26);
    if (*n <= 3) {
	return 0;
    }
    nm1 = *n - 1;
    np1 = *n + 1;
    ns2 = *n / 2;
    dt = pi / (doublereal) nm1;
    fk = 0.;
    i__1 = ns2;
    for (k = 2; k <= i__1; ++k) {
	fk += 1.;
	wsave[k] = sin(fk * dt) * 2.;
/* L101: */
    }
    fk = 0.;
    i__1 = ns2;
    for (k = 2; k <= i__1; ++k) {
	kc = np1 - k;
	fk += 1.;
	wsave[kc] = cos(fk * dt) * 2.;
/* L102: */
    }
    vrffti_(&nm1, &wsave[*n + 1]);
    return 0;
} /* vcosti_ */

/* Subroutine */ int vradb2_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_dim3, ch_offset,
	     i__1, i__2, i__3;

    /* Local variables */
    static integer i__, k, m, ic, idp2;


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * 3);
    cc -= cc_offset;

    /* Function Body */
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k + ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + (((k << 
		    1) + 1) * cc_dim2 + 1) * cc_dim1] + cc[m + (*ido + ((k << 
		    1) + 2) * cc_dim2) * cc_dim1];
	    ch[m + ((k + (ch_dim3 << 1)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    (k << 1) + 1) * cc_dim2 + 1) * cc_dim1] - cc[m + (*ido + (
		    (k << 1) + 2) * cc_dim2) * cc_dim1];
/* L1001: */
	}
/* L101: */
    }
    if ((i__1 = *ido - 2) < 0) {
	goto L107;
    } else if (i__1 == 0) {
	goto L105;
    } else {
	goto L102;
    }
L102:
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + ((k << 1) + 1) * cc_dim2) * cc_dim1] + 
			cc[m + (ic - 1 + ((k << 1) + 2) * cc_dim2) * cc_dim1];
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + ((k << 1) + 1) * cc_dim2) * cc_dim1] - cc[m + (
			ic + ((k << 1) + 2) * cc_dim2) * cc_dim1];
		ch[m + (i__ - 1 + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] =
			 wa1[i__ - 2] * (cc[m + (i__ - 1 + ((k << 1) + 1) * 
			cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((k << 1) + 2) 
			* cc_dim2) * cc_dim1]) - wa1[i__ - 1] * (cc[m + (i__ 
			+ ((k << 1) + 1) * cc_dim2) * cc_dim1] + cc[m + (ic + 
			((k << 1) + 2) * cc_dim2) * cc_dim1]);
		ch[m + (i__ + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = 
			wa1[i__ - 2] * (cc[m + (i__ + ((k << 1) + 1) * 
			cc_dim2) * cc_dim1] + cc[m + (ic + ((k << 1) + 2) * 
			cc_dim2) * cc_dim1]) + wa1[i__ - 1] * (cc[m + (i__ - 
			1 + ((k << 1) + 1) * cc_dim2) * cc_dim1] - cc[m + (ic 
			- 1 + ((k << 1) + 2) * cc_dim2) * cc_dim1]);
/* L1002: */
	    }
/* L103: */
	}
/* L104: */
    }
    if (*ido % 2 == 1) {
	return 0;
    }
L105:
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (*ido + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (*ido 
		    + ((k << 1) + 1) * cc_dim2) * cc_dim1] + cc[m + (*ido + ((
		    k << 1) + 1) * cc_dim2) * cc_dim1];
	    ch[m + (*ido + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = -(cc[
		    m + (((k << 1) + 2) * cc_dim2 + 1) * cc_dim1] + cc[m + (((
		    k << 1) + 2) * cc_dim2 + 1) * cc_dim1]);
/* L1003: */
	}
/* L106: */
    }
L107:
    return 0;
} /* vradb2_ */

/* Subroutine */ int vradb3_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_dim3, ch_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, k, m, ic;
    static doublereal arg;
    static integer idp2;
    static doublereal taui, taur;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_offset = 1 + cc_dim1 * (1 + (cc_dim2 << 2));
    cc -= cc_offset;

    /* Function Body */
    arg = pimach_(&c_b26) * 2. / 3.;
    taur = cos(arg);
    taui = sin(arg);
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k + ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k * 3 
		    + 1) * cc_dim2 + 1) * cc_dim1] + cc[m + (*ido + (k * 3 + 
		    2) * cc_dim2) * cc_dim1] * 2.;
	    ch[m + ((k + (ch_dim3 << 1)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    k * 3 + 1) * cc_dim2 + 1) * cc_dim1] + taur * 2. * cc[m + 
		    (*ido + (k * 3 + 2) * cc_dim2) * cc_dim1] - taui * 2. * 
		    cc[m + ((k * 3 + 3) * cc_dim2 + 1) * cc_dim1];
	    ch[m + ((k + ch_dim3 * 3) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k *
		     3 + 1) * cc_dim2 + 1) * cc_dim1] + taur * 2. * cc[m + (*
		    ido + (k * 3 + 2) * cc_dim2) * cc_dim1] + taui * 2. * cc[
		    m + ((k * 3 + 3) * cc_dim2 + 1) * cc_dim1];
/* L1001: */
	}
/* L101: */
    }
    if (*ido == 1) {
	return 0;
    }
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k * 3 + 1) * cc_dim2) * cc_dim1] + (cc[
			m + (i__ - 1 + (k * 3 + 3) * cc_dim2) * cc_dim1] + cc[
			m + (ic - 1 + (k * 3 + 2) * cc_dim2) * cc_dim1]);
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + (k * 3 + 1) * cc_dim2) * cc_dim1] + (cc[m + (
			i__ + (k * 3 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic 
			+ (k * 3 + 2) * cc_dim2) * cc_dim1]);
		ch[m + (i__ - 1 + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] =
			 wa1[i__ - 2] * (cc[m + (i__ - 1 + (k * 3 + 1) * 
			cc_dim2) * cc_dim1] + taur * (cc[m + (i__ - 1 + (k * 
			3 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			3 + 2) * cc_dim2) * cc_dim1]) - taui * (cc[m + (i__ + 
			(k * 3 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k * 
			3 + 2) * cc_dim2) * cc_dim1])) - wa1[i__ - 1] * (cc[m 
			+ (i__ + (k * 3 + 1) * cc_dim2) * cc_dim1] + taur * (
			cc[m + (i__ + (k * 3 + 3) * cc_dim2) * cc_dim1] - cc[
			m + (ic + (k * 3 + 2) * cc_dim2) * cc_dim1]) + taui * 
			(cc[m + (i__ - 1 + (k * 3 + 3) * cc_dim2) * cc_dim1] 
			- cc[m + (ic - 1 + (k * 3 + 2) * cc_dim2) * cc_dim1]))
			;
		ch[m + (i__ + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = 
			wa1[i__ - 2] * (cc[m + (i__ + (k * 3 + 1) * cc_dim2) *
			 cc_dim1] + taur * (cc[m + (i__ + (k * 3 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 3 + 2) * 
			cc_dim2) * cc_dim1]) + taui * (cc[m + (i__ - 1 + (k * 
			3 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			3 + 2) * cc_dim2) * cc_dim1])) + wa1[i__ - 1] * (cc[m 
			+ (i__ - 1 + (k * 3 + 1) * cc_dim2) * cc_dim1] + taur 
			* (cc[m + (i__ - 1 + (k * 3 + 3) * cc_dim2) * cc_dim1]
			 + cc[m + (ic - 1 + (k * 3 + 2) * cc_dim2) * cc_dim1])
			 - taui * (cc[m + (i__ + (k * 3 + 3) * cc_dim2) * 
			cc_dim1] + cc[m + (ic + (k * 3 + 2) * cc_dim2) * 
			cc_dim1]));
		ch[m + (i__ - 1 + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = 
			wa2[i__ - 2] * (cc[m + (i__ - 1 + (k * 3 + 1) * 
			cc_dim2) * cc_dim1] + taur * (cc[m + (i__ - 1 + (k * 
			3 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			3 + 2) * cc_dim2) * cc_dim1]) + taui * (cc[m + (i__ + 
			(k * 3 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k * 
			3 + 2) * cc_dim2) * cc_dim1])) - wa2[i__ - 1] * (cc[m 
			+ (i__ + (k * 3 + 1) * cc_dim2) * cc_dim1] + taur * (
			cc[m + (i__ + (k * 3 + 3) * cc_dim2) * cc_dim1] - cc[
			m + (ic + (k * 3 + 2) * cc_dim2) * cc_dim1]) - taui * 
			(cc[m + (i__ - 1 + (k * 3 + 3) * cc_dim2) * cc_dim1] 
			- cc[m + (ic - 1 + (k * 3 + 2) * cc_dim2) * cc_dim1]))
			;
		ch[m + (i__ + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = wa2[
			i__ - 2] * (cc[m + (i__ + (k * 3 + 1) * cc_dim2) * 
			cc_dim1] + taur * (cc[m + (i__ + (k * 3 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 3 + 2) * 
			cc_dim2) * cc_dim1]) - taui * (cc[m + (i__ - 1 + (k * 
			3 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			3 + 2) * cc_dim2) * cc_dim1])) + wa2[i__ - 1] * (cc[m 
			+ (i__ - 1 + (k * 3 + 1) * cc_dim2) * cc_dim1] + taur 
			* (cc[m + (i__ - 1 + (k * 3 + 3) * cc_dim2) * cc_dim1]
			 + cc[m + (ic - 1 + (k * 3 + 2) * cc_dim2) * cc_dim1])
			 + taui * (cc[m + (i__ + (k * 3 + 3) * cc_dim2) * 
			cc_dim1] + cc[m + (ic + (k * 3 + 2) * cc_dim2) * 
			cc_dim1]));
/* L1002: */
	    }
/* L102: */
	}
/* L103: */
    }
    return 0;
} /* vradb3_ */

/* Subroutine */ int vradb4_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2, doublereal *wa3)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_dim3, ch_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, k, m, ic, idp2;
    static doublereal sqrt2;


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa3;
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * 5);
    cc -= cc_offset;

    /* Function Body */
    sqrt2 = sqrt(2.);
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k + ch_dim3 * 3) * ch_dim2 + 1) * ch_dim1] = cc[m + (((k 
		    << 2) + 1) * cc_dim2 + 1) * cc_dim1] + cc[m + (*ido + ((k 
		    << 2) + 4) * cc_dim2) * cc_dim1] - (cc[m + (*ido + ((k << 
		    2) + 2) * cc_dim2) * cc_dim1] + cc[m + (*ido + ((k << 2) 
		    + 2) * cc_dim2) * cc_dim1]);
	    ch[m + ((k + ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + (((k << 
		    2) + 1) * cc_dim2 + 1) * cc_dim1] + cc[m + (*ido + ((k << 
		    2) + 4) * cc_dim2) * cc_dim1] + (cc[m + (*ido + ((k << 2) 
		    + 2) * cc_dim2) * cc_dim1] + cc[m + (*ido + ((k << 2) + 2)
		     * cc_dim2) * cc_dim1]);
	    ch[m + ((k + (ch_dim3 << 2)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    (k << 2) + 1) * cc_dim2 + 1) * cc_dim1] - cc[m + (*ido + (
		    (k << 2) + 4) * cc_dim2) * cc_dim1] + (cc[m + (((k << 2) 
		    + 3) * cc_dim2 + 1) * cc_dim1] + cc[m + (((k << 2) + 3) * 
		    cc_dim2 + 1) * cc_dim1]);
	    ch[m + ((k + (ch_dim3 << 1)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    (k << 2) + 1) * cc_dim2 + 1) * cc_dim1] - cc[m + (*ido + (
		    (k << 2) + 4) * cc_dim2) * cc_dim1] - (cc[m + (((k << 2) 
		    + 3) * cc_dim2 + 1) * cc_dim1] + cc[m + (((k << 2) + 3) * 
		    cc_dim2 + 1) * cc_dim1]);
/* L1001: */
	}
/* L101: */
    }
    if ((i__1 = *ido - 2) < 0) {
	goto L107;
    } else if (i__1 == 0) {
	goto L105;
    } else {
	goto L102;
    }
L102:
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + ((k << 2) + 1) * cc_dim2) * cc_dim1] + 
			cc[m + (ic - 1 + ((k << 2) + 4) * cc_dim2) * cc_dim1] 
			+ (cc[m + (i__ - 1 + ((k << 2) + 3) * cc_dim2) * 
			cc_dim1] + cc[m + (ic - 1 + ((k << 2) + 2) * cc_dim2) 
			* cc_dim1]);
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + ((k << 2) + 1) * cc_dim2) * cc_dim1] - cc[m + (
			ic + ((k << 2) + 4) * cc_dim2) * cc_dim1] + (cc[m + (
			i__ + ((k << 2) + 3) * cc_dim2) * cc_dim1] - cc[m + (
			ic + ((k << 2) + 2) * cc_dim2) * cc_dim1]);
		ch[m + (i__ - 1 + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] =
			 wa1[i__ - 2] * (cc[m + (i__ - 1 + ((k << 2) + 1) * 
			cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((k << 2) + 4) 
			* cc_dim2) * cc_dim1] - (cc[m + (i__ + ((k << 2) + 3) 
			* cc_dim2) * cc_dim1] + cc[m + (ic + ((k << 2) + 2) * 
			cc_dim2) * cc_dim1])) - wa1[i__ - 1] * (cc[m + (i__ + 
			((k << 2) + 1) * cc_dim2) * cc_dim1] + cc[m + (ic + ((
			k << 2) + 4) * cc_dim2) * cc_dim1] + (cc[m + (i__ - 1 
			+ ((k << 2) + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + ((k << 2) + 2) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = 
			wa1[i__ - 2] * (cc[m + (i__ + ((k << 2) + 1) * 
			cc_dim2) * cc_dim1] + cc[m + (ic + ((k << 2) + 4) * 
			cc_dim2) * cc_dim1] + (cc[m + (i__ - 1 + ((k << 2) + 
			3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((k << 2) 
			+ 2) * cc_dim2) * cc_dim1])) + wa1[i__ - 1] * (cc[m + 
			(i__ - 1 + ((k << 2) + 1) * cc_dim2) * cc_dim1] - cc[
			m + (ic - 1 + ((k << 2) + 4) * cc_dim2) * cc_dim1] - (
			cc[m + (i__ + ((k << 2) + 3) * cc_dim2) * cc_dim1] + 
			cc[m + (ic + ((k << 2) + 2) * cc_dim2) * cc_dim1]));
		ch[m + (i__ - 1 + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = 
			wa2[i__ - 2] * (cc[m + (i__ - 1 + ((k << 2) + 1) * 
			cc_dim2) * cc_dim1] + cc[m + (ic - 1 + ((k << 2) + 4) 
			* cc_dim2) * cc_dim1] - (cc[m + (i__ - 1 + ((k << 2) 
			+ 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + ((k << 
			2) + 2) * cc_dim2) * cc_dim1])) - wa2[i__ - 1] * (cc[
			m + (i__ + ((k << 2) + 1) * cc_dim2) * cc_dim1] - cc[
			m + (ic + ((k << 2) + 4) * cc_dim2) * cc_dim1] - (cc[
			m + (i__ + ((k << 2) + 3) * cc_dim2) * cc_dim1] - cc[
			m + (ic + ((k << 2) + 2) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = wa2[
			i__ - 2] * (cc[m + (i__ + ((k << 2) + 1) * cc_dim2) * 
			cc_dim1] - cc[m + (ic + ((k << 2) + 4) * cc_dim2) * 
			cc_dim1] - (cc[m + (i__ + ((k << 2) + 3) * cc_dim2) * 
			cc_dim1] - cc[m + (ic + ((k << 2) + 2) * cc_dim2) * 
			cc_dim1])) + wa2[i__ - 1] * (cc[m + (i__ - 1 + ((k << 
			2) + 1) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + ((k 
			<< 2) + 4) * cc_dim2) * cc_dim1] - (cc[m + (i__ - 1 + 
			((k << 2) + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 
			+ ((k << 2) + 2) * cc_dim2) * cc_dim1]));
		ch[m + (i__ - 1 + (k + (ch_dim3 << 2)) * ch_dim2) * ch_dim1] =
			 wa3[i__ - 2] * (cc[m + (i__ - 1 + ((k << 2) + 1) * 
			cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((k << 2) + 4) 
			* cc_dim2) * cc_dim1] + (cc[m + (i__ + ((k << 2) + 3) 
			* cc_dim2) * cc_dim1] + cc[m + (ic + ((k << 2) + 2) * 
			cc_dim2) * cc_dim1])) - wa3[i__ - 1] * (cc[m + (i__ + 
			((k << 2) + 1) * cc_dim2) * cc_dim1] + cc[m + (ic + ((
			k << 2) + 4) * cc_dim2) * cc_dim1] - (cc[m + (i__ - 1 
			+ ((k << 2) + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + ((k << 2) + 2) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k + (ch_dim3 << 2)) * ch_dim2) * ch_dim1] = 
			wa3[i__ - 2] * (cc[m + (i__ + ((k << 2) + 1) * 
			cc_dim2) * cc_dim1] + cc[m + (ic + ((k << 2) + 4) * 
			cc_dim2) * cc_dim1] - (cc[m + (i__ - 1 + ((k << 2) + 
			3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((k << 2) 
			+ 2) * cc_dim2) * cc_dim1])) + wa3[i__ - 1] * (cc[m + 
			(i__ - 1 + ((k << 2) + 1) * cc_dim2) * cc_dim1] - cc[
			m + (ic - 1 + ((k << 2) + 4) * cc_dim2) * cc_dim1] + (
			cc[m + (i__ + ((k << 2) + 3) * cc_dim2) * cc_dim1] + 
			cc[m + (ic + ((k << 2) + 2) * cc_dim2) * cc_dim1]));
/* L1002: */
	    }
/* L103: */
	}
/* L104: */
    }
    if (*ido % 2 == 1) {
	return 0;
    }
L105:
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (*ido + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (*ido 
		    + ((k << 2) + 1) * cc_dim2) * cc_dim1] + cc[m + (*ido + ((
		    k << 2) + 3) * cc_dim2) * cc_dim1] + (cc[m + (*ido + ((k 
		    << 2) + 1) * cc_dim2) * cc_dim1] + cc[m + (*ido + ((k << 
		    2) + 3) * cc_dim2) * cc_dim1]);
	    ch[m + (*ido + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = sqrt2 
		    * (cc[m + (*ido + ((k << 2) + 1) * cc_dim2) * cc_dim1] - 
		    cc[m + (*ido + ((k << 2) + 3) * cc_dim2) * cc_dim1] - (cc[
		    m + (((k << 2) + 2) * cc_dim2 + 1) * cc_dim1] + cc[m + (((
		    k << 2) + 4) * cc_dim2 + 1) * cc_dim1]));
	    ch[m + (*ido + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = cc[m + ((
		    (k << 2) + 4) * cc_dim2 + 1) * cc_dim1] - cc[m + (((k << 
		    2) + 2) * cc_dim2 + 1) * cc_dim1] + (cc[m + (((k << 2) + 
		    4) * cc_dim2 + 1) * cc_dim1] - cc[m + (((k << 2) + 2) * 
		    cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + (k + (ch_dim3 << 2)) * ch_dim2) * ch_dim1] = 
		    -sqrt2 * (cc[m + (*ido + ((k << 2) + 1) * cc_dim2) * 
		    cc_dim1] - cc[m + (*ido + ((k << 2) + 3) * cc_dim2) * 
		    cc_dim1] + (cc[m + (((k << 2) + 2) * cc_dim2 + 1) * 
		    cc_dim1] + cc[m + (((k << 2) + 4) * cc_dim2 + 1) * 
		    cc_dim1]));
/* L1003: */
	}
/* L106: */
    }
L107:
    return 0;
} /* vradb4_ */

/* Subroutine */ int vradb5_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2, doublereal *wa3, doublereal *wa4)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_dim3, ch_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, k, m, ic;
    static doublereal arg, ti11, ti12, tr11, tr12;
    static integer idp2;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa4;
    --wa3;
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * 6);
    cc -= cc_offset;

    /* Function Body */
    arg = pimach_(&c_b26) * 2. / 5.;
    tr11 = cos(arg);
    ti11 = sin(arg);
    tr12 = cos(arg * 2.);
    ti12 = sin(arg * 2.);
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k + ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k * 5 
		    + 1) * cc_dim2 + 1) * cc_dim1] + cc[m + (*ido + (k * 5 + 
		    2) * cc_dim2) * cc_dim1] * 2. + cc[m + (*ido + (k * 5 + 4)
		     * cc_dim2) * cc_dim1] * 2.;
	    ch[m + ((k + (ch_dim3 << 1)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    k * 5 + 1) * cc_dim2 + 1) * cc_dim1] + tr11 * 2. * cc[m + 
		    (*ido + (k * 5 + 2) * cc_dim2) * cc_dim1] + tr12 * 2. * 
		    cc[m + (*ido + (k * 5 + 4) * cc_dim2) * cc_dim1] - (ti11 *
		     2. * cc[m + ((k * 5 + 3) * cc_dim2 + 1) * cc_dim1] + 
		    ti12 * 2. * cc[m + ((k * 5 + 5) * cc_dim2 + 1) * cc_dim1])
		    ;
	    ch[m + ((k + ch_dim3 * 3) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k *
		     5 + 1) * cc_dim2 + 1) * cc_dim1] + tr12 * 2. * cc[m + (*
		    ido + (k * 5 + 2) * cc_dim2) * cc_dim1] + tr11 * 2. * cc[
		    m + (*ido + (k * 5 + 4) * cc_dim2) * cc_dim1] - (ti12 * 
		    2. * cc[m + ((k * 5 + 3) * cc_dim2 + 1) * cc_dim1] - ti11 
		    * 2. * cc[m + ((k * 5 + 5) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + ((k + (ch_dim3 << 2)) * ch_dim2 + 1) * ch_dim1] = cc[m + ((
		    k * 5 + 1) * cc_dim2 + 1) * cc_dim1] + tr12 * 2. * cc[m + 
		    (*ido + (k * 5 + 2) * cc_dim2) * cc_dim1] + tr11 * 2. * 
		    cc[m + (*ido + (k * 5 + 4) * cc_dim2) * cc_dim1] + (ti12 *
		     2. * cc[m + ((k * 5 + 3) * cc_dim2 + 1) * cc_dim1] - 
		    ti11 * 2. * cc[m + ((k * 5 + 5) * cc_dim2 + 1) * cc_dim1])
		    ;
	    ch[m + ((k + ch_dim3 * 5) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k *
		     5 + 1) * cc_dim2 + 1) * cc_dim1] + tr11 * 2. * cc[m + (*
		    ido + (k * 5 + 2) * cc_dim2) * cc_dim1] + tr12 * 2. * cc[
		    m + (*ido + (k * 5 + 4) * cc_dim2) * cc_dim1] + (ti11 * 
		    2. * cc[m + ((k * 5 + 3) * cc_dim2 + 1) * cc_dim1] + ti12 
		    * 2. * cc[m + ((k * 5 + 5) * cc_dim2 + 1) * cc_dim1]);
/* L1001: */
	}
/* L101: */
    }
    if (*ido == 1) {
	return 0;
    }
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k * 5 + 1) * cc_dim2) * cc_dim1] + (cc[
			m + (i__ - 1 + (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[
			m + (ic - 1 + (k * 5 + 2) * cc_dim2) * cc_dim1]) + (
			cc[m + (i__ - 1 + (k * 5 + 5) * cc_dim2) * cc_dim1] + 
			cc[m + (ic - 1 + (k * 5 + 4) * cc_dim2) * cc_dim1]);
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + (k * 5 + 1) * cc_dim2) * cc_dim1] + (cc[m + (
			i__ + (k * 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic 
			+ (k * 5 + 2) * cc_dim2) * cc_dim1]) + (cc[m + (i__ + 
			(k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 
			5 + 4) * cc_dim2) * cc_dim1]);
		ch[m + (i__ - 1 + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] =
			 wa1[i__ - 2] * (cc[m + (i__ - 1 + (k * 5 + 1) * 
			cc_dim2) * cc_dim1] + tr11 * (cc[m + (i__ - 1 + (k * 
			5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]) - (ti11 * (cc[
			m + (i__ + (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[
			m + (i__ + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 4) * cc_dim2) * cc_dim1]))) - wa1[i__ 
			- 1] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr11 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) + (ti11 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ + (k + (ch_dim3 << 1)) * ch_dim2) * ch_dim1] = 
			wa1[i__ - 2] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) *
			 cc_dim1] + tr11 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) + (ti11 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]))) + wa1[i__ - 
			1] * (cc[m + (i__ - 1 + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr11 * (cc[m + (i__ - 1 + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ - 1 + (k * 
			5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 4) * cc_dim2) * cc_dim1]) - (ti11 * (cc[m + (i__ 
			+ (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ 
			+ (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ - 1 + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = 
			wa2[i__ - 2] * (cc[m + (i__ - 1 + (k * 5 + 1) * 
			cc_dim2) * cc_dim1] + tr12 * (cc[m + (i__ - 1 + (k * 
			5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]) - (ti12 * (cc[
			m + (i__ + (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[
			m + (i__ + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 4) * cc_dim2) * cc_dim1]))) - wa2[i__ 
			- 1] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr12 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) + (ti12 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ + (k + ch_dim3 * 3) * ch_dim2) * ch_dim1] = wa2[
			i__ - 2] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr12 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) + (ti12 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]))) + wa2[i__ - 
			1] * (cc[m + (i__ - 1 + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr12 * (cc[m + (i__ - 1 + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ - 1 + (k * 
			5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 4) * cc_dim2) * cc_dim1]) - (ti12 * (cc[m + (i__ 
			+ (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ 
			+ (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ - 1 + (k + (ch_dim3 << 2)) * ch_dim2) * ch_dim1] =
			 wa3[i__ - 2] * (cc[m + (i__ - 1 + (k * 5 + 1) * 
			cc_dim2) * cc_dim1] + tr12 * (cc[m + (i__ - 1 + (k * 
			5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]) + (ti12 * (cc[
			m + (i__ + (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[
			m + (i__ + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 4) * cc_dim2) * cc_dim1]))) - wa3[i__ 
			- 1] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr12 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) - (ti12 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ + (k + (ch_dim3 << 2)) * ch_dim2) * ch_dim1] = 
			wa3[i__ - 2] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) *
			 cc_dim1] + tr12 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) - (ti12 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]))) + wa3[i__ - 
			1] * (cc[m + (i__ - 1 + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr12 * (cc[m + (i__ - 1 + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr11 * (cc[m + (i__ - 1 + (k * 
			5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 4) * cc_dim2) * cc_dim1]) + (ti12 * (cc[m + (i__ 
			+ (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 2) * cc_dim2) * cc_dim1]) - ti11 * (cc[m + (i__ 
			+ (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ - 1 + (k + ch_dim3 * 5) * ch_dim2) * ch_dim1] = 
			wa4[i__ - 2] * (cc[m + (i__ - 1 + (k * 5 + 1) * 
			cc_dim2) * cc_dim1] + tr11 * (cc[m + (i__ - 1 + (k * 
			5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]) + (ti11 * (cc[
			m + (i__ + (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[
			m + (i__ + (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + 
			(ic + (k * 5 + 4) * cc_dim2) * cc_dim1]))) - wa4[i__ 
			- 1] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr11 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) - (ti11 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1])));
		ch[m + (i__ + (k + ch_dim3 * 5) * ch_dim2) * ch_dim1] = wa4[
			i__ - 2] * (cc[m + (i__ + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr11 * (cc[m + (i__ + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ + (k * 5 + 
			5) * cc_dim2) * cc_dim1] - cc[m + (ic + (k * 5 + 4) * 
			cc_dim2) * cc_dim1]) - (ti11 * (cc[m + (i__ - 1 + (k *
			 5 + 3) * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + (k * 
			5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ - 
			1 + (k * 5 + 5) * cc_dim2) * cc_dim1] - cc[m + (ic - 
			1 + (k * 5 + 4) * cc_dim2) * cc_dim1]))) + wa4[i__ - 
			1] * (cc[m + (i__ - 1 + (k * 5 + 1) * cc_dim2) * 
			cc_dim1] + tr11 * (cc[m + (i__ - 1 + (k * 5 + 3) * 
			cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 5 + 2) * 
			cc_dim2) * cc_dim1]) + tr12 * (cc[m + (i__ - 1 + (k * 
			5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic - 1 + (k * 
			5 + 4) * cc_dim2) * cc_dim1]) + (ti11 * (cc[m + (i__ 
			+ (k * 5 + 3) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 2) * cc_dim2) * cc_dim1]) + ti12 * (cc[m + (i__ 
			+ (k * 5 + 5) * cc_dim2) * cc_dim1] + cc[m + (ic + (k 
			* 5 + 4) * cc_dim2) * cc_dim1])));
/* L1002: */
	    }
/* L102: */
	}
/* L103: */
    }
    return 0;
} /* vradb5_ */

/* Subroutine */ int vradbg_(integer *mp, integer *ido, integer *ip, integer *
	l1, integer *idl1, doublereal *cc, doublereal *c1, doublereal *c2, 
	doublereal *ch, doublereal *ch2, integer *mdimc, doublereal *wa)
{
    /* System generated locals */
    integer ch_dim1, ch_dim2, ch_dim3, ch_offset, cc_dim1, cc_dim2, cc_dim3, 
	    cc_offset, c1_dim1, c1_dim2, c1_dim3, c1_offset, c2_dim1, c2_dim2,
	     c2_offset, ch2_dim1, ch2_dim2, ch2_offset, i__1, i__2, i__3, 
	    i__4;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, j, k, l, m, j2, ic, jc, lc, ik, is;
    static doublereal dc2, ai1, ai2, ar1, ar2, ds2;
    static integer nbd;
    static doublereal dcp, arg, dsp, tpi, ar1h, ar2h;
    static integer idp2, ipp2, idij, ipph;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa;
    ch2_dim1 = *mdimc;
    ch2_dim2 = *idl1;
    ch2_offset = 1 + ch2_dim1 * (1 + ch2_dim2);
    ch2 -= ch2_offset;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    c2_dim1 = *mdimc;
    c2_dim2 = *idl1;
    c2_offset = 1 + c2_dim1 * (1 + c2_dim2);
    c2 -= c2_offset;
    c1_dim1 = *mdimc;
    c1_dim2 = *ido;
    c1_dim3 = *l1;
    c1_offset = 1 + c1_dim1 * (1 + c1_dim2 * (1 + c1_dim3));
    c1 -= c1_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *ip;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    tpi = pimach_(&c_b26) * 2.;
    arg = tpi / (doublereal) (*ip);
    dcp = cos(arg);
    dsp = sin(arg);
    idp2 = *ido + 2;
    nbd = (*ido - 1) / 2;
    ipp2 = *ip + 2;
    ipph = (*ip + 1) / 2;
    if (*ido < *l1) {
	goto L103;
    }
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1];
/* L1001: */
	    }
/* L101: */
	}
/* L102: */
    }
    goto L106;
L103:
    i__1 = *ido;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1];
/* L1004: */
	    }
/* L104: */
	}
/* L105: */
    }
L106:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	j2 = j + j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + (
			*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1] + 
			cc[m + (*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) * 
			cc_dim1];
		ch[m + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m + 
			((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) * cc_dim1] + 
			cc[m + ((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) * 
			cc_dim1];
/* L1007: */
	    }
/* L107: */
	}
/* L108: */
    }
    if (*ido == 1) {
	goto L116;
    }
    if (nbd < *l1) {
	goto L112;
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		ic = idp2 - i__;
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = cc[m + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3) *
			     cc_dim2) * cc_dim1] + cc[m + (ic - 1 + ((j << 1) 
			    - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * ch_dim1]
			     = cc[m + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3) 
			    * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((j << 1)
			     - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    cc[m + (i__ + ((j << 1) - 1 + k * cc_dim3) * 
			    cc_dim2) * cc_dim1] - cc[m + (ic + ((j << 1) - 2 
			    + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] = 
			    cc[m + (i__ + ((j << 1) - 1 + k * cc_dim3) * 
			    cc_dim2) * cc_dim1] + cc[m + (ic + ((j << 1) - 2 
			    + k * cc_dim3) * cc_dim2) * cc_dim1];
/* L1009: */
		}
/* L109: */
	    }
/* L110: */
	}
/* L111: */
    }
    goto L116;
L112:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = cc[m + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3) *
			     cc_dim2) * cc_dim1] + cc[m + (ic - 1 + ((j << 1) 
			    - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * ch_dim1]
			     = cc[m + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3) 
			    * cc_dim2) * cc_dim1] - cc[m + (ic - 1 + ((j << 1)
			     - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    cc[m + (i__ + ((j << 1) - 1 + k * cc_dim3) * 
			    cc_dim2) * cc_dim1] - cc[m + (ic + ((j << 1) - 2 
			    + k * cc_dim3) * cc_dim2) * cc_dim1];
		    ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] = 
			    cc[m + (i__ + ((j << 1) - 1 + k * cc_dim3) * 
			    cc_dim2) * cc_dim1] + cc[m + (ic + ((j << 1) - 2 
			    + k * cc_dim3) * cc_dim2) * cc_dim1];
/* L1013: */
		}
/* L113: */
	    }
/* L114: */
	}
/* L115: */
    }
L116:
    ar1 = 1.;
    ai1 = 0.;
    i__1 = ipph;
    for (l = 2; l <= i__1; ++l) {
	lc = ipp2 - l;
	ar1h = dcp * ar1 - dsp * ai1;
	ai1 = dcp * ai1 + dsp * ar1;
	ar1 = ar1h;
	i__2 = *idl1;
	for (ik = 1; ik <= i__2; ++ik) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		c2[m + (ik + l * c2_dim2) * c2_dim1] = ch2[m + (ik + ch2_dim2)
			 * ch2_dim1] + ar1 * ch2[m + (ik + (ch2_dim2 << 1)) * 
			ch2_dim1];
		c2[m + (ik + lc * c2_dim2) * c2_dim1] = ai1 * ch2[m + (ik + *
			ip * ch2_dim2) * ch2_dim1];
/* L1017: */
	    }
/* L117: */
	}
	dc2 = ar1;
	ds2 = ai1;
	ar2 = ar1;
	ai2 = ai1;
	i__2 = ipph;
	for (j = 3; j <= i__2; ++j) {
	    jc = ipp2 - j;
	    ar2h = dc2 * ar2 - ds2 * ai2;
	    ai2 = dc2 * ai2 + ds2 * ar2;
	    ar2 = ar2h;
	    i__3 = *idl1;
	    for (ik = 1; ik <= i__3; ++ik) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    c2[m + (ik + l * c2_dim2) * c2_dim1] += ar2 * ch2[m + (ik 
			    + j * ch2_dim2) * ch2_dim1];
		    c2[m + (ik + lc * c2_dim2) * c2_dim1] += ai2 * ch2[m + (
			    ik + jc * ch2_dim2) * ch2_dim1];
/* L1018: */
		}
/* L118: */
	    }
/* L119: */
	}
/* L120: */
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	i__2 = *idl1;
	for (ik = 1; ik <= i__2; ++ik) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch2[m + (ik + ch2_dim2) * ch2_dim1] += ch2[m + (ik + j * 
			ch2_dim2) * ch2_dim1];
/* L1021: */
	    }
/* L121: */
	}
/* L122: */
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m + (
			(k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] - c1[m + (
			(k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1];
		ch[m + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m + 
			((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] + c1[m + 
			((k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1];
/* L1023: */
	    }
/* L123: */
	}
/* L124: */
    }
    if (*ido == 1) {
	goto L132;
    }
    if (nbd < *l1) {
	goto L128;
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) *
			     c1_dim1] - c1[m + (i__ + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * ch_dim1]
			     = c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) 
			    * c1_dim1] + c1[m + (i__ + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * 
			    c1_dim1] + c1[m + (i__ - 1 + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] = 
			    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * 
			    c1_dim1] - c1[m + (i__ - 1 + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
/* L1025: */
		}
/* L125: */
	    }
/* L126: */
	}
/* L127: */
    }
    goto L132;
L128:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) *
			     c1_dim1] - c1[m + (i__ + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * ch_dim1]
			     = c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) 
			    * c1_dim1] + c1[m + (i__ + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * 
			    c1_dim1] + c1[m + (i__ - 1 + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] = 
			    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * 
			    c1_dim1] - c1[m + (i__ - 1 + (k + jc * c1_dim3) * 
			    c1_dim2) * c1_dim1];
/* L1029: */
		}
/* L129: */
	    }
/* L130: */
	}
/* L131: */
    }
L132:
    if (*ido == 1) {
	return 0;
    }
    i__1 = *idl1;
    for (ik = 1; ik <= i__1; ++ik) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    c2[m + (ik + c2_dim2) * c2_dim1] = ch2[m + (ik + ch2_dim2) * 
		    ch2_dim1];
/* L1033: */
	}
/* L133: */
    }
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		c1[m + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m + (
			(k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
/* L1034: */
	    }
/* L134: */
	}
/* L135: */
    }
    if (nbd > *l1) {
	goto L139;
    }
    is = -(*ido);
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	is += *ido;
	idij = is;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    idij += 2;
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1] 
			    = wa[idij - 1] * ch[m + (i__ - 1 + (k + j * 
			    ch_dim3) * ch_dim2) * ch_dim1] - wa[idij] * ch[m 
			    + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] = 
			    wa[idij - 1] * ch[m + (i__ + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] + wa[idij] * ch[m + (i__ - 1 
			    + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
/* L1036: */
		}
/* L136: */
	    }
/* L137: */
	}
/* L138: */
    }
    goto L143;
L139:
    is = -(*ido);
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	is += *ido;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    idij = is;
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		idij += 2;
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1] 
			    = wa[idij - 1] * ch[m + (i__ - 1 + (k + j * 
			    ch_dim3) * ch_dim2) * ch_dim1] - wa[idij] * ch[m 
			    + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] = 
			    wa[idij - 1] * ch[m + (i__ + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] + wa[idij] * ch[m + (i__ - 1 
			    + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
/* L1040: */
		}
/* L140: */
	    }
/* L141: */
	}
/* L142: */
    }
L143:
    return 0;
} /* vradbg_ */

/* Subroutine */ int vradf2_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1)
{
    /* System generated locals */
    integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_dim3, cc_offset,
	     i__1, i__2, i__3;

    /* Local variables */
    static integer i__, k, m, ic, idp2;


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * 3);
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *l1;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (((k << 1) + 1) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1];
	    ch[m + (*ido + ((k << 1) + 2) * ch_dim2) * ch_dim1] = cc[m + ((k 
		    + cc_dim3) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1];
/* L1001: */
	}
/* L101: */
    }
    if ((i__1 = *ido - 2) < 0) {
	goto L107;
    } else if (i__1 == 0) {
	goto L105;
    } else {
	goto L102;
    }
L102:
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ + ((k << 1) + 1) * ch_dim2) * ch_dim1] = cc[m + (
			i__ + (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa1[i__ 
			- 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) *
			 cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1]);
		ch[m + (ic + ((k << 1) + 2) * ch_dim2) * ch_dim1] = wa1[i__ - 
			2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - cc[m + (i__ + (
			k + cc_dim3) * cc_dim2) * cc_dim1];
		ch[m + (i__ - 1 + ((k << 1) + 1) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + (
			wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) 
			* cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (
			k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1]);
		ch[m + (ic - 1 + ((k << 1) + 2) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] - (
			wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) 
			* cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (
			k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1]);
/* L1003: */
	    }
/* L103: */
	}
/* L104: */
    }
    if (*ido % 2 == 1) {
	return 0;
    }
L105:
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (((k << 1) + 2) * ch_dim2 + 1) * ch_dim1] = -cc[m + (*ido 
		    + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1];
	    ch[m + (*ido + ((k << 1) + 1) * ch_dim2) * ch_dim1] = cc[m + (*
		    ido + (k + cc_dim3) * cc_dim2) * cc_dim1];
/* L1006: */
	}
/* L106: */
    }
L107:
    return 0;
} /* vradf2_ */

/* Subroutine */ int vradf3_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2)
{
    /* System generated locals */
    integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_dim3, cc_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, k, m, ic;
    static doublereal arg;
    static integer idp2;
    static doublereal taui, taur;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_offset = 1 + ch_dim1 * (1 + (ch_dim2 << 2));
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *l1;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    arg = pimach_(&c_b26) * 2. / 3.;
    taur = cos(arg);
    taui = sin(arg);
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k * 3 + 1) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + (cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + 
		    cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + ((k * 3 + 3) * ch_dim2 + 1) * ch_dim1] = taui * (cc[m + ((
		    k + cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + 
		    (cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + (k * 3 + 2) * ch_dim2) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + taur * (cc[m + ((k + 
		    (cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + 
		    cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
/* L1001: */
	}
/* L101: */
    }
    if (*ido == 1) {
	return 0;
    }
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k * 3 + 1) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa1[
			i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k 
			+ (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[i__ - 
			2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 *
			 3) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k * 3 + 1) * ch_dim2) * ch_dim1] = cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa1[i__ - 2] 
			* cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1]
			 - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) 
			* cc_dim2) * cc_dim1]));
		ch[m + (i__ - 1 + (k * 3 + 3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + taur *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1])) + taui * (wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa2[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (ic - 1 + (k * 3 + 2) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + taur *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1])) - taui * (wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa2[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k * 3 + 3) * ch_dim2) * ch_dim1] = cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + taur * (wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1])) + taui * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1]));
		ch[m + (ic + (k * 3 + 2) * ch_dim2) * ch_dim1] = taui * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1])) - (cc[m + (i__ + (k + 
			cc_dim3) * cc_dim2) * cc_dim1] + taur * (wa1[i__ - 2] 
			* cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa2[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1]
			 - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) 
			* cc_dim2) * cc_dim1])));
/* L1002: */
	    }
/* L102: */
	}
/* L103: */
    }
    return 0;
} /* vradf3_ */

/* Subroutine */ int vradf4_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2, doublereal *wa3)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_dim3, cc_offset, ch_dim1, ch_dim2, ch_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, k, m, ic, idp2;
    static doublereal hsqt2;


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa3;
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * 5);
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *l1;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    hsqt2 = sqrt(2.) / 2.;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (((k << 2) + 1) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] + (cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + cc_dim3 
		    * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + ((k << 2) + 4) * ch_dim2) * ch_dim1] = cc[m + ((k 
		    + cc_dim3) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + 
		    cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1] - (cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + ((k << 2) + 2) * ch_dim2) * ch_dim1] = cc[m + ((k 
		    + cc_dim3) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + 
		    cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1];
	    ch[m + (((k << 2) + 3) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k + (
		    cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1];
/* L1001: */
	}
/* L101: */
    }
    if ((i__1 = *ido - 2) < 0) {
	goto L107;
    } else if (i__1 == 0) {
	goto L105;
    } else {
	goto L102;
    }
L102:
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + ((k << 2) + 1) * ch_dim2) * ch_dim1] = wa1[
			i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k 
			+ (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa3[i__ - 
			2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1]) + (cc[m + (i__ 
			- 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa2[i__ 
			- 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (ic - 1 + ((k << 2) + 4) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + (
			wa2[i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1]) - (wa1[i__ - 2] 
			* cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] * cc[m + 
			(i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] 
			+ wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 << 2)) * 
			cc_dim2) * cc_dim1]));
		ch[m + (i__ + ((k << 2) + 1) * ch_dim2) * ch_dim1] = wa1[i__ 
			- 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) *
			 cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] *
			 cc[m + (i__ + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1]) + (cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa2[i__ - 2] 
			* cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (ic + ((k << 2) + 4) * ch_dim2) * ch_dim1] = wa1[i__ - 
			2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] *
			 cc[m + (i__ + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1]) - (cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa2[i__ - 2] 
			* cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (i__ - 1 + ((k << 2) + 3) * ch_dim2) * ch_dim1] = wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa3[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 2)) * 
			cc_dim2) * cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1]) + (cc[m 
			+ (i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] - (
			wa2[i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (ic - 1 + ((k << 2) + 2) * ch_dim2) * ch_dim1] = cc[m 
			+ (i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] - (
			wa2[i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1]) - (wa1[i__ - 2] 
			* cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa3[i__ - 2] *
			 cc[m + (i__ + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1]));
		ch[m + (i__ + ((k << 2) + 3) * ch_dim2) * ch_dim1] = wa3[i__ 
			- 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * 
			cc_dim2) * cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k 
			+ (cc_dim3 << 2)) * cc_dim2) * cc_dim1] - (wa1[i__ - 
			2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2)
			 * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1]) + (cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] 
			* cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1]));
		ch[m + (ic + ((k << 2) + 2) * ch_dim2) * ch_dim1] = wa3[i__ - 
			2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1]) - (cc[m + (i__ + (k + 
			cc_dim3) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] * cc[m 
			+ (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - 
			wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1]));
/* L1003: */
	    }
/* L103: */
	}
/* L104: */
    }
    if (*ido % 2 == 1) {
	return 0;
    }
L105:
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + (*ido + ((k << 2) + 1) * ch_dim2) * ch_dim1] = hsqt2 * (cc[
		    m + (*ido + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - 
		    cc[m + (*ido + (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1])
		     + cc[m + (*ido + (k + cc_dim3) * cc_dim2) * cc_dim1];
	    ch[m + (*ido + ((k << 2) + 3) * ch_dim2) * ch_dim1] = cc[m + (*
		    ido + (k + cc_dim3) * cc_dim2) * cc_dim1] - hsqt2 * (cc[m 
		    + (*ido + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - cc[
		    m + (*ido + (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1]);
	    ch[m + (((k << 2) + 2) * ch_dim2 + 1) * ch_dim1] = -hsqt2 * (cc[m 
		    + (*ido + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + cc[
		    m + (*ido + (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1]) - 
		    cc[m + (*ido + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1];
	    ch[m + (((k << 2) + 4) * ch_dim2 + 1) * ch_dim1] = -hsqt2 * (cc[m 
		    + (*ido + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + cc[
		    m + (*ido + (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1]) + 
		    cc[m + (*ido + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1];
/* L1006: */
	}
/* L106: */
    }
L107:
    return 0;
} /* vradf4_ */

/* Subroutine */ int vradf5_(integer *mp, integer *ido, integer *l1, 
	doublereal *cc, doublereal *ch, integer *mdimc, doublereal *wa1, 
	doublereal *wa2, doublereal *wa3, doublereal *wa4)
{
    /* System generated locals */
    integer cc_dim1, cc_dim2, cc_dim3, cc_offset, ch_dim1, ch_dim2, ch_offset,
	     i__1, i__2, i__3;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, k, m, ic;
    static doublereal arg, ti11, ti12, tr11, tr12;
    static integer idp2;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa4;
    --wa3;
    --wa2;
    --wa1;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * 6);
    ch -= ch_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *l1;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    arg = pimach_(&c_b26) * 2. / 5.;
    tr11 = cos(arg);
    ti11 = sin(arg);
    tr12 = cos(arg * 2.);
    ti12 = sin(arg * 2.);
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch[m + ((k * 5 + 1) * ch_dim2 + 1) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + (cc[m + ((k + 
		    cc_dim3 * 5) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]) + (cc[m + ((k + 
		    (cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + 
		    cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + (k * 5 + 2) * ch_dim2) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + tr11 * (cc[m + ((k + 
		    cc_dim3 * 5) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]) + tr12 * (cc[m 
		    + ((k + (cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] + cc[m 
		    + ((k + cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + ((k * 5 + 3) * ch_dim2 + 1) * ch_dim1] = ti11 * (cc[m + ((
		    k + cc_dim3 * 5) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + 
		    (cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]) + ti12 * (cc[m 
		    + ((k + (cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] - cc[m 
		    + ((k + cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + (*ido + (k * 5 + 4) * ch_dim2) * ch_dim1] = cc[m + ((k + 
		    cc_dim3) * cc_dim2 + 1) * cc_dim1] + tr12 * (cc[m + ((k + 
		    cc_dim3 * 5) * cc_dim2 + 1) * cc_dim1] + cc[m + ((k + (
		    cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]) + tr11 * (cc[m 
		    + ((k + (cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] + cc[m 
		    + ((k + cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
	    ch[m + ((k * 5 + 5) * ch_dim2 + 1) * ch_dim1] = ti12 * (cc[m + ((
		    k + cc_dim3 * 5) * cc_dim2 + 1) * cc_dim1] - cc[m + ((k + 
		    (cc_dim3 << 1)) * cc_dim2 + 1) * cc_dim1]) - ti11 * (cc[m 
		    + ((k + (cc_dim3 << 2)) * cc_dim2 + 1) * cc_dim1] - cc[m 
		    + ((k + cc_dim3 * 3) * cc_dim2 + 1) * cc_dim1]);
/* L1001: */
	}
/* L101: */
    }
    if (*ido == 1) {
	return 0;
    }
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + (i__ - 1 + (k * 5 + 1) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa1[
			i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k 
			+ (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[i__ - 
			2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * cc_dim2) * 
			cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k + cc_dim3 *
			 5) * cc_dim2) * cc_dim1])) + (wa2[i__ - 2] * cc[m + (
			i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] + 
			wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + (wa3[i__ - 2] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] + wa3[
			i__ - 1] * cc[m + (i__ + (k + (cc_dim3 << 2)) * 
			cc_dim2) * cc_dim1]));
		ch[m + (i__ + (k * 5 + 1) * ch_dim2) * ch_dim1] = cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + (wa1[i__ - 2] 
			* cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) * cc_dim1]
			 - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) 
			* cc_dim2) * cc_dim1])) + (wa2[i__ - 2] * cc[m + (i__ 
			+ (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - wa2[i__ - 
			1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] + (wa3[i__ - 2] * cc[m + (i__ + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 1] * 
			cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1]));
		ch[m + (i__ - 1 + (k * 5 + 3) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + tr11 *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1]) + tr12 * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] + wa3[i__ - 2] * 
			cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1]) + ti11 * (wa1[i__ - 2] *
			 cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa4[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) * cc_dim1]
			 - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) 
			* cc_dim2) * cc_dim1])) + ti12 * (wa2[i__ - 2] * cc[m 
			+ (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - 
			wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] - (wa3[i__ - 2] * cc[m + (i__ + (
			k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 
			1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1]));
		ch[m + (ic - 1 + (k * 5 + 2) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + tr11 *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1]) + tr12 * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] + wa3[i__ - 2] * 
			cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1]) - (ti11 * (wa1[i__ - 2] 
			* cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa4[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) * cc_dim1]
			 - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) 
			* cc_dim2) * cc_dim1])) + ti12 * (wa2[i__ - 2] * cc[m 
			+ (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - 
			wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] - (wa3[i__ - 2] * cc[m + (i__ + (
			k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 
			1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1])));
		ch[m + (i__ + (k * 5 + 3) * ch_dim2) * ch_dim1] = cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + tr11 * (wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) 
			* cc_dim1] - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 5) * cc_dim2) * cc_dim1])) + tr12 * (wa2[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] * 
			cc[m + (i__ + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1])) + (ti11 * (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1])) + ti12 * (wa3[i__ - 2] 
			* cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] * cc[m + 
			(i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] + 
			wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1])));
		ch[m + (ic + (k * 5 + 2) * ch_dim2) * ch_dim1] = ti11 * (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1])) + ti12 * (wa3[i__ - 2] 
			* cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] * cc[m + 
			(i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] + 
			wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1])) - (cc[m + (i__ + (k + cc_dim3) *
			 cc_dim2) * cc_dim1] + tr11 * (wa1[i__ - 2] * cc[m + (
			i__ + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - 
			wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) 
			* cc_dim2) * cc_dim1] + (wa4[i__ - 2] * cc[m + (i__ + 
			(k + cc_dim3 * 5) * cc_dim2) * cc_dim1] - wa4[i__ - 1]
			 * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * cc_dim2) * 
			cc_dim1])) + tr12 * (wa2[i__ - 2] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] - wa2[i__ - 1] * 
			cc[m + (i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] + (wa3[i__ - 2] * cc[m + (i__ + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 1] * 
			cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1])));
		ch[m + (i__ - 1 + (k * 5 + 5) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + tr12 *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1])) + tr11 * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1])) + (ti12 * (wa1[i__ - 2]
			 * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa4[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) * cc_dim1]
			 - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) 
			* cc_dim2) * cc_dim1])) - ti11 * (wa2[i__ - 2] * cc[m 
			+ (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - 
			wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] - (wa3[i__ - 2] * cc[m + (i__ + (
			k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 
			1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1])));
		ch[m + (ic - 1 + (k * 5 + 4) * ch_dim2) * ch_dim1] = cc[m + (
			i__ - 1 + (k + cc_dim3) * cc_dim2) * cc_dim1] + tr12 *
			 (wa1[i__ - 2] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)
			) * cc_dim2) * cc_dim1] + wa1[i__ - 1] * cc[m + (i__ 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1])) + tr11 * (wa2[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] + wa2[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1])) - (ti12 * (wa1[i__ - 2]
			 * cc[m + (i__ + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 1)) * cc_dim2) * cc_dim1] - (wa4[i__ - 2] *
			 cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) * cc_dim1]
			 - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) 
			* cc_dim2) * cc_dim1])) - ti11 * (wa2[i__ - 2] * cc[m 
			+ (i__ + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] - 
			wa2[i__ - 1] * cc[m + (i__ - 1 + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1] - (wa3[i__ - 2] * cc[m + (i__ + (
			k + (cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 
			1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2)
			 * cc_dim1])));
		ch[m + (i__ + (k * 5 + 5) * ch_dim2) * ch_dim1] = cc[m + (i__ 
			+ (k + cc_dim3) * cc_dim2) * cc_dim1] + tr12 * (wa1[
			i__ - 2] * cc[m + (i__ + (k + (cc_dim3 << 1)) * 
			cc_dim2) * cc_dim1] - wa1[i__ - 1] * cc[m + (i__ - 1 
			+ (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] + (wa4[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 5) * cc_dim2) 
			* cc_dim1] - wa4[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 5) * cc_dim2) * cc_dim1])) + tr11 * (wa2[
			i__ - 2] * cc[m + (i__ + (k + cc_dim3 * 3) * cc_dim2) 
			* cc_dim1] - wa2[i__ - 1] * cc[m + (i__ - 1 + (k + 
			cc_dim3 * 3) * cc_dim2) * cc_dim1] + (wa3[i__ - 2] * 
			cc[m + (i__ + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] - wa3[i__ - 1] * cc[m + (i__ - 1 + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1])) + (ti12 * (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1])) - ti11 * (wa3[i__ - 2] 
			* cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] * cc[m + 
			(i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] + 
			wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1])));
		ch[m + (ic + (k * 5 + 4) * ch_dim2) * ch_dim1] = ti12 * (wa4[
			i__ - 2] * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * 
			cc_dim2) * cc_dim1] + wa4[i__ - 1] * cc[m + (i__ + (k 
			+ cc_dim3 * 5) * cc_dim2) * cc_dim1] - (wa1[i__ - 2] *
			 cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) * cc_dim2) * 
			cc_dim1] + wa1[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 1)) * cc_dim2) * cc_dim1])) - ti11 * (wa3[i__ - 2] 
			* cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1] + wa3[i__ - 1] * cc[m + (i__ + (k + (cc_dim3 
			<< 2)) * cc_dim2) * cc_dim1] - (wa2[i__ - 2] * cc[m + 
			(i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * cc_dim1] + 
			wa2[i__ - 1] * cc[m + (i__ + (k + cc_dim3 * 3) * 
			cc_dim2) * cc_dim1])) - (cc[m + (i__ + (k + cc_dim3) *
			 cc_dim2) * cc_dim1] + tr12 * (wa1[i__ - 2] * cc[m + (
			i__ + (k + (cc_dim3 << 1)) * cc_dim2) * cc_dim1] - 
			wa1[i__ - 1] * cc[m + (i__ - 1 + (k + (cc_dim3 << 1)) 
			* cc_dim2) * cc_dim1] + (wa4[i__ - 2] * cc[m + (i__ + 
			(k + cc_dim3 * 5) * cc_dim2) * cc_dim1] - wa4[i__ - 1]
			 * cc[m + (i__ - 1 + (k + cc_dim3 * 5) * cc_dim2) * 
			cc_dim1])) + tr11 * (wa2[i__ - 2] * cc[m + (i__ + (k 
			+ cc_dim3 * 3) * cc_dim2) * cc_dim1] - wa2[i__ - 1] * 
			cc[m + (i__ - 1 + (k + cc_dim3 * 3) * cc_dim2) * 
			cc_dim1] + (wa3[i__ - 2] * cc[m + (i__ + (k + (
			cc_dim3 << 2)) * cc_dim2) * cc_dim1] - wa3[i__ - 1] * 
			cc[m + (i__ - 1 + (k + (cc_dim3 << 2)) * cc_dim2) * 
			cc_dim1])));
/* L1002: */
	    }
/* L102: */
	}
/* L103: */
    }
    return 0;
} /* vradf5_ */

/* Subroutine */ int vradfg_(integer *mp, integer *ido, integer *ip, integer *
	l1, integer *idl1, doublereal *cc, doublereal *c1, doublereal *c2, 
	doublereal *ch, doublereal *ch2, integer *mdimc, doublereal *wa)
{
    /* System generated locals */
    integer ch_dim1, ch_dim2, ch_dim3, ch_offset, cc_dim1, cc_dim2, cc_dim3, 
	    cc_offset, c1_dim1, c1_dim2, c1_dim3, c1_offset, c2_dim1, c2_dim2,
	     c2_offset, ch2_dim1, ch2_dim2, ch2_offset, i__1, i__2, i__3, 
	    i__4;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, j, k, l, m, j2, ic, jc, lc, ik, is;
    static doublereal dc2, ai1, ai2, ar1, ar2, ds2;
    static integer nbd;
    static doublereal dcp, arg, dsp, tpi, ar1h, ar2h;
    static integer idp2, ipp2, idij, ipph;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa;
    ch2_dim1 = *mdimc;
    ch2_dim2 = *idl1;
    ch2_offset = 1 + ch2_dim1 * (1 + ch2_dim2);
    ch2 -= ch2_offset;
    ch_dim1 = *mdimc;
    ch_dim2 = *ido;
    ch_dim3 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
    ch -= ch_offset;
    c2_dim1 = *mdimc;
    c2_dim2 = *idl1;
    c2_offset = 1 + c2_dim1 * (1 + c2_dim2);
    c2 -= c2_offset;
    c1_dim1 = *mdimc;
    c1_dim2 = *ido;
    c1_dim3 = *l1;
    c1_offset = 1 + c1_dim1 * (1 + c1_dim2 * (1 + c1_dim3));
    c1 -= c1_offset;
    cc_dim1 = *mdimc;
    cc_dim2 = *ido;
    cc_dim3 = *ip;
    cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
    cc -= cc_offset;

    /* Function Body */
    tpi = pimach_(&c_b26) * 2.;
    arg = tpi / (doublereal) (*ip);
    dcp = cos(arg);
    dsp = sin(arg);
    ipph = (*ip + 1) / 2;
    ipp2 = *ip + 2;
    idp2 = *ido + 2;
    nbd = (*ido - 1) / 2;
    if (*ido == 1) {
	goto L119;
    }
    i__1 = *idl1;
    for (ik = 1; ik <= i__1; ++ik) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    ch2[m + (ik + ch2_dim2) * ch2_dim1] = c2[m + (ik + c2_dim2) * 
		    c2_dim1];
/* L1001: */
	}
/* L101: */
    }
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch[m + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m + (
			(k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1];
/* L1002: */
	    }
/* L102: */
	}
/* L103: */
    }
    if (nbd > *l1) {
	goto L107;
    }
    is = -(*ido);
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	is += *ido;
	idij = is;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    idij += 2;
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = wa[idij - 1] * c1[m + (i__ - 1 + (k + j * 
			    c1_dim3) * c1_dim2) * c1_dim1] + wa[idij] * c1[m 
			    + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    wa[idij - 1] * c1[m + (i__ + (k + j * c1_dim3) * 
			    c1_dim2) * c1_dim1] - wa[idij] * c1[m + (i__ - 1 
			    + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
/* L1004: */
		}
/* L104: */
	    }
/* L105: */
	}
/* L106: */
    }
    goto L111;
L107:
    is = -(*ido);
    i__1 = *ip;
    for (j = 2; j <= i__1; ++j) {
	is += *ido;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    idij = is;
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		idij += 2;
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1] 
			    = wa[idij - 1] * c1[m + (i__ - 1 + (k + j * 
			    c1_dim3) * c1_dim2) * c1_dim1] + wa[idij] * c1[m 
			    + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
		    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] = 
			    wa[idij - 1] * c1[m + (i__ + (k + j * c1_dim3) * 
			    c1_dim2) * c1_dim1] - wa[idij] * c1[m + (i__ - 1 
			    + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
/* L1008: */
		}
/* L108: */
	    }
/* L109: */
	}
/* L110: */
    }
L111:
    if (nbd < *l1) {
	goto L115;
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1] 
			    = ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) *
			     ch_dim1] + ch[m + (i__ - 1 + (k + jc * ch_dim3) *
			     ch_dim2) * ch_dim1];
		    c1[m + (i__ - 1 + (k + jc * c1_dim3) * c1_dim2) * c1_dim1]
			     = ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] = 
			    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] + ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + jc * c1_dim3) * c1_dim2) * c1_dim1] = 
			    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1];
/* L1012: */
		}
/* L112: */
	    }
/* L113: */
	}
/* L114: */
    }
    goto L121;
L115:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    c1[m + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1] 
			    = ch[m + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) *
			     ch_dim1] + ch[m + (i__ - 1 + (k + jc * ch_dim3) *
			     ch_dim2) * ch_dim1];
		    c1[m + (i__ - 1 + (k + jc * c1_dim3) * c1_dim2) * c1_dim1]
			     = ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] = 
			    ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] + ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    c1[m + (i__ + (k + jc * c1_dim3) * c1_dim2) * c1_dim1] = 
			    ch[m + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1];
/* L1016: */
		}
/* L116: */
	    }
/* L117: */
	}
/* L118: */
    }
    goto L121;
L119:
    i__1 = *idl1;
    for (ik = 1; ik <= i__1; ++ik) {
	i__2 = *mp;
	for (m = 1; m <= i__2; ++m) {
	    c2[m + (ik + c2_dim2) * c2_dim1] = ch2[m + (ik + ch2_dim2) * 
		    ch2_dim1];
/* L1020: */
	}
/* L120: */
    }
L121:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		c1[m + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m + (
			(k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] + ch[m + (
			(k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1];
		c1[m + ((k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m + 
			((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] - ch[m 
			+ ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
/* L1022: */
	    }
/* L122: */
	}
/* L123: */
    }

    ar1 = 1.;
    ai1 = 0.;
    i__1 = ipph;
    for (l = 2; l <= i__1; ++l) {
	lc = ipp2 - l;
	ar1h = dcp * ar1 - dsp * ai1;
	ai1 = dcp * ai1 + dsp * ar1;
	ar1 = ar1h;
	i__2 = *idl1;
	for (ik = 1; ik <= i__2; ++ik) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch2[m + (ik + l * ch2_dim2) * ch2_dim1] = c2[m + (ik + 
			c2_dim2) * c2_dim1] + ar1 * c2[m + (ik + (c2_dim2 << 
			1)) * c2_dim1];
		ch2[m + (ik + lc * ch2_dim2) * ch2_dim1] = ai1 * c2[m + (ik + 
			*ip * c2_dim2) * c2_dim1];
/* L1024: */
	    }
/* L124: */
	}
	dc2 = ar1;
	ds2 = ai1;
	ar2 = ar1;
	ai2 = ai1;
	i__2 = ipph;
	for (j = 3; j <= i__2; ++j) {
	    jc = ipp2 - j;
	    ar2h = dc2 * ar2 - ds2 * ai2;
	    ai2 = dc2 * ai2 + ds2 * ar2;
	    ar2 = ar2h;
	    i__3 = *idl1;
	    for (ik = 1; ik <= i__3; ++ik) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    ch2[m + (ik + l * ch2_dim2) * ch2_dim1] += ar2 * c2[m + (
			    ik + j * c2_dim2) * c2_dim1];
		    ch2[m + (ik + lc * ch2_dim2) * ch2_dim1] += ai2 * c2[m + (
			    ik + jc * c2_dim2) * c2_dim1];
/* L1025: */
		}
/* L125: */
	    }
/* L126: */
	}
/* L127: */
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	i__2 = *idl1;
	for (ik = 1; ik <= i__2; ++ik) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		ch2[m + (ik + ch2_dim2) * ch2_dim1] += c2[m + (ik + j * 
			c2_dim2) * c2_dim1];
/* L1028: */
	    }
/* L128: */
	}
/* L129: */
    }

    if (*ido < *l1) {
	goto L132;
    }
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		cc[m + (i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1] = ch[m 
			+ (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1];
/* L1030: */
	    }
/* L130: */
	}
/* L131: */
    }
    goto L135;
L132:
    i__1 = *ido;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		cc[m + (i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1] = ch[m 
			+ (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1];
/* L1033: */
	    }
/* L133: */
	}
/* L134: */
    }
L135:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	j2 = j + j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *mp;
	    for (m = 1; m <= i__3; ++m) {
		cc[m + (*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1] = 
			ch[m + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
		cc[m + ((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) * cc_dim1] = ch[
			m + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1];
/* L1036: */
	    }
/* L136: */
	}
/* L137: */
    }
    if (*ido == 1) {
	return 0;
    }
    if (nbd < *l1) {
	goto L141;
    }
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	j2 = j + j;
	i__2 = *l1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *ido;
	    for (i__ = 3; i__ <= i__3; i__ += 2) {
		ic = idp2 - i__;
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    cc[m + (i__ - 1 + (j2 - 1 + k * cc_dim3) * cc_dim2) * 
			    cc_dim1] = ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] + ch[m + (i__ - 1 + (k + jc * 
			    ch_dim3) * ch_dim2) * ch_dim1];
		    cc[m + (ic - 1 + (j2 - 2 + k * cc_dim3) * cc_dim2) * 
			    cc_dim1] = ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] - ch[m + (i__ - 1 + (k + jc * 
			    ch_dim3) * ch_dim2) * ch_dim1];
		    cc[m + (i__ + (j2 - 1 + k * cc_dim3) * cc_dim2) * cc_dim1]
			     = ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] + ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    cc[m + (ic + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1] 
			    = ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1];
/* L1038: */
		}
/* L138: */
	    }
/* L139: */
	}
/* L140: */
    }
    return 0;
L141:
    i__1 = ipph;
    for (j = 2; j <= i__1; ++j) {
	jc = ipp2 - j;
	j2 = j + j;
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    i__3 = *l1;
	    for (k = 1; k <= i__3; ++k) {
		i__4 = *mp;
		for (m = 1; m <= i__4; ++m) {
		    cc[m + (i__ - 1 + (j2 - 1 + k * cc_dim3) * cc_dim2) * 
			    cc_dim1] = ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] + ch[m + (i__ - 1 + (k + jc * 
			    ch_dim3) * ch_dim2) * ch_dim1];
		    cc[m + (ic - 1 + (j2 - 2 + k * cc_dim3) * cc_dim2) * 
			    cc_dim1] = ch[m + (i__ - 1 + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1] - ch[m + (i__ - 1 + (k + jc * 
			    ch_dim3) * ch_dim2) * ch_dim1];
		    cc[m + (i__ + (j2 - 1 + k * cc_dim3) * cc_dim2) * cc_dim1]
			     = ch[m + (i__ + (k + j * ch_dim3) * ch_dim2) * 
			    ch_dim1] + ch[m + (i__ + (k + jc * ch_dim3) * 
			    ch_dim2) * ch_dim1];
		    cc[m + (ic + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1] 
			    = ch[m + (i__ + (k + jc * ch_dim3) * ch_dim2) * 
			    ch_dim1] - ch[m + (i__ + (k + j * ch_dim3) * 
			    ch_dim2) * ch_dim1];
/* L1042: */
		}
/* L142: */
	    }
/* L143: */
	}
/* L144: */
    }
    return 0;
} /* vradfg_ */

/* Subroutine */ int vrfftb_(integer *m, integer *n, doublereal *r__, 
	doublereal *rt, integer *mdimr, doublereal *wsave)
{
    /* System generated locals */
    integer r_dim1, r_offset, rt_dim1, rt_offset;

    /* Local variables */
    extern /* Subroutine */ int vrftb1_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, doublereal *);

/* ***BEGIN PROLOGUE  VRFFTB */
/* ***DATE WRITTEN   850801   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A1 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM, */
/*             FOURIER SYNTHESIS, BACKWARD TRANSFORM, MULTIPLE SEQUENCES */
/* ***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST) */
/* ***PURPOSE  Backward real periodic transform, M sequences. */
/* ***DESCRIPTION */

/*  Subroutine VRFFTB computes the synthesis (backward transform) of a */
/*  number of real periodic sequences from their Fourier coefficients. */
/*  Specifically, for each set of independent Fourier coefficients */
/*  F(K), the corresponding real periodic sequence is computed. */

/*  The array WSAVE which is used by subroutine VRFFTB must be */
/*  initialized by calling subroutine VRFFTI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sets of coefficients. */

/*  N       the length of the sequences of coefficients to be */
/*          transformed.  The method is most efficient when N is a */
/*          product of small primes, however n may be any positive */
/*          integer. */

/*  R       areal two-dimensional array of size MDIMX x N containing the */
/*          coefficients to be transformed.  Each set of coefficients */
/*          F(K), K\0,1,..,N-1, is stored as a ROW of R.  Specifically, */
/*          the I-th set of independent Fourier coefficients is stored */

/*                R(I,1) = REAL( F(I,0) ), */

/*                R(I,2*K) = REAL( F(I,K) ) */

/*                R(I,2*K+1) = IMAG( F(I,K) ) */

/*                   for K = 1, 2, . . . , M-1, */

/*                and, when N is even, */

/*                R(I,N) = REAL( F(I,N/2) ). */

/*  RT      a real two-dimensional work array of size MDIMX x N. */

/*  MDIMR   the row (or first) dimension of the arrays R and RT exactly */
/*          as they appear in the calling program.  This parameter is */
/*          used to specify the variable dimension of these arrays. */

/*  WSAVE   a real one-dimensional work array which must be dimensioned */
/*          at least N+15.  The WSAVE array must be initialized by */
/*          calling subroutine VRFFTI.  A different WSAVE array must be */
/*          used for each different value of N.  This initialization does */
/*          not have to be repeated so long as N remains unchanged.  The */
/*          same WSAVE array may be used by VRFFTB and VRFFTB. */

/*  Output Parameters */

/*  R       contains M real periodic sequences corresponding to the given */
/*          coefficients.  Specifically, the I-th row of R contains the */
/*          real periodic sequence corresponding to the I-th set of */
/*          independent Fourier coefficients F(I,K) stored as */

/*               R(I,J) = X(I,J-1) ,   J = 1, 2, . . . , N, where */

/*               X(I,J) = SQRT(1/N)* F(I,0) + (-1)**J*F(I,N/2) */
/*                        + 2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N) */
/*                        - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  , */

/*                 when N is even, and */

/*               X(I,J) = SQRT(1/N)* F(I,0) + */
/*                        2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N) */
/*                        - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  , */

/*                 when N is odd. */

/*  WSAVE   contains results which must not be destroyed between calls */
/*          to VRFFTF or VRFFTB. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VRFFTF followed immediately by a call of */
/*           of VRFFTB will return the original sequences R.  Thus, */
/*           VRFFTB is the correctly normalized inverse of VRFFTF. */

/*  ----------------------------------------------------------------- */

/*  VRFFTB is a straightforward extension of the subprogram RFFTB to */
/*  handle M simultaneous sequences.  RFFTB was originally developed */
/*  by P. N. Swarztrauber of NCAR. */


/*              * * * * * * * * * * * * * * * * * * * * * */
/*              *                                       * */
/*              *         PROGRAM SPECIFICATIONS        * */
/*              *                                       * */
/*              * * * * * * * * * * * * * * * * * * * * * */


/*     DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15) */
/*     ARGUMENTS */

/*     LATEST          AUGUST 1, 1985 */
/*     REVISION */

/*     SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3, */
/*     REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2, */
/*                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH */

/*     SPECIAL         NONE */
/*     CONDITIONS */

/*     COMMON          NONE */
/*     BLOCKS */

/*     I/O             NONE */

/*     PRECISION       SINGLE */

/*     SPECIALIST      ROLAND SWEET */

/*     LANGUAGE        FORTRAN */

/*     HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE */
/*                     NATIONAL BUREAU OF STANDARDS (BOULDER). */

/*     ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION */
/*                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM. */

/*     PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77. */
/*                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN */
/*                     THE FUNCTION PIMACH. */

/*     REQUIRED        COS,SIN */
/*     RESIDENT */
/*     ROUTINES */


/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFTB1 */
/* ***END PROLOGUE  VRFFTB */

/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wsave;
    rt_dim1 = *mdimr;
    rt_offset = 1 + rt_dim1;
    rt -= rt_offset;
    r_dim1 = *mdimr;
    r_offset = 1 + r_dim1;
    r__ -= r_offset;

    /* Function Body */
    if (*n == 1) {
	return 0;
    }
    vrftb1_(m, n, &r__[r_offset], &rt[rt_offset], mdimr, &wsave[1], &wsave[*n 
	    + 1]);
    return 0;
} /* vrfftb_ */

/* Subroutine */ int vrfftf_(integer *m, integer *n, doublereal *r__, 
	doublereal *rt, integer *mdimr, doublereal *wsave)
{
    /* System generated locals */
    integer r_dim1, r_offset, rt_dim1, rt_offset;

    /* Local variables */
    extern /* Subroutine */ int vrftf1_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, doublereal *);

/* ***BEGIN PROLOGUE  VRFFTF */
/* ***DATE WRITTEN   850801   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A1 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM, */
/*             FOURIER ANALYSIS, FORWARD TRANSFORM, MULTIPLE SEQUENCES */
/* ***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST) */
/* ***PURPOSE  Forward real periodic transform, M sequences. */
/* ***DESCRIPTION */

/*  Subroutine VRFFTF computes the Fourier coefficients (forward */
/*  transform) of a number of real periodic sequences.  Specifically, */
/*  for each sequence the subroutine claculates the independent */
/*  Fourier coefficients described below at output parameter R. */

/*  The array WSAVE which is used by subroutine VRFFTF must be */
/*  initialized by calling subroutine VRFFTI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequences to be transformed.  The method */
/*          is most efficient when N is a product of small primes, */
/*          however n may be any positive integer. */

/*  R       areal two-dimensional array of size MDIMX x N containing the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of R.  Thus, the I-th sequence to be transformed, */
/*          X(I,J), J=0,1,...,N-1, is stored as */

/*               R(I,J) = X(I,J-1) , J=1, 2, . . . , N. */

/*  RT      a real two-dimensional work array of size MDIMX x N. */

/*  MDIMR   the row (or first) dimension of the arrays R and RT exactly */
/*          as they appear in the calling program.  This parameter is */
/*          used to specify the variable dimension of these arrays. */

/*  WSAVE   a real one-dimensional work array which must be dimensioned */
/*          at least N+15.  The WSAVE array must be initialized by */
/*          calling subroutine VRFFTI.  A different WSAVE array must be */
/*          used for each different value of N.  This initialization does */
/*          not have to be repeated so long as N remains unchanged.  The */
/*          same WSAVE array may be used by VRFFTF and VRFFTB. */

/*  Output Parameters */

/*  R       contains the Fourier coefficients F(K) for each of the M */
/*          input sequences.  Specifically, row I of R, R(I,J), */
/*          J=1,2,..,N, contains the independent Fourier coefficients */
/*          F(I,K), for the I-th input sequence stored as */

/*             R(I,1) = REAL( F(I,0) ), */
/*                    = SQRT(1/N)*SUM(J=0,N-1)[ X(I,J) ], */

/*             R(I,2*K) = REAL( F(I,K) ) */
/*                      = SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*COS(2J*K*PI/N)] */

/*             R(I,2*K+1) = IMAG( F(I,K) ) */
/*                        =-SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*SIN(2J*K*PI/N)] */

/*                   for K = 1, 2, . . . , M-1, */

/*              and, when N is even, */

/*              R(I,N) = REAL( F(I,N/2) ). */
/*                     = SQRT(1/N)*SUM(J=0,N-1)[ (-1)**J*X(I,J) ]. */

/*  WSAVE   contains results which must not be destroyed between calls */
/*          to VRFFTF or VRFFTB. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VRFFTF followed immediately by a call of */
/*           of VRFFTB will return the original sequences R.  Thus, */
/*           VRFFTB is the correctly normalized inverse of VRFFTF. */

/*  ----------------------------------------------------------------- */

/*  VRFFTF is a straightforward extension of the subprogram RFFTF to */
/*  handle M simultaneous sequences.  RFFTF was originally developed */
/*  by P. N. Swarztrauber of NCAR. */


/*              * * * * * * * * * * * * * * * * * * * * * */
/*              *                                       * */
/*              *         PROGRAM SPECIFICATIONS        * */
/*              *                                       * */
/*              * * * * * * * * * * * * * * * * * * * * * */


/*     DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15) */
/*     ARGUMENTS */

/*     LATEST          AUGUST 1, 1985 */
/*     REVISION */

/*     SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3, */
/*     REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2, */
/*                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH */

/*     SPECIAL         NONE */
/*     CONDITIONS */

/*     COMMON          NONE */
/*     BLOCKS */

/*     I/O             NONE */

/*     PRECISION       SINGLE */

/*     SPECIALIST      ROLAND SWEET */

/*     LANGUAGE        FORTRAN */

/*     HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE */
/*                     NATIONAL BUREAU OF STANDARDS (BOULDER). */

/*     ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION */
/*                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM. */

/*     PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77. */
/*                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN */
/*                     THE FUNCTION PIMACH. */

/*     REQUIRED        COS,SIN */
/*     RESIDENT */
/*     ROUTINES */


/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFTF1 */
/* ***END PROLOGUE  VRFFTF */

/*     VRFFTPK, VERSION 1, AUGUST 1985 */

/* ***FIRST EXECUTABLE STATEMENT  VRFFTF */
    /* Parameter adjustments */
    --wsave;
    rt_dim1 = *mdimr;
    rt_offset = 1 + rt_dim1;
    rt -= rt_offset;
    r_dim1 = *mdimr;
    r_offset = 1 + r_dim1;
    r__ -= r_offset;

    /* Function Body */
    if (*n == 1) {
	return 0;
    }
    vrftf1_(m, n, &r__[r_offset], &rt[rt_offset], mdimr, &wsave[1], &wsave[*n 
	    + 1]);
    return 0;
} /* vrfftf_ */

/* Subroutine */ int vrffti_(integer *n, doublereal *wsave)
{
    extern /* Subroutine */ int vrfti1_(integer *, doublereal *, doublereal *)
	    ;

/* ***BEGIN PROLOGUE  VRFFTI */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A1 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM, */
/*             MULTIPLE SEQUENCES */
/* ***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST) */
/* ***PURPOSE  Initialization for VRFFTF and VRFFTB. */
/* ***DESCRIPTION */

/*  Subroutine VRFFTI initializes the array WSAVE which is used in */
/*  both VRFFTF and VRFFTB.  The prime factorization of N together with */
/*  a tabulation of certain trigonometric functions are computed and */
/*  stored in the array WSAVE. */

/*  Input Parameter */

/*  N       the length of the sequence to be transformed.  There is no */
/*          restriction on N. */

/*  Output Parameter */

/*  WSAVE   a work array which must be dimensioned at least N+15. */
/*          The same work array can be used for both VRFFTF and VRFFTB */
/*          as long as N remains unchanged.  Different WSAVE arrays */
/*          are required for different values of N.  The contents of */
/*          WSAVE must not be changed between calls of VRFFTF or VRFFTB. */


/*              * * * * * * * * * * * * * * * * * * * * * */
/*              *                                       * */
/*              *         PROGRAM SPECIFICATIONS        * */
/*              *                                       * */
/*              * * * * * * * * * * * * * * * * * * * * * */


/*     DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15) */
/*     ARGUMENTS */

/*     LATEST          AUGUST 1, 1985 */
/*     REVISION */

/*     SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3, */
/*     REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2, */
/*                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH */

/*     SPECIAL         NONE */
/*     CONDITIONS */

/*     COMMON          NONE */
/*     BLOCKS */

/*     I/O             NONE */

/*     PRECISION       SINGLE */

/*     SPECIALIST      ROLAND SWEET */

/*     LANGUAGE        FORTRAN */

/*     HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE */
/*                     NATIONAL BUREAU OF STANDARDS (BOULDER). */

/*     ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION */
/*                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM. */

/*     PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77. */
/*                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN */
/*                     THE FUNCTION PIMACH. */

/*     REQUIRED        COS,SIN */
/*     RESIDENT */
/*     ROUTINES */


/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFTI1 */
/* ***END PROLOGUE  VRFFTI */

/*     VRFFTPK, VERSION 1, AUGUST 1985 */

/* ***FIRST EXECUTABLE STATEMENT  VRFFTI */
    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    if (*n == 1) {
	return 0;
    }
    vrfti1_(n, &wsave[1], &wsave[*n + 1]);
    return 0;
} /* vrffti_ */

/* Subroutine */ int vrftb1_(integer *m, integer *n, doublereal *c__, 
	doublereal *ch, integer *mdimc, doublereal *wa, doublereal *fac)
{
    /* System generated locals */
    integer ch_dim1, ch_offset, c_dim1, c_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, 
	    idl1;
    static doublereal scale;
    extern /* Subroutine */ int vradb2_(integer *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *), vradb3_(
	    integer *, integer *, integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, doublereal *), vradb4_(integer *, 
	    integer *, integer *, doublereal *, doublereal *, integer *, 
	    doublereal *, doublereal *, doublereal *), vradb5_(integer *, 
	    integer *, integer *, doublereal *, doublereal *, integer *, 
	    doublereal *, doublereal *, doublereal *, doublereal *), vradbg_(
	    integer *, integer *, integer *, integer *, integer *, doublereal 
	    *, doublereal *, doublereal *, doublereal *, doublereal *, 
	    integer *, doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa;
    ch_dim1 = *mdimc;
    ch_offset = 1 + ch_dim1;
    ch -= ch_offset;
    c_dim1 = *mdimc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
    --fac;

    /* Function Body */
    nf = (integer) fac[2];
    na = 0;
    l1 = 1;
    iw = 1;
    i__1 = nf;
    for (k1 = 1; k1 <= i__1; ++k1) {
	ip = (integer) fac[k1 + 2];
	l2 = ip * l1;
	ido = *n / l2;
	idl1 = ido * l1;
	if (ip != 4) {
	    goto L103;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	if (na != 0) {
	    goto L101;
	}
	vradb4_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3]);
	goto L102;
L101:
	vradb4_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3]);
L102:
	na = 1 - na;
	goto L115;
L103:
	if (ip != 2) {
	    goto L106;
	}
	if (na != 0) {
	    goto L104;
	}
	vradb2_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw]);
	goto L105;
L104:
	vradb2_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw]);
L105:
	na = 1 - na;
	goto L115;
L106:
	if (ip != 3) {
	    goto L109;
	}
	ix2 = iw + ido;
	if (na != 0) {
	    goto L107;
	}
	vradb3_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2]);
	goto L108;
L107:
	vradb3_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2]);
L108:
	na = 1 - na;
	goto L115;
L109:
	if (ip != 5) {
	    goto L112;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	ix4 = ix3 + ido;
	if (na != 0) {
	    goto L110;
	}
	vradb5_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3], &wa[ix4]);
	goto L111;
L110:
	vradb5_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3], &wa[ix4]);
L111:
	na = 1 - na;
	goto L115;
L112:
	if (na != 0) {
	    goto L113;
	}
	vradbg_(m, &ido, &ip, &l1, &idl1, &c__[c_offset], &c__[c_offset], &
		c__[c_offset], &ch[ch_offset], &ch[ch_offset], mdimc, &wa[iw])
		;
	goto L114;
L113:
	vradbg_(m, &ido, &ip, &l1, &idl1, &ch[ch_offset], &ch[ch_offset], &ch[
		ch_offset], &c__[c_offset], &c__[c_offset], mdimc, &wa[iw]);
L114:
	if (ido == 1) {
	    na = 1 - na;
	}
L115:
	l1 = l2;
	iw += (ip - 1) * ido;
/* L116: */
    }
    scale = sqrt(1. / *n);
    if (na == 0) {
	goto L118;
    }
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *m;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    c__[i__ + j * c_dim1] = scale * ch[i__ + j * ch_dim1];
/* L117: */
	}
    }
    return 0;
L118:
    i__2 = *n;
    for (j = 1; j <= i__2; ++j) {
	i__1 = *m;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    c__[i__ + j * c_dim1] = scale * c__[i__ + j * c_dim1];
/* L119: */
	}
    }
    return 0;
} /* vrftb1_ */

/* Subroutine */ int vrftf1_(integer *m, integer *n, doublereal *c__, 
	doublereal *ch, integer *mdimc, doublereal *wa, doublereal *fac)
{
    /* System generated locals */
    integer ch_dim1, ch_offset, c_dim1, c_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido,
	     idl1;
    static doublereal scale;
    extern /* Subroutine */ int vradf2_(integer *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *), vradf3_(
	    integer *, integer *, integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, doublereal *), vradf4_(integer *, 
	    integer *, integer *, doublereal *, doublereal *, integer *, 
	    doublereal *, doublereal *, doublereal *), vradf5_(integer *, 
	    integer *, integer *, doublereal *, doublereal *, integer *, 
	    doublereal *, doublereal *, doublereal *, doublereal *), vradfg_(
	    integer *, integer *, integer *, integer *, integer *, doublereal 
	    *, doublereal *, doublereal *, doublereal *, doublereal *, 
	    integer *, doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa;
    ch_dim1 = *mdimc;
    ch_offset = 1 + ch_dim1;
    ch -= ch_offset;
    c_dim1 = *mdimc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
    --fac;

    /* Function Body */
    nf = (integer) fac[2];
    na = 1;
    l2 = *n;
    iw = *n;
    i__1 = nf;
    for (k1 = 1; k1 <= i__1; ++k1) {
	kh = nf - k1;
	ip = (integer) fac[kh + 3];
	l1 = l2 / ip;
	ido = *n / l2;
	idl1 = ido * l1;
	iw -= (ip - 1) * ido;
	na = 1 - na;
	if (ip != 4) {
	    goto L102;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	if (na != 0) {
	    goto L101;
	}
	vradf4_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3]);
	goto L110;
L101:
	vradf4_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3]);
	goto L110;
L102:
	if (ip != 2) {
	    goto L104;
	}
	if (na != 0) {
	    goto L103;
	}
	vradf2_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw]);
	goto L110;
L103:
	vradf2_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw]);
	goto L110;
L104:
	if (ip != 3) {
	    goto L106;
	}
	ix2 = iw + ido;
	if (na != 0) {
	    goto L105;
	}
	vradf3_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2]);
	goto L110;
L105:
	vradf3_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2]);
	goto L110;
L106:
	if (ip != 5) {
	    goto L108;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	ix4 = ix3 + ido;
	if (na != 0) {
	    goto L107;
	}
	vradf5_(m, &ido, &l1, &c__[c_offset], &ch[ch_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3], &wa[ix4]);
	goto L110;
L107:
	vradf5_(m, &ido, &l1, &ch[ch_offset], &c__[c_offset], mdimc, &wa[iw], 
		&wa[ix2], &wa[ix3], &wa[ix4]);
	goto L110;
L108:
	if (ido == 1) {
	    na = 1 - na;
	}
	if (na != 0) {
	    goto L109;
	}
	vradfg_(m, &ido, &ip, &l1, &idl1, &c__[c_offset], &c__[c_offset], &
		c__[c_offset], &ch[ch_offset], &ch[ch_offset], mdimc, &wa[iw])
		;
	na = 1;
	goto L110;
L109:
	vradfg_(m, &ido, &ip, &l1, &idl1, &ch[ch_offset], &ch[ch_offset], &ch[
		ch_offset], &c__[c_offset], &c__[c_offset], mdimc, &wa[iw]);
	na = 0;
L110:
	l2 = l1;
/* L111: */
    }
    scale = sqrt(1. / *n);
    if (na == 1) {
	goto L113;
    }
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *m;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    c__[i__ + j * c_dim1] = scale * ch[i__ + j * ch_dim1];
/* L112: */
	}
    }
    return 0;
L113:
    i__2 = *n;
    for (j = 1; j <= i__2; ++j) {
	i__1 = *m;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    c__[i__ + j * c_dim1] = scale * c__[i__ + j * c_dim1];
/* L114: */
	}
    }
    return 0;
} /* vrftf1_ */

/* Subroutine */ int vrfti1_(integer *n, doublereal *wa, doublereal *fac)
{
    /* Initialized data */

    static integer ntryh[4] = { 4,2,3,5 };

    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static integer i__, j, k1, l1, l2, ib;
    static doublereal fi;
    static integer ld, ii, nf, ip, nl, is, nq, nr;
    static doublereal arg;
    static integer ido, ipm;
    static doublereal tpi;
    static integer nfm1;
    static doublereal argh;
    static integer ntry;
    static doublereal argld;
    extern doublereal pimach_(doublereal *);


/*     VRFFTPK, VERSION 1, AUGUST 1985 */

    /* Parameter adjustments */
    --wa;
    --fac;

    /* Function Body */
    nl = *n;
    nf = 0;
    j = 0;
L101:
    ++j;
    if (j - 4 <= 0) {
	goto L102;
    } else {
	goto L103;
    }
L102:
    ntry = ntryh[j - 1];
    goto L104;
L103:
    ntry += 2;
L104:
    nq = nl / ntry;
    nr = nl - ntry * nq;
    if (nr != 0) {
	goto L101;
    } else {
	goto L105;
    }
L105:
    ++nf;
    fac[nf + 2] = (doublereal) ntry;
    nl = nq;
    if (ntry != 2) {
	goto L107;
    }
    if (nf == 1) {
	goto L107;
    }
    i__1 = nf;
    for (i__ = 2; i__ <= i__1; ++i__) {
	ib = nf - i__ + 2;
	fac[ib + 2] = fac[ib + 1];
/* L106: */
    }
    fac[3] = 2.;
L107:
    if (nl != 1) {
	goto L104;
    }
    fac[1] = (doublereal) (*n);
    fac[2] = (doublereal) nf;
    tpi = pimach_(&c_b26) * 2.;
    argh = tpi / (doublereal) (*n);
    is = 0;
    nfm1 = nf - 1;
    l1 = 1;
    if (nfm1 == 0) {
	return 0;
    }
    i__1 = nfm1;
    for (k1 = 1; k1 <= i__1; ++k1) {
	ip = (integer) fac[k1 + 2];
	ld = 0;
	l2 = l1 * ip;
	ido = *n / l2;
	ipm = ip - 1;
	i__2 = ipm;
	for (j = 1; j <= i__2; ++j) {
	    ld += l1;
	    i__ = is;
	    argld = (doublereal) ld * argh;
	    fi = 0.;
	    i__3 = ido;
	    for (ii = 3; ii <= i__3; ii += 2) {
		i__ += 2;
		fi += 1.;
		arg = fi * argld;
		wa[i__ - 1] = cos(arg);
		wa[i__] = sin(arg);
/* L108: */
	    }
	    is += ido;
/* L109: */
	}
	l1 = l2;
/* L110: */
    }
    return 0;
} /* vrfti1_ */

/* Subroutine */ int vsinqb_(integer *m, integer *n, doublereal *x, 
	doublereal *xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Local variables */
    static integer j, k, kc, ns2;
    static doublereal xhold;
    extern /* Subroutine */ int vcosqb_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VSINQB */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Normalized inverse of VSINQF. */
/* ***DESCRIPTION */

/*  Subroutine VSINQB computes the backward fast Fourier sine transform */
/*  of M quarter wave sequences.  That is, sine series representations */
/*  with only odd wave numbers.  The transform is defined below at output */
/*  parameter X. */

/*  The array WSAVE which is used by subroutine VSINQB must be */
/*  initialized by calling subroutine VSINQI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequences to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  X       an array of size at least X(MDIMX,N) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N. */

/*  XT      a work array of size at least XT(MDIMX,N). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array which must be dimensioned at least 2*N+15 */
/*          in the program that calls VSINQB.  The WSAVE array must be */
/*          initialized by calling subroutine VSINQI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged. */

/*  Output Parameters */

/*  X       for I=1,...,N and J=1,...,M */

/*               X(I)= the sum from K=1 to K=N of */

/*                 4*X(K)*SIN((2k-1)*I*PI/(2*N)) /SQRT(4*N) */


/*  WSAVE   contains initialization calculations which must not */
/*          be destroyed between calls of VSINQB or VSINQF. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VSINQF followed immediately by a call of */
/*           of VSINQB will return the original sequences X.  Thus, */
/*           VSINQB is the correctly normalized inverse VSINQF. */

/*  ----------------------------------------------------------------- */

/*  VSINQB is a straightforward extension of the subprogram SINQB to */
/*  handle M simultaneous sequences.  SINQB was originally developed */
/*  by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VCOSQB */
/* ***END PROLOGUE  VSINQB */
/* ***FIRST EXECUTABLE STATEMENT  VSINQB */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n <= 1) {
	goto L900;
    }

/*  CASE  N .GT. 1 */

/*     ... PREPROCESSING */

    ns2 = *n / 2;
    i__1 = *n;
    for (k = 2; k <= i__1; k += 2) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + k * x_dim1] = -x[j + k * x_dim1];
/* L210: */
	}
    }

/*     ... COSINE QUARTER WAVE TRANSFORM */

    vcosqb_(m, n, &x[x_offset], &xt[xt_offset], mdimx, &wsave[1]);

/*     ... POSTPROCESSING */

    i__2 = ns2;
    for (k = 1; k <= i__2; ++k) {
	kc = *n - k;
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    xhold = x[j + k * x_dim1];
	    x[j + k * x_dim1] = x[j + (kc + 1) * x_dim1];
	    x[j + (kc + 1) * x_dim1] = xhold;
/* L220: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vsinqb_ */

/* Subroutine */ int vsinqf_(integer *m, integer *n, doublereal *x, 
	doublereal *xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Local variables */
    static integer j, k, kc, ns2;
    static doublereal xhold;
    extern /* Subroutine */ int vcosqf_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VSINQF */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Forward sine transform, odd wave numbers, M sequences. */
/* ***DESCRIPTION */

/*  Subroutine VSINQF computes the forward fast Fourier sine transform */
/*  of M quarter wave sequences.  That is, sine series representations */
/*  with only odd wave numbers.  The transform is defined below at output */
/*  parameter X. */

/*  The array WSAVE which is used by subroutine VSINQF must be */
/*  initialized by calling subroutine VSINQI(N,WSAVE). */


/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequences to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  X       an array of size at least X(MDIMX,N) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N. */

/*  XT      a work array of size at least XT(MDIMX,N). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array which must be dimensioned at least 2*N+15 */
/*          in the program that calls VSINQF.  The WSAVE array must be */
/*          initialized by calling subroutine VSINQI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged. */

/*  Output Parameters */

/*  X       For I=1,...,N and J=1,...,M */

/*               X(I) = ( (-1)**(I-1)*X(N) */

/*                        + the sum from K=1 to K=N-1 of */

/*                        2*X(K)*SIN((2*I-1)*K*PI/(2*N)) )/SQRT(4*N) */

/*  WSAVE   contains initialization calculations which must not */
/*          be destroyed between calls of VSINQF or VSINQB. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VSINQF followed immediately by a call of */
/*           of VSINQB will return the original sequences X.  Thus, */
/*           VSINQB is the correctly normalized inverse VSINQF. */

/*  ----------------------------------------------------------------- */

/*  VSINQF is a straightforward extension of the subprogram SINQF to */
/*  handle M simultaneous sequences.  SINQF was originally developed */
/*  by P. N. Swarztrauber of NCAR. */

/* ***ROUTINES CALLED  VCOSQF */
/* ***END PROLOGUE  VSINQF */
/* ***FIRST EXECUTABLE STATEMENT  VSINQF */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n <= 1) {
	goto L900;
    }

/*     ... PREPROCESSING */

    ns2 = *n / 2;
    i__1 = ns2;
    for (k = 1; k <= i__1; ++k) {
	kc = *n - k;
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xhold = x[j + k * x_dim1];
	    x[j + k * x_dim1] = x[j + (kc + 1) * x_dim1];
	    x[j + (kc + 1) * x_dim1] = xhold;
/* L110: */
	}
    }

/*     ... COSINE QUARTER WAVE TRANSFORM */

    vcosqf_(m, n, &x[x_offset], &xt[xt_offset], mdimx, &wsave[1]);

/*     ... POSTPROCESSING */

    i__2 = *n;
    for (k = 2; k <= i__2; k += 2) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + k * x_dim1] = -x[j + k * x_dim1];
/* L120: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vsinqf_ */

/* Subroutine */ int vsinqi_(integer *n, doublereal *wsave)
{
    extern /* Subroutine */ int vcosqi_(integer *, doublereal *);

/* ***BEGIN PROLOGUE  VSINQI */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE */
/*             NUMBERS, MULTIPLE SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Initialize for VSINQF and VSINQB. */
/* ***DESCRIPTION */

/*  Subroutine VSINQI initializes the array WSAVE which is used in */
/*  both VSINQF and VSINQB.  The prime factorization of N together with */
/*  a tabulation of the trigonometric functions are computed and */
/*  stored in WSAVE. */

/*  Input Parameter */

/*  N       the length of the sequence to be transformed.  The method */
/*          is most efficient when N is a product of small primes. */

/*  Output Parameter */

/*  WSAVE   a work array which must be dimensioned at least 3*N+15. */
/*          The same work array can be used for both VSINQF and VSINQB */
/*          as long as N remains unchanged.  Different WSAVE arrays */
/*          are required for different values of N.  The contents of */

/*          WSAVE must not be changed between calls of VSINQF or VSINQB. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VCOSQI */
/* ***END PROLOGUE  VSINQI */
/* ***FIRST EXECUTABLE STATEMENT  VSINQI */
    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    vcosqi_(n, &wsave[1]);
    return 0;
} /* vsinqi_ */

/* Subroutine */ int vsint_(integer *m, integer *n, doublereal *x, doublereal 
	*xt, integer *mdimx, doublereal *wsave)
{
    /* System generated locals */
    integer x_dim1, x_offset, xt_dim1, xt_offset, i__1, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j, k;
    static doublereal t1, t2;
    static integer kc, nf;
    static doublereal xh;
    static integer np1, ns2, modn;
    static doublereal scale, sqrth;
    extern /* Subroutine */ int vrfftf_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);

/* ***BEGIN PROLOGUE  VSINT */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, MULTIPLE */
/*             SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F., (NIST) */
/* ***PURPOSE  Sine transform of one or more real, odd sequences. */
/* ***DESCRIPTION */

/*  Subroutine VSINT computes the discrete Fourier sine transform */
/*  of M odd sequences X(J,I), J=1,...,M.  The transform is defined */
/*  below at output parameter X. */

/*  The array WSAVE which is used by subroutine VSINT must be */
/*  initialized by calling subroutine VSINTI(N,WSAVE). */

/*  Input Parameters */

/*  M       the number of sequences to be transformed. */

/*  N       the length of the sequence to be transformed.  The method */
/*          is most efficient when N+1 is the product of small primes. */

/*  X       an array of size at least X(MDIMX,N+1) which contains the */
/*          the sequences to be transformed.  The sequences are stored */
/*          in the ROWS of X.  Thus, the Jth sequence is stored in */
/*          X(J,I), I=1,..,N.  The extra column of X is used as work */
/*          storage. */

/*  XT      a work array of size at least XT(MDIMX,N+1). */

/*  MDIMX   the first dimension of the array X exactly as it appears in */
/*          the calling program. */

/*  WSAVE   a work array with dimension at least INT(2.5*N+15) */
/*          in the program that calls VSINT.  The WSAVE array must be */
/*          initialized by calling subroutine VSINTI(N,WSAVE), and a */
/*          different WSAVE array must be used for each different */
/*          value of N.  This initialization does not have to be */
/*          repeated so long as N remains unchanged. */

/*  Output Parameters */

/*  X       for I=1,...,N and J=1,...,M */

/*               X(J,I)= the sum from k=1 to k=N */

/*                    2*X(J,K)*SIN(K*I*PI/(N+1))/SQRT(2*(N+1)) */

/*  WSAVE   contains initialization calculations which must not be */
/*          destroyed between calls of VSINT. */

/*  ----------------------------------------------------------------- */

/*  NOTE  -  A call of VSINT followed immediately by another call */
/*           of VSINT will return the original sequences X.  Thus, */
/*           VSINT is the correctly normalized inverse of itself. */

/*  ----------------------------------------------------------------- */

/*  VSINT is a straightforward extension of the subprogram SINT to */
/*  handle M simultaneous sequences.  The scaling of the sequences */
/*  computed by VSINT is different than that of SINT.  SINT was */
/*  originally developed by P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFFTF */
/* ***END PROLOGUE  VSINT */
/* ***FIRST EXECUTABLE STATEMENT  SINT */
    /* Parameter adjustments */
    xt_dim1 = *mdimx;
    xt_offset = 1 + xt_dim1;
    xt -= xt_offset;
    x_dim1 = *mdimx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --wsave;

    /* Function Body */
    if (*m <= 0) {
	goto L900;
    }
    if (*n <= 1) {
	goto L900;
    }
    if (*n > 2) {
	goto L300;
    }

/*  CASE   N = 2 */

    sqrth = sqrt(.5);
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	xh = sqrth * (x[j + x_dim1] + x[j + (x_dim1 << 1)]);
	x[j + (x_dim1 << 1)] = sqrth * (x[j + x_dim1] - x[j + (x_dim1 << 1)]);
	x[j + x_dim1] = xh;
/* L201: */
    }
    goto L900;

/*  CASE   N .GT. 2 */

/*     ... PREPROCESSING */

L300:
    np1 = *n + 1;
    ns2 = *n / 2;
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	xt[j + xt_dim1] = 0.;
/* L301: */
    }
    i__1 = ns2;
    for (k = 1; k <= i__1; ++k) {
	kc = np1 - k;
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    t1 = x[j + k * x_dim1] - x[j + kc * x_dim1];
	    t2 = wsave[k] * (x[j + k * x_dim1] + x[j + kc * x_dim1]);
	    xt[j + (k + 1) * xt_dim1] = t1 + t2;
	    xt[j + (kc + 1) * xt_dim1] = t2 - t1;
/* L310: */
	}
    }
    modn = *n % 2;
    if (modn != 0) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    xt[j + (ns2 + 2) * xt_dim1] = x[j + (ns2 + 1) * x_dim1] * 4.;
/* L320: */
	}
    }

/*     ... REAL PERIODIC TRANSFORM */

    nf = ns2 + 1;
    vrfftf_(m, &np1, &xt[xt_offset], &x[x_offset], mdimx, &wsave[nf]);

/*     ... POSTPROCESSING */

    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	x[j + x_dim1] = xt[j + xt_dim1] * .5;
/* L330: */
    }
    i__2 = *n;
    for (i__ = 3; i__ <= i__2; i__ += 2) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + (i__ - 1) * x_dim1] = -xt[j + i__ * xt_dim1];
/* L340: */
	}
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + i__ * x_dim1] = x[j + (i__ - 2) * x_dim1] + xt[j + (i__ - 1)
		     * xt_dim1];
/* L345: */
	}
/* L350: */
    }
    if (modn == 0) {
	i__2 = *m;
	for (j = 1; j <= i__2; ++j) {
	    x[j + *n * x_dim1] = -xt[j + (*n + 1) * xt_dim1];
/* L360: */
	}
    }

/*     ... NORMALIZATION */

    scale = sqrt(.5);
    i__2 = *n;
    for (i__ = 1; i__ <= i__2; ++i__) {
	i__1 = *m;
	for (j = 1; j <= i__1; ++j) {
	    x[j + i__ * x_dim1] = scale * x[j + i__ * x_dim1];
/* L370: */
	}
    }

/*  EXIT */

L900:
    return 0;
} /* vsint_ */

/* Subroutine */ int vsinti_(integer *n, doublereal *wsave)
{
    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    double sin(doublereal);

    /* Local variables */
    static integer k, kf;
    static doublereal fk, dt, pi;
    static integer ks, np1, ns2;
    extern doublereal pimach_(doublereal *);
    extern /* Subroutine */ int vrffti_(integer *, doublereal *);

/* ***BEGIN PROLOGUE  VSINTI */
/* ***DATE WRITTEN   860701   (YYMMDD) */
/* ***REVISION DATE  900509   (YYMMDD) */
/* ***CATEGORY NO.  J1A3 */
/* ***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, MULTIPLE */
/*             SEQUENCES */
/* ***AUTHOR  BOISVERT, R. F. (NIST) */
/* ***PURPOSE  Initialize for VSINT. */
/* ***DESCRIPTION */

/*  Subroutine VSINTI initializes the array WSAVE which is used in */
/*  subroutine SINT.  The prime factorization of N together with */
/*  a tabulation of the trigonometric functions are computed and */
/*  stored in WSAVE. */

/*  Input Parameter */

/*  N       the length of the sequence to be transformed.  The method */
/*          is most efficient when N+1 is a product of small primes. */

/*  Output Parameter */

/*  WSAVE   a work array with at least INT(2.5*N+15) locations. */
/*          Different WSAVE arrays are required for different values */
/*          of N.  The contents of WSAVE must not be changed between */
/*          calls of VSINT. */

/*  ----------------------------------------------------------------- */

/*  VSINTI is a straightforward extension of the subprogram SINTI to */
/*  handle M simultaneous sequences.  SINTI was originally developed */
/*  P. N. Swarztrauber of NCAR. */

/* ***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel */
/*               Computations, (G. Rodrigue, ed.), Academic Press, 1982, */
/*               pp. 51-83. */
/* ***ROUTINES CALLED  VRFFTI */
/* ***END PROLOGUE  VSINTI */
/* ***FIRST EXECUTABLE STATEMENT  SINTI */
    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    pi = pimach_(&c_b26);
    if (*n <= 1) {
	return 0;
    }
    np1 = *n + 1;
    ns2 = *n / 2;
    dt = pi / (doublereal) np1;
    ks = 1;
    kf = ks + ns2 - 1;
    fk = 0.;
    i__1 = kf;
    for (k = ks; k <= i__1; ++k) {
	fk += 1.;
	wsave[k] = sin(fk * dt) * 2.;
/* L101: */
    }
    vrffti_(&np1, &wsave[kf + 1]);
    return 0;
} /* vsinti_ */

