/* harm.f -- translated by f2c (version of 23 April 1993  18:34:30).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "f2c.h"

/* Table of constant values */

static integer c__2 = 2;


/* Subroutine */ int harm_(real *a, integer *m, integer *inv, real *s, 
	integer *ifset, integer *iferr)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6;
    static integer equiv_2[3];

    /* Builtin functions */
    integer pow_ii(integer *, integer *);
    double sqrt(doublereal), asin(doublereal), sin(doublereal), cos(
	    doublereal);

    /* Local variables */
    real ak3_0__, ak3_1__;
    integer idif, jdif, kbit, ntsq, minn1, minn2, minn3;
    real root2;
    integer n1vnt, n2vnt, n3vnt, i, j, k, l;
#define n (equiv_2)
    integer ntvn3, ntvn2, ntvn1;
    real r, t;
    integer jjdif;
    real w[2], theta;
    integer ilast, jlast, klast, jstep, i2, k1, k2, m1;
#define n1 (equiv_2)
#define n2 (equiv_2 + 1)
#define n3 (equiv_2 + 2)
    integer m2, m3, k3, i3;
    real w2[2], w3[2];
    integer lm1exp, jstep2, id, jc, jd;
    real fn;
    integer il, mi, kl, jj, np[3], mt, nt, nx, jc1, lfirst, jj1, il1, jj3, 
	    jj2, mtlexp, ip1, ip2, ip3, jp3, jp2, jp1;
    real awi;
    integer mev, nnn;
    real awr;
    integer mtt, jjd1, jjd2, jjd3, igo1, igo2, igo3, ipp1, ipp2, ipp3, jpp3, 
	    jpp2, m1mt, m2mt, m3mt, jpp1, ntv2;
    real ak0_0__, ak0_1__, ak1_0__, ak1_1__, ak2_0__, ak2_1__;


/*        SUBROUTINE HARM */

/*        PURPOSE */
/*           PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX */
/*           THREE DIMENSIONAL ARRAY */

/*        USAGE */
/*           CALL HARM (A,M,INV,S,IFSET,IFERR) */

/*        NOTE */
/*           THESE COMMENT CARDS HAVE BEEN CORRECTED FOR SIGN ERROR */
/*           ORIGINALLY PRESENT IN  IFSET  THAT USED TO SPECIFY THE */
/*           INCORRECT DIRECTION OF THE FFT. */

/*        DESCRIPTION OF PARAMETERS */
/*           A     - AS INPUT, A CONTAINS THE COMPLEX, 3-DIMENSIONAL */
/*                   ARRAY TO BE TRANSFORMED.  THE REAL PART OF */
/*                   A(I1,I2,I3) IS STORED IN VECTOR FASHION IN A CELL */
/*                   WITH INDEX 2*(I3*N1*N2 + I2*N1 + I1) + 1 WHERE */
/*                   NI = 2**M(I), I=1,2,3 AND I1 = 0,1,...,N1-1 ETC. */
/*                   THE IMAGINARY PART IS IN THE CELL IMMEDIATELY */
/*                   FOLLOWING.  NOTE THAT THE SUBSCRIPT I1 INCREASES */
/*                   MOST RAPIDLY AND I3 INCREASES LEAST RAPIDLY. */
/*                   AS OUTPUT, A CONTAINS THE COMPLEX FOURIER */
/*                   TRANSFORM.  THE NUMBER OF CORE LOCATIONS OF */
/*                   ARRAY A IS 2*(N1*N2*N3) */
/*           M     - A THREE CELL VECTOR WHICH DETERMINES THE SIZES */
/*                   OF THE 3 DIMENSIONS OF THE ARRAY A.   THE SIZE, */
/*                   NI, OF THE I DIMENSION OF A IS 2**M(I), I = 1,2,3 */
/*           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION */
/*                   OF DIMENSION ONE FOURTH OF THE LARGEST OF N1,N2,N3 */
/*           S     - A VECTOR WORK AREA FOR SINE TABLES WITH DIMENSION */
/*                   THE SAME AS INV */
/*           IFSET - AN OPTION PARAMETER WITH THE FOLLOWING SETTINGS */
/*                      0    SET UP SINE AND INV TABLES ONLY */
/*                     -1    SET UP SINE AND INV TABLES ONLY AND */
/*                           CALCULATE FOURIER TRANSFORM */
/*                      1    SET UP SINE AND INV TABLES ONLY AND */
/*                           CALCULATE INVERSE FOURIER TRANSFORM (FOR */
/*                           THE MEANING OF INVERSE SEE THE EQUATIONS */
/*                           UNDER METHOD BELOW) */
/*                     -2    CALCULATE FOURIER TRANSFORM ONLY (ASSUME */
/*                           SINE AND INV TABLES EXIST) */
/*                      2    CALCULATE INVERSE FOURIER TRANSFORM ONLY */
/*                           (ASSUME SINE AND INV TABLES EXIST) */
/*           IFERR - ERROR INDICATOR.   WHEN IFSET IS 0,+1,-1, */
/*                   IFERR = 1 MEANS THE MAXIMUM M(I) IS GREATER THAN */
/*                  20 , I=1,2,3 */
/*                   " WHEN IFSET IS 2,-2 , IFERR = 1 */
/*                   MEANS THAT THE SINE AND INV TABLES ARE NOT LARGE */
/*                  ENOUGH OR HAVE NOT BEEN COMPUTED ." ..removed as not s
tandard*/
/*                   IF ON RETURN IFERR = 0 THEN NONE OF THE ABOVE */
/*                   CONDITIONS ARE PRESENT */

