#include "f2c.h"

/* Subroutine */ int clartg_(complex *f, complex *g, real *cs, complex *sn, 

	complex *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   

    =======   

    CLARTG 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 CROTG, 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   

            The first component of vector to be rotated.   

    G       (input) COMPLEX   

            The second component of vector to be rotated.   

    CS      (output) REAL   

            The cosine of the rotation.   

    SN      (output) COMPLEX   

            The sine of the rotation.   

    R       (output) COMPLEX   

            The nonzero component of the rotated vector.   

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

       [ 25 or 38 ops for main paths ] */

    /* System generated locals */

    real r__1, r__2;

    doublereal d__1;

    complex q__1, q__2, q__3;

    /* Builtin functions */

    void r_cnjg(complex *, complex *);

    double c_abs(complex *), r_imag(complex *), sqrt(doublereal);

    /* Local variables */

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

    static complex fs, gs, ss;

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

	*cs = 1.f;

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

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

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

	*cs = 0.f;

	r_cnjg(&q__2, g);

	d__1 = c_abs(g);

	q__1.r = q__2.r / d__1, q__1.i = q__2.i / d__1;

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

	d__1 = c_abs(g);

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

/*         SN = ONE   

           R = G */

    } else {

	f1 = (r__1 = f->r, dabs(r__1)) + (r__2 = r_imag(f), dabs(r__2));

	g1 = (r__1 = g->r, dabs(r__1)) + (r__2 = r_imag(g), dabs(r__2));

	if (f1 >= g1) {

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

	    gs.r = q__1.r, gs.i = q__1.i;

/* Computing 2nd power */

	    r__1 = gs.r;

/* Computing 2nd power */

	    r__2 = r_imag(&gs);

	    g2 = r__1 * r__1 + r__2 * r__2;

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

	    fs.r = q__1.r, fs.i = q__1.i;

/* Computing 2nd power */

	    r__1 = fs.r;

/* Computing 2nd power */

	    r__2 = r_imag(&fs);

	    f2 = r__1 * r__1 + r__2 * r__2;

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

	    *cs = 1.f / d;

	    r_cnjg(&q__3, &gs);

	    q__2.r = q__3.r * fs.r - q__3.i * fs.i, q__2.i = q__3.r * fs.i + 

		    q__3.i * fs.r;

	    d__1 = *cs / f2;

	    q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;

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

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

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

	} else {

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

	    fs.r = q__1.r, fs.i = q__1.i;

/* Computing 2nd power */

	    r__1 = fs.r;

/* Computing 2nd power */

	    r__2 = r_imag(&fs);

	    f2 = r__1 * r__1 + r__2 * r__2;

	    fa = sqrt(f2);

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

	    gs.r = q__1.r, gs.i = q__1.i;

/* Computing 2nd power */

	    r__1 = gs.r;

/* Computing 2nd power */

	    r__2 = r_imag(&gs);

	    g2 = r__1 * r__1 + r__2 * r__2;

	    ga = sqrt(g2);

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

	    di = 1.f / d;

	    *cs = fa / ga * di;

	    r_cnjg(&q__3, &gs);

	    q__2.r = q__3.r * fs.r - q__3.i * fs.i, q__2.i = q__3.r * fs.i + 

		    q__3.i * fs.r;

	    d__1 = fa * ga;

	    q__1.r = q__2.r / d__1, q__1.i = q__2.i / d__1;

	    ss.r = q__1.r, ss.i = q__1.i;

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

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

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

		    ss.r;

	    q__1.r = d * q__2.r, q__1.i = d * q__2.i;

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

	}

    }

    return 0;

/*     End of CLARTG */

} /* clartg_ */