/* dlanv2.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "blaswrap.h"
/* Table of constant values */
static doublereal c_b4 = 1.;
/* Subroutine */ int dlanv2_(doublereal *a, doublereal *b, doublereal *c__,
doublereal *d__, doublereal *rt1r, doublereal *rt1i, doublereal *rt2r,
doublereal *rt2i, doublereal *cs, doublereal *sn)
{
/* System generated locals */
doublereal d__1, d__2;
/* Builtin functions */
double d_sign(doublereal *, doublereal *), sqrt(doublereal);
/* Local variables */
doublereal p, z__, aa, bb, cc, dd, cs1, sn1, sab, sac, eps, tau, temp,
scale, bcmax, bcmis, sigma;
extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
/* -- LAPACK driver routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric */
/* matrix in standard form: */
/* [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ] */
/* [ C D ] [ SN CS ] [ CC DD ] [-SN CS ] */
/* where either */
/* 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or */
/* 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex */
/* conjugate eigenvalues. */
/* Arguments */
/* ========= */
/* A (input/output) DOUBLE PRECISION */
/* B (input/output) DOUBLE PRECISION */
/* C (input/output) DOUBLE PRECISION */
/* D (input/output) DOUBLE PRECISION */
/* On entry, the elements of the input matrix. */
/* On exit, they are overwritten by the elements of the */
/* standardised Schur form. */
/* RT1R (output) DOUBLE PRECISION */
/* RT1I (output) DOUBLE PRECISION */
/* RT2R (output) DOUBLE PRECISION */
/* RT2I (output) DOUBLE PRECISION */
/* The real and imaginary parts of the eigenvalues. If the */
/* eigenvalues are a complex conjugate pair, RT1I > 0. */
/* CS (output) DOUBLE PRECISION */
/* SN (output) DOUBLE PRECISION */
/* Parameters of the rotation matrix. */
/* Further Details */
/* =============== */
/* Modified by V. Sima, Research Institute for Informatics, Bucharest, */
/* Romania, to reduce the risk of cancellation errors, */
/* when computing real eigenvalues, and to ensure, if possible, that */
/* abs(RT1R) >= abs(RT2R). */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
eps = dlamch_("P");
if (*c__ == 0.) {
*cs = 1.;
*sn = 0.;
goto L10;
} else if (*b == 0.) {
/* Swap rows and columns */
*cs = 0.;
*sn = 1.;
temp = *d__;
*d__ = *a;
*a = temp;
*b = -(*c__);
*c__ = 0.;
goto L10;
} else if (*a - *d__ == 0. && d_sign(&c_b4, b) != d_sign(&c_b4, c__)) {
*cs = 1.;
*sn = 0.;
goto L10;
} else {
temp = *a - *d__;
p = temp * .5;
/* Computing MAX */
d__1 = abs(*b), d__2 = abs(*c__);
bcmax = max(d__1,d__2);
/* Computing MIN */
d__1 = abs(*b), d__2 = abs(*c__);
bcmis = min(d__1,d__2) * d_sign(&c_b4, b) * d_sign(&c_b4, c__);
/* Computing MAX */
d__1 = abs(p);
scale = max(d__1,bcmax);
z__ = p / scale * p + bcmax / scale * bcmis;
/* If Z is of the order of the machine accuracy, postpone the */
/* decision on the nature of eigenvalues */
if (z__ >= eps * 4.) {
/* Real eigenvalues. Compute A and D. */
d__1 = sqrt(scale) * sqrt(z__);
z__ = p + d_sign(&d__1, &p);
*a = *d__ + z__;
*d__ -= bcmax / z__ * bcmis;
/* Compute B and the rotation matrix */
tau = dlapy2_(c__, &z__);
*cs = z__ / tau;
*sn = *c__ / tau;
*b -= *c__;
*c__ = 0.;
} else {
/* Complex eigenvalues, or real (almost) equal eigenvalues. */
/* Make diagonal elements equal. */
sigma = *b + *c__;
tau = dlapy2_(&sigma, &temp);
*cs = sqrt((abs(sigma) / tau + 1.) * .5);
*sn = -(p / (tau * *cs)) * d_sign(&c_b4, &sigma);
/* Compute [ AA BB ] = [ A B ] [ CS -SN ] */
/* [ CC DD ] [ C D ] [ SN CS ] */
aa = *a * *cs + *b * *sn;
bb = -(*a) * *sn + *b * *cs;
cc = *c__ * *cs + *d__ * *sn;
dd = -(*c__) * *sn + *d__ * *cs;
/* Compute [ A B ] = [ CS SN ] [ AA BB ] */
/* [ C D ] [-SN CS ] [ CC DD ] */
*a = aa * *cs + cc * *sn;
*b = bb * *cs + dd * *sn;
*c__ = -aa * *sn + cc * *cs;
*d__ = -bb * *sn + dd * *cs;
temp = (*a + *d__) * .5;
*a = temp;
*d__ = temp;
if (*c__ != 0.) {
if (*b != 0.) {
if (d_sign(&c_b4, b) == d_sign(&c_b4, c__)) {
/* Real eigenvalues: reduce to upper triangular form */
sab = sqrt((abs(*b)));
sac = sqrt((abs(*c__)));
d__1 = sab * sac;
p = d_sign(&d__1, c__);
tau = 1. / sqrt((d__1 = *b + *c__, abs(d__1)));
*a = temp + p;
*d__ = temp - p;
*b -= *c__;
*c__ = 0.;
cs1 = sab * tau;
sn1 = sac * tau;
temp = *cs * cs1 - *sn * sn1;
*sn = *cs * sn1 + *sn * cs1;
*cs = temp;
}
} else {
*b = -(*c__);
*c__ = 0.;
temp = *cs;
*cs = -(*sn);
*sn = temp;
}
}
}
}
L10:
/* Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I). */
*rt1r = *a;
*rt2r = *d__;
if (*c__ == 0.) {
*rt1i = 0.;
*rt2i = 0.;
} else {
*rt1i = sqrt((abs(*b))) * sqrt((abs(*c__)));
*rt2i = -(*rt1i);
}
return 0;
/* End of DLANV2 */
} /* dlanv2_ */