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

#include "f2c.h"

/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/* Subroutine */ int fourn_(real *data, integer *nn, integer *ndim, integer *
	isign)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6;
    doublereal d__1;

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

    /* Local variables */
    static integer n, i1, i2, i3, k1, k2;
    static doublereal wi, wr;
    static integer ip1, ip2, ip3;
    static doublereal wpi, wpr;
    static integer ifp1, ifp2, idim, ibit, nrem, ntot, i2rev, i3rev;
    static doublereal theta;
    static real tempi, tempr;
    static integer nprev;
    static doublereal wtemp;

/*       do 111 i=1,10 */
/*       write(*,*) data(i) */
/* 111    CONTINUE */
/*        do 222 i=1,10 */
/*        write(*,*) nn(i) */
/* 222    CONTINUE */
/*       write(*,*) NDIM,ISIGN */
    /* Parameter adjustments */
    --data;
    --nn;

    /* Function Body */
    ntot = 1;
    i__1 = *ndim;
    for (idim = 1; idim <= i__1; ++idim) {
	ntot *= nn[idim];
/* L11: */
    }
    nprev = 1;
    i__1 = *ndim;
    for (idim = 1; idim <= i__1; ++idim) {
	n = nn[idim];
	nrem = ntot / (n * nprev);
	ip1 = nprev << 1;
	ip2 = ip1 * n;
	ip3 = ip2 * nrem;
	i2rev = 1;
	i__2 = ip2;
	i__3 = ip1;
	for (i2 = 1; i__3 < 0 ? i2 >= i__2 : i2 <= i__2; i2 += i__3) {
	    if (i2 < i2rev) {
		i__4 = i2 + ip1 - 2;
		for (i1 = i2; i1 <= i__4; i1 += 2) {
		    i__5 = ip3;
		    i__6 = ip2;
		    for (i3 = i1; i__6 < 0 ? i3 >= i__5 : i3 <= i__5; i3 += 
			    i__6) {
			i3rev = i2rev + i3 - i2;
			tempr = data[i3];
			tempi = data[i3 + 1];
			data[i3] = data[i3rev];
			data[i3 + 1] = data[i3rev + 1];
			data[i3rev] = tempr;
			data[i3rev + 1] = tempi;
/* L12: */
		    }
/* L13: */
		}
	    }
	    ibit = ip2 / 2;
L1:
	    if (ibit >= ip1 && i2rev > ibit) {
		i2rev -= ibit;
		ibit /= 2;
		goto L1;
	    }
	    i2rev += ibit;
/* L14: */
	}
	ifp1 = ip1;
L2:
	if (ifp1 < ip2) {
	    ifp2 = ifp1 << 1;
	    theta = *isign * 6.28318530717959 / (ifp2 / ip1);
/* Computing 2nd power */
	    d__1 = sin(theta * .5);
	    wpr = d__1 * d__1 * -2.;
	    wpi = sin(theta);
	    wr = 1.;
	    wi = 0.;
	    i__3 = ifp1;
	    i__2 = ip1;
	    for (i3 = 1; i__2 < 0 ? i3 >= i__3 : i3 <= i__3; i3 += i__2) {
		i__4 = i3 + ip1 - 2;
		for (i1 = i3; i1 <= i__4; i1 += 2) {
		    i__6 = ip3;
		    i__5 = ifp2;
		    for (i2 = i1; i__5 < 0 ? i2 >= i__6 : i2 <= i__6; i2 += 
			    i__5) {
			k1 = i2;
			k2 = k1 + ifp1;
			tempr = (real) wr * data[k2] - (real) wi * data[k2 + 
				1];
			tempi = (real) wr * data[k2 + 1] + (real) wi * data[
				k2];
			data[k2] = data[k1] - tempr;
			data[k2 + 1] = data[k1 + 1] - tempi;
			data[k1] += tempr;
			data[k1 + 1] += tempi;
/* L15: */
		    }
/* L16: */
		}
		wtemp = wr;
		wr = wr * wpr - wi * wpi + wr;
		wi = wi * wpr + wtemp * wpi + wi;
/* L17: */
	    }
	    ifp1 = ifp2;
	    goto L2;
	}
	nprev = n * nprev;
/* L18: */
    }
    return 0;
} /* fourn_ */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/* Subroutine */ int svbksb_(real *u, real *w, real *v, integer *m, integer *
	n, integer *mp, integer *np, real *b, real *x)
{
    /* System generated locals */
    integer u_dim1, u_offset, v_dim1, v_offset, i__1, i__2;

    /* Local variables */
    static integer i__, j;
    static real s;
    static integer jj;
    static real tmp[100];

/*       write(*,*) 'u(1,1)=',u(1,1),'w(1)=',w(1),'v(1,1)=',v(1,1) */
/*       write(*,*) 'm=',m,'n=',n,'mp=',mp,'np=',np */
/*       write(*,*) 'b(1)=',b(1),'x(1)=',x(1) */
    /* Parameter adjustments */
    --b;
    --x;
    v_dim1 = *np;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --w;
    u_dim1 = *mp;
    u_offset = 1 + u_dim1;
    u -= u_offset;

    /* Function Body */
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	s = 0.f;
	if (w[j] != 0.f) {
	    i__2 = *m;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		s += u[i__ + j * u_dim1] * b[i__];
/* L11: */
	    }
	    s /= w[j];
	}
	tmp[j - 1] = s;