/*        REMARKS */
/*           THIS SUBROUTINE IS TO BE USED FOR COMPLEX, 3-DIMENSIONAL */
/*           ARRAYS IN WHICH EACH DIMENSION IS A POWER OF 2.  THE */
/*           MAXIMUM M(I) MUST NOT BE LESS THAN 3 OR GREATER THAN 20, */
/*           I = 1,2,3 */

/*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED */
/*           NONE */

/*        METHOD */
/*           FOR IFSET = -1, OR -2, THE FOURIER TRANSFORM OF COMPLEX */
/*           ARRAY A IS OBTAINED. */

/*                  N1-1   N2-1   N3-1                L1   L2   L3 */
/*     X(J1,J2,J3)=SUM    SUM    SUM    A(K1,K2,K3)*W1  *W2  *W3 */
/*                  K1=0   K2=0   K3=0 */

/*                  WHERE WI IS THE N(I) ROOT OF UNITY AND L1=K1*J1, */
/*                        L2=K2*J2, L3=K3*J3 */


/*           FOR IFSET = +1, OR +2, THE INVERSE FOURIER TRANSFORM A OF */
/*           COMPLEX ARRAY X IS OBTAINED. */

/*     A(K1,K2,K3)= */
/*               1      N1-1   N2-1   N3-1                -L1  -L2  -L3 */
/*           -------- *SUM    SUM    SUM    X(J1,J2,J3)*W1  *W2  *W3 */
/*           N1*N2*N3   J1=0   J2=0   J3=0 */


/*           SEE J.W. COOLEY AND J.W. TUKEY, @AN ALGORITHM FOR THE */
/*           MACHINE CALCULATION OF COMPLEX FOURIER SERIES@, */
/*           MATHEMATICS OF COMPUTATIONS, VOL. 19 (APR. 1965), P. 297. */

/************************************************************************
*******/
/*                   THIS ROUTINE IS PUBLIC DOMAIN                        
   **/
/************************************************************************
*******/


/* * */
/* *    The following are parameters which may be overwritten. */
/* *    The original coding in fact assumed that these variables */
/* *    would be retained */
    /* Parameter adjustments */
    --s;
    --inv;
    --m;
    --a;

    /* Function Body */
/* L869: */
/* L879: */
/* Computing MAX */
    i__1 = max(m[1],m[2]);
    mt = max(i__1,m[3]) - 2;
    mt = max(2,mt);
    nt = pow_ii(&c__2, &mt);
/* * */
/* L10: */
    if (abs(*ifset) - 1 <= 0) {
	goto L900;
    } else {
	goto L12;
    }
L12:
/* Computing MAX */
    i__1 = max(m[1],m[2]);
    mtt = max(i__1,m[3]) - 2;
    root2 = sqrt((float)2.);
    if (mtt - mt <= 0) {
	goto L14;
    } else {
	goto L13;
    }
L13:
    *iferr = 1;
    return 0;
L14:
    *iferr = 0;
    m1 = m[1];
    m2 = m[2];
    m3 = m[3];
    *n1 = pow_ii(&c__2, &m1);
    *n2 = pow_ii(&c__2, &m2);
    *n3 = pow_ii(&c__2, &m3);
