zlartg.c

SuperLU is a general purpose library for the direct solution of large, sparse, nonsymmetric systems

#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_ */