dlaev2.c
Go to the documentation of this file.
00001 /* dlaev2.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int dlaev2_(doublereal *a, doublereal *b, doublereal *c__, 
00017         doublereal *rt1, doublereal *rt2, doublereal *cs1, doublereal *sn1)
00018 {
00019     /* System generated locals */
00020     doublereal d__1;
00021 
00022     /* Builtin functions */
00023     double sqrt(doublereal);
00024 
00025     /* Local variables */
00026     doublereal ab, df, cs, ct, tb, sm, tn, rt, adf, acs;
00027     integer sgn1, sgn2;
00028     doublereal acmn, acmx;
00029 
00030 
00031 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00032 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00033 /*     November 2006 */
00034 
00035 /*     .. Scalar Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix */
00042 /*     [  A   B  ] */
00043 /*     [  B   C  ]. */
00044 /*  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the */
00045 /*  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right */
00046 /*  eigenvector for RT1, giving the decomposition */
00047 
00048 /*     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ] */
00049 /*     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ]. */
00050 
00051 /*  Arguments */
00052 /*  ========= */
00053 
00054 /*  A       (input) DOUBLE PRECISION */
00055 /*          The (1,1) element of the 2-by-2 matrix. */
00056 
00057 /*  B       (input) DOUBLE PRECISION */
00058 /*          The (1,2) element and the conjugate of the (2,1) element of */
00059 /*          the 2-by-2 matrix. */
00060 
00061 /*  C       (input) DOUBLE PRECISION */
00062 /*          The (2,2) element of the 2-by-2 matrix. */
00063 
00064 /*  RT1     (output) DOUBLE PRECISION */
00065 /*          The eigenvalue of larger absolute value. */
00066 
00067 /*  RT2     (output) DOUBLE PRECISION */
00068 /*          The eigenvalue of smaller absolute value. */
00069 
00070 /*  CS1     (output) DOUBLE PRECISION */
00071 /*  SN1     (output) DOUBLE PRECISION */
00072 /*          The vector (CS1, SN1) is a unit right eigenvector for RT1. */
00073 
00074 /*  Further Details */
00075 /*  =============== */
00076 
00077 /*  RT1 is accurate to a few ulps barring over/underflow. */
00078 
00079 /*  RT2 may be inaccurate if there is massive cancellation in the */
00080 /*  determinant A*C-B*B; higher precision or correctly rounded or */
00081 /*  correctly truncated arithmetic would be needed to compute RT2 */
00082 /*  accurately in all cases. */
00083 
00084 /*  CS1 and SN1 are accurate to a few ulps barring over/underflow. */
00085 
00086 /*  Overflow is possible only if RT1 is within a factor of 5 of overflow. */
00087 /*  Underflow is harmless if the input data is 0 or exceeds */
00088 /*     underflow_threshold / macheps. */
00089 
00090 /* ===================================================================== */
00091 
00092 /*     .. Parameters .. */
00093 /*     .. */
00094 /*     .. Local Scalars .. */
00095 /*     .. */
00096 /*     .. Intrinsic Functions .. */
00097 /*     .. */
00098 /*     .. Executable Statements .. */
00099 
00100 /*     Compute the eigenvalues */
00101 
00102     sm = *a + *c__;
00103     df = *a - *c__;
00104     adf = abs(df);
00105     tb = *b + *b;
00106     ab = abs(tb);
00107     if (abs(*a) > abs(*c__)) {
00108         acmx = *a;
00109         acmn = *c__;
00110     } else {
00111         acmx = *c__;
00112         acmn = *a;
00113     }
00114     if (adf > ab) {
00115 /* Computing 2nd power */
00116         d__1 = ab / adf;
00117         rt = adf * sqrt(d__1 * d__1 + 1.);
00118     } else if (adf < ab) {
00119 /* Computing 2nd power */
00120         d__1 = adf / ab;
00121         rt = ab * sqrt(d__1 * d__1 + 1.);
00122     } else {
00123 
00124 /*        Includes case AB=ADF=0 */
00125 
00126         rt = ab * sqrt(2.);
00127     }
00128     if (sm < 0.) {
00129         *rt1 = (sm - rt) * .5;
00130         sgn1 = -1;
00131 
00132 /*        Order of execution important. */
00133 /*        To get fully accurate smaller eigenvalue, */
00134 /*        next line needs to be executed in higher precision. */
00135 
00136         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
00137     } else if (sm > 0.) {
00138         *rt1 = (sm + rt) * .5;
00139         sgn1 = 1;
00140 
00141 /*        Order of execution important. */
00142 /*        To get fully accurate smaller eigenvalue, */
00143 /*        next line needs to be executed in higher precision. */
00144 
00145         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
00146     } else {
00147 
00148 /*        Includes case RT1 = RT2 = 0 */
00149 
00150         *rt1 = rt * .5;
00151         *rt2 = rt * -.5;
00152         sgn1 = 1;
00153     }
00154 
00155 /*     Compute the eigenvector */
00156 
00157     if (df >= 0.) {
00158         cs = df + rt;
00159         sgn2 = 1;
00160     } else {
00161         cs = df - rt;
00162         sgn2 = -1;
00163     }
00164     acs = abs(cs);
00165     if (acs > ab) {
00166         ct = -tb / cs;
00167         *sn1 = 1. / sqrt(ct * ct + 1.);
00168         *cs1 = ct * *sn1;
00169     } else {
00170         if (ab == 0.) {
00171             *cs1 = 1.;
00172             *sn1 = 0.;
00173         } else {
00174             tn = -cs / tb;
00175             *cs1 = 1. / sqrt(tn * tn + 1.);
00176             *sn1 = tn * *cs1;
00177         }
00178     }
00179     if (sgn1 == sgn2) {
00180         tn = *cs1;
00181         *cs1 = -(*sn1);
00182         *sn1 = tn;
00183     }
00184     return 0;
00185 
00186 /*     End of DLAEV2 */
00187 
00188 } /* dlaev2_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:46