/* L16: */
    if (*ifset >= 0) {
	goto L20;
    } else {
	goto L18;
    }
L18:
    nx = *n1 * *n2 * *n3;
    fn = (real) nx;
/*     / may be faster than * */
    nnn = nx << 1;
    i__1 = nnn;
    for (i = 1; i <= i__1; i += 2) {
	a[i] /= fn;
/* L19: */
	a[i + 1] = -(doublereal)a[i + 1] / fn;
    }
L20:
    np[0] = *n1 << 1;
    np[1] = np[0] * *n2;
    np[2] = np[1] * *n3;
    for (id = 1; id <= 3; ++id) {
	il = np[2] - np[id - 1];
	il1 = il + 1;
	mi = m[id];
	if (mi <= 0) {
	    goto L250;
	} else {
	    goto L30;
	}
L30:
	idif = np[id - 1];
	kbit = np[id - 1];
	mev = mi / 2 << 1;
	if (mi - mev <= 0) {
	    goto L60;
	} else {
	    goto L40;
	}

/*     M IS ODD. DO L=1 CASE */
L40:
	kbit /= 2;
	kl = kbit - 2;
	i__1 = il1;
	i__2 = idif;
	for (i = 1; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
	    klast = kl + i;
	    i__3 = klast;
	    for (k = i; k <= i__3; k += 2) {
		k1 = k + kbit;

/*     DO ONE STEP WITH L=1,J=0 */
/*     A(K)=A(K)+A(K1) */
/*     A(K1)=A(K)-A(K1) */

/*     NOTE: USE AK0_0 AS A little unwrapping */
/*      T=A(K1) */
/*      AK0_0=A(K) */
/*      A(K1)=AK0_0-T */
/*      A(K)=AK0_0+T */

/*      T=A(K1+1) */
/*      AK0_0=A(K+1) */
/*      A(K1+1)=AK0_0-T */
/*   50 A(K+1)=AK0_0+T */
		ak0_0__ = a[k];
		ak0_1__ = a[k + 1];
		ak1_0__ = a[k1];
		ak1_1__ = a[k1 + 1];
		a[k] = ak0_0__ + ak1_0__;
		a[k + 1] = ak0_1__ + ak1_1__;
		a[k1] = ak0_0__ - ak1_0__;
		a[k1 + 1] = ak0_1__ - ak1_1__;
/* L50: */
	    }
	}
	if (mi - 1 <= 0) {
	    goto L250;
	} else {
	    goto L52;
	}
L52:
	lfirst = 3;

/*     DEF - JLAST = 2**(L-2) -1 */
	jlast = 1;
	goto L70;

/*     M IS EVEN */
L60:
	lfirst = 2;
	jlast = 0;
L70:
	i__3 = mi;
	for (l = lfirst; l <= i__3; l += 2) {
	    jjdif = kbit;
	    kbit /= 4;
	    kl = kbit - 2;

/*     DO FOR J=0 */
/*      CALL CLOCK(TIME1) */
	    i__2 = il1;
	    i__1 = idif;
	    for (i = 1; i__1 < 0 ? i >= i__2 : i <= i__2; i += i__1) {
		klast = i + kl;
		i__4 = klast;
		for (k = i; k <= i__4; k += 2) {
		    k1 = k + kbit;
		    k2 = k1 + kbit;
		    k3 = k2 + kbit;

/*     DO TWO STEPS WITH J=0 */
/*     A(K)=A(K)+A(K2) */
/*     A(K2)=A(K)-A(K2) */
/*     A(K1)=A(K1)+A(K3) */
/*     A(K3)=A(K1)-A(K3) */

/*     A(K)=A(K)+A(K1) */
/*     A(K1)=A(K)-A(K1) */
/*     A(K2)=A(K2)+A(K3)*I */
/*     A(K3)=A(K2)-A(K3)*I */


/*     unwind variables */
		    ak0_0__ = a[k];
		    ak0_1__ = a[k + 1];
		    ak1_0__ = a[k1];
		    ak1_1__ = a[k1 + 1];
		    ak2_0__ = a[k2];
		    ak2_1__ = a[k2 + 1];
		    ak3_0__ = a[k3];
		    ak3_1__ = a[k3 + 1];

		    t = ak2_0__;
		    ak2_0__ = ak0_0__ - t;
		    ak0_0__ += t;
		    t = ak2_1__;
		    ak2_1__ = ak0_1__ - t;
		    ak0_1__ += t;

		    t = ak3_0__;
		    ak3_0__ = ak1_0__ - t;
		    ak1_0__ += t;
		    t = ak3_1__;
		    ak3_1__ = ak1_1__ - t;
		    ak1_1__ += t;

/*     small change in the order here */
/*      T=AK1_0 */
/*      R=AK0_0 */
/*      T=AK1_1 */
/*      R=AK0_1 */
		    a[k] = ak0_0__ + ak1_0__;
		    a[k + 1] = ak0_1__ + ak1_1__;
		    a[k1] = ak0_0__ - ak1_0__;
		    a[k1 + 1] = ak0_1__ - ak1_1__;

/*      R=-AK3_1 */
/*      T = AK3_0 */
/*      T2=AK2_0 */
/*      T3=AK2_1 */
		    a[k2] = ak2_0__ - ak3_1__;
		    a[k2 + 1] = ak2_1__ + ak3_0__;
		    a[k3] = ak2_0__ + ak3_1__;
		    a[k3 + 1] = ak2_1__ - ak3_0__;
/* L80: */
		}
	    }
	    if (jlast <= 0) {
		goto L235;
	    } else {
		goto L82;
	    }
L82:
	    jj = jjdif + 1;

/*     DO FOR J=1 */
	    ilast = il + jj;
	    i__4 = ilast;
	    i__1 = idif;
	    for (i = jj; i__1 < 0 ? i >= i__4 : i <= i__4; i += i__1) {
		klast = kl + i;
		i__2 = klast;
		for (k = i; k <= i__2; k += 2) {
		    k1 = k + kbit;
		    k2 = k1 + kbit;
		    k3 = k2 + kbit;

/*     LETTING W=(1+I)/ROOT2,W3=(-1+I)/ROOT2,W2=I, */
/*     A(K)=A(K)+A(K2)*I */
/*     A(K2)=A(K)-A(K2)*I */
/*     A(K1)=A(K1)*W+A(K3)*W3 */
/*     A(K3)=A(K1)*W-A(K3)*W3 */

/*     A(K)=A(K)+A(K1) */
/*     A(K1)=A(K)-A(K1) */
/*     A(K2)=A(K2)+A(K3)*I */
/*     A(K3)=A(K2)-A(K3)*I */


/*     unwind variables */
		    ak0_0__ = a[k];
		    ak0_1__ = a[k + 1];
		    ak1_0__ = a[k1];
		    ak1_1__ = a[k1 + 1];
		    ak2_0__ = a[k2];
		    ak2_1__ = a[k2 + 1];
		    ak3_0__ = a[k3];
		    ak3_1__ = a[k3 + 1];
		    r = -(doublereal)ak2_1__;
		    t = ak2_0__;
		    ak2_0__ = ak0_0__ - r;
		    ak0_0__ += r;
		    ak2_1__ = ak0_1__ - t;
		    ak0_1__ += t;

		    awr = ak1_0__ - ak1_1__;
		    awi = ak1_1__ + ak1_0__;
		    r = -(doublereal)ak3_0__ - ak3_1__;
		    t = ak3_0__ - ak3_1__;
		    ak3_0__ = (awr - r) / root2;
		    ak3_1__ = (awi - t) / root2;
		    ak1_0__ = (awr + r) / root2;
		    ak1_1__ = (awi + t) / root2;
/*     order change here */
/*      T= AK1_0 */
/*      T=AK1_1 */
		    a[k] = ak0_0__ + ak1_0__;
		    a[k + 1] = ak0_1__ + ak1_1__;
		    a[k1] = ak0_0__ - ak1_0__;
		    a[k1 + 1] = ak0_1__ - ak1_1__;

/*      R=-AK3_1 */
/*      T=AK3_0 */
		    a[k2] = ak2_0__ - ak3_1__;
		    a[k2 + 1] = ak2_1__ + ak3_0__;
		    a[k3] = ak2_0__ + ak3_1__;
		    a[k3 + 1] = ak2_1__ - ak3_0__;
/* L85: */
		}
	    }
	    if (jlast - 1 <= 0) {
		goto L235;
	    } else {
		goto L90;
	    }
L90:
	    jj += jjdif;

/*     NOW DO THE REMAINING J@S */
	    i__2 = jlast;
	    for (j = 2; j <= i__2; ++j) {
/*     THIS EQUIVALENCE FORCES SOME VARIABLES INTO PERMANENT R
EGISTER */
/*      EQUIVALENCE (I2,K,I3CC),(K1,I2C,I3),(K2,I2CC,I3C) */

/*     FETCH W@S */
/*     DEF- W=W**INV(J), W2=W**2, W3=W**3 */
/* L96: */
		i = inv[j + 1];
/*   98 IC=NT-I */
		w[0] = s[nt - i];
		w[1] = s[i];
		k = i + i;
		k1 = nt - k;
		if (k1 < 0) {
		    goto L120;
		} else if (k1 == 0) {
		    goto L110;
		} else {
		    goto L100;
		}
/*     2*I IS IN FIRST QUADRANT */
L100:
		w2[0] = s[k1];
		w2[1] = s[k];
		goto L130;
L110:
		w2[0] = (float)0.;
		w2[1] = (float)1.;
		goto L130;

/*     2*I IS IN SECOND QUADRANT */
L120:
		k2 = k1 + nt;
		k1 = -k1;
		w2[0] = -(doublereal)s[k1];
		w2[1] = s[k2];
L130:
		k1 = i + k;
		k2 = nt - k1;
		if (k2 < 0) {
		    goto L160;
		} else if (k2 == 0) {
		    goto L150;
		} else {
		    goto L140;
		}

/*     I3 IN FIRST QUADRANT */
L140:
		w3[0] = s[k2];
		w3[1] = s[k1];
		goto L200;
L150:
		w3[0] = (float)0.;
		w3[1] = (float)1.;
		goto L200;

L160:
		k = k2 + nt;
		if (k < 0) {
		    goto L190;
		} else if (k == 0) {
		    goto L180;
		} else {
		    goto L170;
		}

/*     I3 IN SECOND QUADRANT */
L170:
		k2 = -k2;
		w3[0] = -(doublereal)s[k2];
		w3[1] = s[k];
		goto L200;
L180:
		w3[0] = (float)-1.;
		w3[1] = (float)0.;
		goto L200;

/*     3*I IN THIRD QUADRANT */
L190:
		k1 = nt + k;
		k = -k;
		w3[0] = -(doublereal)s[k1];
		w3[1] = -(doublereal)s[k];
L200:
		ilast = il + jj;
		i__1 = ilast;
		i__4 = idif;
		for (i = jj; i__4 < 0 ? i >= i__1 : i <= i__1; i += i__4) {
		    klast = kl + i;
		    i__5 = klast;
		    for (k = i; k <= i__5; k += 2) {
			k1 = k + kbit;
			k2 = k1 + kbit;
/*      K3=K2+KBIT */

/*     DO TWO STEPS WITH J NOT 0 */
/*     A(K)=A(K)+A(K2)*W2 */
/*     A(K2)=A(K)-A(K2)*W2 */
/*     A(K1)=A(K1)*W+A(K3)*W3 */
/*     A(K3)=A(K1)*W-A(K3)*W3 */

/*     A(K)=A(K)+A(K1) */
/*     A(K1)=A(K)-A(K1) */
/*     A(K2)=A(K2)+A(K3)*I */
/*     A(K3)=A(K2)-A(K3)*I */

/*     unwind variables */
			ak0_0__ = a[k];
			ak0_1__ = a[k + 1];
			ak1_0__ = a[k1];
			ak1_1__ = a[k1 + 1];
			ak2_0__ = a[k2];
			ak2_1__ = a[k2 + 1];
			ak3_0__ = a[k2 + kbit];
			ak3_1__ = a[k2 + kbit + 1];
			r = ak2_0__ * w2[0] - ak2_1__ * w2[1];
			t = ak2_0__ * w2[1] + ak2_1__ * w2[0];
			ak2_0__ = ak0_0__ - r;
			ak0_0__ += r;
			ak2_1__ = ak0_1__ - t;
			ak0_1__ += t;

			r = ak3_0__ * w3[0] - ak3_1__ * w3[1];
			t = ak3_0__ * w3[1] + ak3_1__ * w3[0];
			awr = ak1_0__ * w[0] - ak1_1__ * w[1];
			awi = ak1_0__ * w[1] + ak1_1__ * w[0];
			ak3_0__ = awr - r;
			ak3_1__ = awi - t;
			ak1_0__ = awr + r;
			ak1_1__ = awi + t;

/*     order change */
/*     T=ak1_0 */
			a[k] = ak0_0__ + ak1_0__;
			a[k + 1] = ak0_1__ + ak1_1__;
			a[k1] = ak0_0__ - ak1_0__;
/*     T=ak1_1 */
			a[k1 + 1] = ak0_1__ - ak1_1__;
/*     R=-ak3_1 */
/*     T=ak3_0 */
			a[k2] = ak2_0__ - ak3_1__;
			a[k2 + 1] = ak2_1__ + ak3_0__;
			a[k2 + kbit] = ak2_0__ + ak3_1__;
			a[k2 + kbit + 1] = ak2_1__ - ak3_0__;
/* L220: */
		    }
		}
/*     END OF I AND K LOOPS */

/* L230: */
		jj = jjdif + jj;
	    }
/*     END OF J-LOOP */

L235:
	    jlast = (jlast << 2) + 3;
/* L240: */
	}
/*     END OF  L  LOOP */
L250:
	;
    }
