#include "f2c.h"

/* Subroutine */ int zlartg_(doublecomplex *f, doublecomplex *g, doublereal *

	cs, doublecomplex *sn, doublecomplex *r)

{

/*  -- LAPACK auxiliary routine (version 2.0) --   

       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   

       Courant Institute, Argonne National Lab, and Rice University   

       September 30, 1994   

    Purpose   

    =======   

    ZLARTG generates a plane rotation so that   

       [  CS  SN  ]     [ F ]     [ R ]   

       [  __      ]  .  [   ]  =  [   ]   where CS**2 + |SN|**2 = 1.   

       [ -SN  CS  ]     [ G ]     [ 0 ]   

    This is a faster version of the BLAS1 routine ZROTG, except for   

    the following differences:   

       F and G are unchanged on return.   

       If G=0, then CS=1 and SN=0.   

       If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any   

          floating point operations.   

    Arguments   

    =========   

    F       (input) COMPLEX*16   

            The first component of vector to be rotated.   

    G       (input) COMPLEX*16   

            The second component of vector to be rotated.   

    CS      (output) DOUBLE PRECISION   

            The cosine of the rotation.   

    SN      (output) COMPLEX*16   

            The sine of the rotation.   

    R       (output) COMPLEX*16   

            The nonzero component of the rotated vector.   

    ===================================================================== 

       [ 25 or 38 ops for main paths ] */

    /* System generated locals */

    doublereal d__1, d__2;

    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */

    void d_cnjg(doublecomplex *, doublecomplex *);

    double z_abs(doublecomplex *), d_imag(doublecomplex *), sqrt(doublereal);

    /* Local variables */

    static doublereal d, f1, f2, g1, g2, fa, ga, di;

    static doublecomplex fs, gs, ss;

    if (g->r == 0. && g->i == 0.) {

	*cs = 1.;

	sn->r = 0., sn->i = 0.;

	r->r = f->r, r->i = f->i;

    } else if (f->r == 0. && f->i == 0.) {

	*cs = 0.;

	d_cnjg(&z__2, g);

	d__1 = z_abs(g);

	z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;

	sn->r = z__1.r, sn->i = z__1.i;

	d__1 = z_abs(g);

	r->r = d__1, r->i = 0.;

/*         SN = ONE   

           R = G */

    } else {

	f1 = (d__1 = f->r, abs(d__1)) + (d__2 = d_imag(f), abs(d__2));

	g1 = (d__1 = g->r, abs(d__1)) + (d__2 = d_imag(g), abs(d__2));

	if (f1 >= g1) {

	    z__1.r = g->r / f1, z__1.i = g->i / f1;

	    gs.r = z__1.r, gs.i = z__1.i;

/* Computing 2nd power */

	    d__1 = gs.r;

/* Computing 2nd power */

	    d__2 = d_imag(&gs);

	    g2 = d__1 * d__1 + d__2 * d__2;

	    z__1.r = f->r / f1, z__1.i = f->i / f1;

	    fs.r = z__1.r, fs.i = z__1.i;

/* Computing 2nd power */

	    d__1 = fs.r;

/* Computing 2nd power */

	    d__2 = d_imag(&fs);

	    f2 = d__1 * d__1 + d__2 * d__2;

	    d = sqrt(g2 / f2 + 1.);

	    *cs = 1. / d;

	    d_cnjg(&z__3, &gs);

	    z__2.r = z__3.r * fs.r - z__3.i * fs.i, z__2.i = z__3.r * fs.i + 

		    z__3.i * fs.r;

	    d__1 = *cs / f2;

	    z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;

	    sn->r = z__1.r, sn->i = z__1.i;

	    z__1.r = d * f->r, z__1.i = d * f->i;

	    r->r = z__1.r, r->i = z__1.i;

	} else {

	    z__1.r = f->r / g1, z__1.i = f->i / g1;

	    fs.r = z__1.r, fs.i = z__1.i;

/* Computing 2nd power */

	    d__1 = fs.r;

/* Computing 2nd power */

	    d__2 = d_imag(&fs);

	    f2 = d__1 * d__1 + d__2 * d__2;

	    fa = sqrt(f2);

	    z__1.r = g->r / g1, z__1.i = g->i / g1;

	    gs.r = z__1.r, gs.i = z__1.i;

/* Computing 2nd power */

	    d__1 = gs.r;

/* Computing 2nd power */

	    d__2 = d_imag(&gs);

	    g2 = d__1 * d__1 + d__2 * d__2;

	    ga = sqrt(g2);

	    d = sqrt(f2 / g2 + 1.);

	    di = 1. / d;

	    *cs = fa / ga * di;

	    d_cnjg(&z__3, &gs);

	    z__2.r = z__3.r * fs.r - z__3.i * fs.i, z__2.i = z__3.r * fs.i + 

		    z__3.i * fs.r;

	    d__1 = fa * ga;

	    z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;

	    ss.r = z__1.r, ss.i = z__1.i;

	    z__1.r = di * ss.r, z__1.i = di * ss.i;

	    sn->r = z__1.r, sn->i = z__1.i;

	    z__2.r = g->r * ss.r - g->i * ss.i, z__2.i = g->r * ss.i + g->i * 

		    ss.r;

	    z__1.r = d * z__2.r, z__1.i = d * z__2.i;

	    r->r = z__1.r, r->i = z__1.i;

	}

    }

    return 0;

/*     End of ZLARTG */

} /* zlartg_ */