/* L12: */
    }
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	s = 0.f;
	i__2 = *n;
	for (jj = 1; jj <= i__2; ++jj) {
	    s += v[j + jj * v_dim1] * tmp[jj - 1];
/* L13: */
	}
	x[j] = s;
/* L14: */
    }
    return 0;
} /* svbksb_ */

/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/* Subroutine */ int svdcmp_(real *a, integer *m, integer *n, integer *mp, 
	integer *np, real *w, real *v)
{
    /* System generated locals */
    integer a_dim1, a_offset, v_dim1, v_offset, i__1, i__2, i__3;
    real r__1, r__2, r__3, r__4;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    static real c__, f, g, h__;
    static integer i__, j, k, l;
    static real s, x, y, z__;
    static integer nm;
    static real rv1[200];
    static integer its;
    static real scale, anorm;

/*       WRITE(*,*) 'CALLING SVDCMP....',M,N,MP,NP */
/*       do 111 i=1,M */
/*       write (*,*) 'A<1,x>',a(i,1) */
/* 111    continue */
/*       do 222 i=1,5 */
/*       write(*,*) 'W<x>:',w(i) */
/* 222    continue */
/*       do 333 i=1,5 */
/*       write(*,*) 'V<1,x>',v(i,1) */
/* 333    continue */
    /* Parameter adjustments */
    v_dim1 = *np;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --w;
    a_dim1 = *mp;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    g = 0.f;
    scale = 0.f;
    anorm = 0.f;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	l = i__ + 1;
	rv1[i__ - 1] = scale * g;
	g = 0.f;
	s = 0.f;
	scale = 0.f;
	if (i__ <= *m) {
	    i__2 = *m;
	    for (k = i__; k <= i__2; ++k) {
		scale += (r__1 = a[k + i__ * a_dim1], dabs(r__1));
/* L11: */
	    }
	    if (scale != 0.f) {
		i__2 = *m;
		for (k = i__; k <= i__2; ++k) {
		    a[k + i__ * a_dim1] /= scale;
		    s += a[k + i__ * a_dim1] * a[k + i__ * a_dim1];
/* L12: */
		}
		f = a[i__ + i__ * a_dim1];
		r__1 = sqrt(s);
		g = -r_sign(&r__1, &f);
		h__ = f * g - s;
		a[i__ + i__ * a_dim1] = f - g;
		if (i__ != *n) {
		    i__2 = *n;
		    for (j = l; j <= i__2; ++j) {
			s = 0.f;
			i__3 = *m;
			for (k = i__; k <= i__3; ++k) {
			    s += a[k + i__ * a_dim1] * a[k + j * a_dim1];
/* L13: */
			}
			f = s / h__;
			i__3 = *m;
			for (k = i__; k <= i__3; ++k) {
			    a[k + j * a_dim1] += f * a[k + i__ * a_dim1];
/* L14: */
			}
/* L15: */
		    }
		}
		i__2 = *m;
		for (k = i__; k <= i__2; ++k) {
		    a[k + i__ * a_dim1] = scale * a[k + i__ * a_dim1];
/* L16: */
		}
	    }
	}
	w[i__] = scale * g;
	g = 0.f;
	s = 0.f;
	scale = 0.f;
	if (i__ <= *m && i__ != *n) {
	    i__2 = *n;
	    for (k = l; k <= i__2; ++k) {
		scale += (r__1 = a[i__ + k * a_dim1], dabs(r__1));
/* L17: */
	    }
	    if (scale != 0.f) {
		i__2 = *n;
		for (k = l; k <= i__2; ++k) {
		    a[i__ + k * a_dim1] /= scale;
		    s += a[i__ + k * a_dim1] * a[i__ + k * a_dim1];
/* L18: */
		}
		f = a[i__ + l * a_dim1];
		r__1 = sqrt(s);
		g = -r_sign(&r__1, &f);
		h__ = f * g - s;
		a[i__ + l * a_dim1] = f - g;
		i__2 = *n;
		for (k = l; k <= i__2; ++k) {
		    rv1[k - 1] = a[i__ + k * a_dim1] / h__;
/* L19: */
		}
		if (i__ != *m) {
		    i__2 = *m;
		    for (j = l; j <= i__2; ++j) {
			s = 0.f;
			i__3 = *n;
			for (k = l; k <= i__3; ++k) {
			    s += a[j + k * a_dim1] * a[i__ + k * a_dim1];
/* L21: */
			}
			i__3 = *n;
			for (k = l; k <= i__3; ++k) {
			    a[j + k * a_dim1] += s * rv1[k - 1];
/* L22: */
			}
/* L23: */
		    }
		}
		i__2 = *n;
		for (k = l; k <= i__2; ++k) {
		    a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1];
/* L24: */
		}
	    }
	}