/*     END OF  ID  LOOP */

/*     WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE */
/*     BIT-REVERSED.  THE FOLLOWING ROUTINE PUTS THEM IN ORDER */
    ntsq = nt * nt;
    m3mt = m3 - mt;
/* L350: */
    if (m3mt >= 0) {
	goto L360;
    } else {
	goto L370;
    }

/*     M3 GR. OR EQ. MT */
L360:
    igo3 = 1;
    n3vnt = *n3 / nt;
    minn3 = nt;
    goto L380;

/*     M3 LESS THAN MT */
L370:
    igo3 = 2;
    n3vnt = 1;
    ntvn3 = nt / *n3;
    minn3 = *n3;
L380:
    jjd3 = ntsq / *n3;
    m2mt = m2 - mt;
/* L450: */
    if (m2mt >= 0) {
	goto L460;
    } else {
	goto L470;
    }

/*     M2 GR. OR EQ. MT */
L460:
    igo2 = 1;
    n2vnt = *n2 / nt;
    minn2 = nt;
    goto L480;

/*     M2 LESS THAN MT */
L470:
    igo2 = 2;
    n2vnt = 1;
    ntvn2 = nt / *n2;
    minn2 = *n2;
L480:
    jjd2 = ntsq / *n2;
    m1mt = m1 - mt;
