/*
cp glassman-fft.f tmp.f
replace '^\*.*$' '' tmp.f
replace '!.*$' '' tmp.f
cat tmp.f | cpp | f2c > glassman-fft.c

the lib -lf2c doesn't seem to be necessary
*/

/*
  -- translated by f2c (version 19950110).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/*int main()  { return 0; }*/

/* Subroutine */ int spcfft_(u, n, isign, work, interp)
complex *u;
integer *n, *isign;
complex *work;
real *interp;
{
    /* System generated locals */
    integer i__1, i__2, i__3;
    real r__1;
    complex q__1;

    /* Local variables */
    static integer a, b, c, i;
    static logical1 scale;
    extern /* Subroutine */ int spcpft_();
    static logical1 inu;

    /* Parameter adjustments */
    --work;
    --u;

    /* Function Body */
    a = 1;
    b = *n;
    c = 1;
    inu = TRUE_;
    if (*isign == 1) {
	scale = TRUE_;
    } else if (*isign == -1) {
	scale = FALSE_;
    }
    while(b > 1) {
	a = c * a;
	c = 2;
	while(b % c != 0) {
	    ++c;
	}
	b /= c;
	if (inu) {
	    spcpft_(&a, &b, &c, &u[1], &work[1], isign);
	} else {
	    spcpft_(&a, &b, &c, &work[1], &u[1], isign);
	}
	inu = ! inu;
    }
    if (! inu) {
	i__1 = *n;
	for (i = 1; i <= i__1; ++i) {
	    i__2 = i;
	    i__3 = i;
	    u[i__2].r = work[i__3].r, u[i__2].i = work[i__3].i;
	}
    }
    if (scale) {
	i__1 = *n;
	for (i = 1; i <= i__1; ++i) {
	    i__2 = i;
	    i__3 = i;
	    r__1 = *n / *interp;
	    q__1.r = u[i__3].r / r__1, q__1.i = u[i__3].i / r__1;
	    u[i__2].r = q__1.r, u[i__2].i = q__1.i;
	}
    }
    return 0;
} /* spcfft_ */

/* Subroutine */ int spcpft_(a, b, c, uin, uout, isign)
integer *a, *b, *c;
complex *uin, *uout;
integer *isign;
{
    /* System generated locals */
    integer uin_dim1, uin_dim2, uin_offset, uout_dim1, uout_dim2, uout_offset,
	     i__1, i__2, i__3, i__4, i__5;
    doublereal d__1, d__2;
    complex q__1, q__2;
    doublecomplex z__1;

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

    /* Local variables */
    static doublereal angle;
    static complex omega, delta;
    static integer ia, ib, ic, jc, jcr;
    static complex sum;

    /* Parameter adjustments */
    uout_dim1 = *b;
    uout_dim2 = *a;
    uout_offset = uout_dim1 * (uout_dim2 + 1) + 1;
    uout -= uout_offset;
    uin_dim1 = *b;
    uin_dim2 = *c;
    uin_offset = uin_dim1 * (uin_dim2 + 1) + 1;
    uin -= uin_offset;

    /* Function Body */
    angle = (float)6.28318530717958 / (real) (*a * *c);
    omega.r = (float)1., omega.i = (float)0.;
    if (*isign == 1) {
	d__1 = cos(angle);
	d__2 = sin(angle);
	z__1.r = d__1, z__1.i = d__2;
	delta.r = z__1.r, delta.i = z__1.i;
    } else {
	d__1 = cos(angle);
	d__2 = -sin(angle);
	z__1.r = d__1, z__1.i = d__2;
	delta.r = z__1.r, delta.i = z__1.i;
    }
    i__1 = *c;
    for (ic = 1; ic <= i__1; ++ic) {
	i__2 = *a;
	for (ia = 1; ia <= i__2; ++ia) {
	    i__3 = *b;
	    for (ib = 1; ib <= i__3; ++ib) {
		i__4 = ib + (*c + ia * uin_dim2) * uin_dim1;
		sum.r = uin[i__4].r, sum.i = uin[i__4].i;
		i__4 = *c;
		for (jcr = 2; jcr <= i__4; ++jcr) {
		    jc = *c + 1 - jcr;
		    i__5 = ib + (jc + ia * uin_dim2) * uin_dim1;
		    q__2.r = omega.r * sum.r - omega.i * sum.i, q__2.i = 
			    omega.r * sum.i + omega.i * sum.r;
		    q__1.r = uin[i__5].r + q__2.r, q__1.i = uin[i__5].i + 
			    q__2.i;
		    sum.r = q__1.r, sum.i = q__1.i;
		}
		i__4 = ib + (ia + ic * uout_dim2) * uout_dim1;
		uout[i__4].r = sum.r, uout[i__4].i = sum.i;
	    }
	    q__1.r = delta.r * omega.r - delta.i * omega.i, q__1.i = delta.r *
		     omega.i + delta.i * omega.r;
	    omega.r = q__1.r, omega.i = q__1.i;
	}
    }
    return 0;
} /* spcpft_ */