/* Computing MAX */
	r__3 = anorm, r__4 = (r__1 = w[i__], dabs(r__1)) + (r__2 = rv1[i__ - 
		1], dabs(r__2));
	anorm = dmax(r__3,r__4);
/* L25: */
    }
    for (i__ = *n; i__ >= 1; --i__) {
	if (i__ < *n) {
	    if (g != 0.f) {
		i__1 = *n;
		for (j = l; j <= i__1; ++j) {
		    v[j + i__ * v_dim1] = a[i__ + j * a_dim1] / a[i__ + l * 
			    a_dim1] / g;
/* L26: */
		}
		i__1 = *n;
		for (j = l; j <= i__1; ++j) {
		    s = 0.f;
		    i__2 = *n;
		    for (k = l; k <= i__2; ++k) {
			s += a[i__ + k * a_dim1] * v[k + j * v_dim1];
/* L27: */
		    }
		    i__2 = *n;
		    for (k = l; k <= i__2; ++k) {
			v[k + j * v_dim1] += s * v[k + i__ * v_dim1];
/* L28: */
		    }
/* L29: */
		}
	    }
	    i__1 = *n;
	    for (j = l; j <= i__1; ++j) {
		v[i__ + j * v_dim1] = 0.f;
		v[j + i__ * v_dim1] = 0.f;
/* L31: */
	    }
	}
	v[i__ + i__ * v_dim1] = 1.f;
	g = rv1[i__ - 1];
	l = i__;
/* L32: */
    }
    for (i__ = *n; i__ >= 1; --i__) {
	l = i__ + 1;
	g = w[i__];
	if (i__ < *n) {
	    i__1 = *n;
	    for (j = l; j <= i__1; ++j) {
		a[i__ + j * a_dim1] = 0.f;
/* L33: */
	    }
	}
	if (g != 0.f) {
	    g = 1.f / g;
	    if (i__ != *n) {
		i__1 = *n;
		for (j = l; j <= i__1; ++j) {
		    s = 0.f;
		    i__2 = *m;
		    for (k = l; k <= i__2; ++k) {
			s += a[k + i__ * a_dim1] * a[k + j * a_dim1];
/* L34: */
		    }
		    f = s / a[i__ + i__ * a_dim1] * g;
		    i__2 = *m;
		    for (k = i__; k <= i__2; ++k) {
			a[k + j * a_dim1] += f * a[k + i__ * a_dim1];
/* L35: */
		    }
/* L36: */
		}
	    }
	    i__1 = *m;
	    for (j = i__; j <= i__1; ++j) {
		a[j + i__ * a_dim1] *= g;
/* L37: */
	    }
	} else {
	    i__1 = *m;
	    for (j = i__; j <= i__1; ++j) {
		a[j + i__ * a_dim1] = 0.f;
/* L38: */
	    }
	}
	a[i__ + i__ * a_dim1] += 1.f;
/* L39: */
    }
    for (k = *n; k >= 1; --k) {
	for (its = 1; its <= 30; ++its) {
	    for (l = k; l >= 1; --l) {
		nm = l - 1;
		if ((r__1 = rv1[l - 1], dabs(r__1)) + anorm == anorm) {
		    goto L2;
		}
		if ((r__1 = w[nm], dabs(r__1)) + anorm == anorm) {
		    goto L1;
		}
/* L41: */
	    }
L1:
	    c__ = 0.f;
	    s = 1.f;
	    i__1 = k;
	    for (i__ = l; i__ <= i__1; ++i__) {
		f = s * rv1[i__ - 1];
		if (dabs(f) + anorm != anorm) {
		    g = w[i__];
		    h__ = sqrt(f * f + g * g);
		    w[i__] = h__;
		    h__ = 1.f / h__;
		    c__ = g * h__;
		    s = -(f * h__);
		    i__2 = *m;
		    for (j = 1; j <= i__2; ++j) {
			y = a[j + nm * a_dim1];
			z__ = a[j + i__ * a_dim1];
			a[j + nm * a_dim1] = y * c__ + z__ * s;
			a[j + i__ * a_dim1] = -(y * s) + z__ * c__;
/* L42: */
		    }
		}
/* L43: */
	    }
L2:
	    z__ = w[k];
	    if (l == k) {
		if (z__ < 0.f) {
		    w[k] = -z__;
		    i__1 = *n;
		    for (j = 1; j <= i__1; ++j) {
			v[j + k * v_dim1] = -v[j + k * v_dim1];
/* L44: */
		    }
		}
		goto L3;
	    }
/*          IF (ITS.EQ.30) PAUSE 'No convergence in 30 iterations' */
	    if (its == 30) {
		return 0;
	    }
	    x = w[l];
	    nm = k - 1;
	    y = w[nm];
	    g = rv1[nm - 1];
	    h__ = rv1[k - 1];
	    f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2.f *
		     y);
	    g = sqrt(f * f + 1.f);
	    f = ((x - z__) * (x + z__) + h__ * (y / (f + r_sign(&g, &f)) - 
		    h__)) / x;
	    c__ = 1.f;
	    s = 1.f;
	    i__1 = nm;
	    for (j = l; j <= i__1; ++j) {
		i__ = j + 1;
		g = rv1[i__ - 1];
		y = w[i__];
		h__ = s * g;
		g = c__ * g;
		z__ = sqrt(f * f + h__ * h__);
		rv1[j - 1] = z__;
		if (z__ == 0.f) {
		    c__ = 0.f;
		    s = 0.f;
		} else {
		    c__ = f / z__;
		    s = h__ / z__;
		}
		f = x * c__ + g * s;
		g = -(x * s) + g * c__;
		h__ = y * s;
		y *= c__;
		i__2 = *n;
		for (nm = 1; nm <= i__2; ++nm) {
		    x = v[nm + j * v_dim1];
		    z__ = v[nm + i__ * v_dim1];
		    v[nm + j * v_dim1] = x * c__ + z__ * s;
		    v[nm + i__ * v_dim1] = -(x * s) + z__ * c__;
/* L45: */
		}
		z__ = sqrt(f * f + h__ * h__);
		w[j] = z__;
		if (z__ != 0.f) {
		    z__ = 1.f / z__;
		    c__ = f * z__;
		    s = h__ * z__;
		}
		f = c__ * g + s * y;
		x = -(s * g) + c__ * y;
		i__2 = *m;
		for (nm = 1; nm <= i__2; ++nm) {
		    y = a[nm + j * a_dim1];
		    z__ = a[nm + i__ * a_dim1];
		    a[nm + j * a_dim1] = y * c__ + z__ * s;
		    a[nm + i__ * a_dim1] = -(y * s) + z__ * c__;
/* L46: */
		}
/* L47: */
	    }
	    rv1[l - 1] = 0.f;
	    rv1[k - 1] = f;
	    w[k] = x;
/* L48: */
	}
L3:
/* L49: */
	;
    }
/*     do 222 i=1,NP */
/*     write (*,*) 'W<x>:',W(i) */
/*     222 continue */
    return 0;
} /* svdcmp_ */