/* L550: */
    if (m1mt >= 0) {
	goto L560;
    } else {
	goto L570;
    }

/*     M1 GR. OR EQ. MT */
L560:
    igo1 = 1;
    n1vnt = *n1 / nt;
    minn1 = nt;
    goto L580;

/*     M1 LESS THAN MT */
L570:
    igo1 = 2;
    n1vnt = 1;
    ntvn1 = nt / *n1;
    minn1 = *n1;
L580:
    jjd1 = ntsq / *n1;
/* L600: */
    jj3 = 1;
    j = 1;
    i__3 = n3vnt;
    for (jpp3 = 1; jpp3 <= i__3; ++jpp3) {
	ipp3 = inv[jj3];
	i__2 = minn3;
	for (jp3 = 1; jp3 <= i__2; ++jp3) {
	    switch (igo3) {
		case 1:  goto L610;
		case 2:  goto L620;
	    }
L610:
	    ip3 = inv[jp3] * n3vnt;
	    goto L630;
L620:
	    ip3 = inv[jp3] / ntvn3;
L630:
	    i3 = (ipp3 + ip3) * *n2;
/* L700: */
	    jj2 = 1;
	    i__5 = n2vnt;
	    for (jpp2 = 1; jpp2 <= i__5; ++jpp2) {
		ipp2 = inv[jj2] + i3;
		i__4 = minn2;
		for (jp2 = 1; jp2 <= i__4; ++jp2) {
		    switch (igo2) {
			case 1:  goto L710;
			case 2:  goto L720;
		    }
L710:
		    ip2 = inv[jp2] * n2vnt;
		    goto L730;
L720:
		    ip2 = inv[jp2] / ntvn2;
L730:
		    i2 = (ipp2 + ip2) * *n1;
/* L800: */
		    jj1 = 1;
		    i__1 = n1vnt;
		    for (jpp1 = 1; jpp1 <= i__1; ++jpp1) {
			ipp1 = inv[jj1] + i2;
			i__6 = minn1;
			for (jp1 = 1; jp1 <= i__6; ++jp1) {
			    switch (igo1) {
				case 1:  goto L810;
				case 2:  goto L820;
			    }
L810:
			    ip1 = inv[jp1] * n1vnt;
			    goto L830;
L820:
			    ip1 = inv[jp1] / ntvn1;
L830:
			    i = (ipp1 + ip1 << 1) + 1;
			    if (j - i >= 0) {
				goto L850;
			    } else {
				goto L840;
			    }
L840:
			    t = a[i];
			    a[i] = a[j];
			    a[j] = t;
			    t = a[i + 1];
			    a[i + 1] = a[j + 1];
			    a[j + 1] = t;
L850:
			    j += 2;
			}
/* L860: */
			jj1 += jjd1;
		    }
		}
/*     END OF JPP1 AND JP2 */

/* L870: */
		jj2 += jjd2;
	    }
	}
/*     END OF JPP2 AND JP3 LOOPS */

/* L880: */
	jj3 += jjd3;
    }
/*     END OF JPP3 LOOP */

/* L890: */
    if (*ifset >= 0) {
	goto L895;
    } else {
	goto L891;
    }
L891:
    nnn = nx << 1;
/*      fn=1./nx */
/*      DO 892 I = 1,Nnn,2 */
/*      a(i)=a(i)*fn */
/*  892 A(i+1) = -A(i+1)*fn */
    i__3 = nnn;
    for (i = 2; i <= i__3; i += 2) {
/* L892: */
	a[i] = -(doublereal)a[i];
    }
L895:
    return 0;

/*     THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES. */

L900:
/* Computing MAX */
    i__3 = max(m[1],m[2]);
    mt = max(i__3,m[3]) - 2;
    mt = max(2,mt);
/* L904: */
    if (mt - 18 <= 0) {
	goto L906;
    } else {
	goto L13;
    }
L906:
    *iferr = 0;
    nt = pow_ii(&c__2, &mt);
    ntv2 = nt / 2;

/*     SET UP SIN TABLE */
/*     THETA=PIE/2**(L+1) FOR L=1 */
/*  910 THETA=.7853981634 */
/* L910: */
    theta = asin((float)1.) / (float)2.;

/*     JSTEP=2**(MT-L+1) FOR L=1 */
    jstep = nt;

/*     JDIF=2**(MT-L) FOR L=1 */
    jdif = ntv2;
    s[jdif] = sin(theta);
    i__3 = mt;
    for (l = 2; l <= i__3; ++l) {
	theta /= (float)2.;
	jstep2 = jstep;
	jstep = jdif;
	jdif = jstep / 2;
	s[jdif] = sin(theta);
	jc1 = nt - jdif;
	s[jc1] = cos(theta);
	jlast = nt - jstep2;
	if (jlast - jstep >= 0) {
	    goto L920;
	} else {
	    goto L950;
	}
L920:
	i__5 = jlast;
	i__2 = jstep;
	for (j = jstep; i__2 < 0 ? j >= i__5 : j <= i__5; j += i__2) {
	    jc = nt - j;
	    jd = j + jdif;
/* L940: */
	    s[jd] = s[j] * s[jc1] + s[jdif] * s[jc];
	}
L950:
	;
    }

/*     SET UP INV(J) TABLE */

/* L960: */
    mtlexp = ntv2;

/*     MTLEXP=2**(MT-L). FOR L=1 */
    lm1exp = 1;

/*     LM1EXP=2**(L-1). FOR L=1 */
    inv[1] = 0;
    i__3 = mt;
    for (l = 1; l <= i__3; ++l) {
	inv[lm1exp + 1] = mtlexp;
	i__2 = lm1exp;
	for (j = 2; j <= i__2; ++j) {
	    jj = j + lm1exp;
/* L970: */
	    inv[jj] = inv[j] + mtlexp;
	}
	mtlexp /= 2;
/* L980: */
	lm1exp <<= 1;
    }
/* L982: */
    if (*ifset != 0) {
	goto L12;
    } else {
	goto L895;
    }
} /* harm_ */

#undef n3
#undef n2
#undef n1
#undef n


#ifdef __cplusplus
	}
#endif