git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlaev2.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlaev2_(doublereal *a, doublereal *b, doublereal *c__,
   4         doublereal *rt1, doublereal *rt2, doublereal *cs1, doublereal *sn1)
   5 {
   6     /* System generated locals */
   7     doublereal d__1;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     doublereal ab, df, cs, ct, tb, sm, tn, rt, adf, acs;
  14     integer sgn1, sgn2;
  15     doublereal acmn, acmx;
  16
  17
  18 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  19 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  20 /*     November 2006 */
  21
  22 /*     .. Scalar Arguments .. */
  23 /*     .. */
  24
  25 /*  Purpose */
  26 /*  ======= */
  27
  28 /*  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix */
  29 /*     [  A   B  ] */
  30 /*     [  B   C  ]. */
  31 /*  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the */
  32 /*  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right */
  33 /*  eigenvector for RT1, giving the decomposition */
  34
  35 /*     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ] */
  36 /*     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ]. */
  37
  38 /*  Arguments */
  39 /*  ========= */
  40
  41 /*  A       (input) DOUBLE PRECISION */
  42 /*          The (1,1) element of the 2-by-2 matrix. */
  43
  44 /*  B       (input) DOUBLE PRECISION */
  45 /*          The (1,2) element and the conjugate of the (2,1) element of */
  46 /*          the 2-by-2 matrix. */
  47
  48 /*  C       (input) DOUBLE PRECISION */
  49 /*          The (2,2) element of the 2-by-2 matrix. */
  50
  51 /*  RT1     (output) DOUBLE PRECISION */
  52 /*          The eigenvalue of larger absolute value. */
  53
  54 /*  RT2     (output) DOUBLE PRECISION */
  55 /*          The eigenvalue of smaller absolute value. */
  56
  57 /*  CS1     (output) DOUBLE PRECISION */
  58 /*  SN1     (output) DOUBLE PRECISION */
  59 /*          The vector (CS1, SN1) is a unit right eigenvector for RT1. */
  60
  61 /*  Further Details */
  62 /*  =============== */
  63
  64 /*  RT1 is accurate to a few ulps barring over/underflow. */
  65
  66 /*  RT2 may be inaccurate if there is massive cancellation in the */
  67 /*  determinant A*C-B*B; higher precision or correctly rounded or */
  68 /*  correctly truncated arithmetic would be needed to compute RT2 */
  69 /*  accurately in all cases. */
  70
  71 /*  CS1 and SN1 are accurate to a few ulps barring over/underflow. */
  72
  73 /*  Overflow is possible only if RT1 is within a factor of 5 of overflow. */
  74 /*  Underflow is harmless if the input data is 0 or exceeds */
  75 /*     underflow_threshold / macheps. */
  76
  77 /* ===================================================================== */
  78
  79 /*     .. Parameters .. */
  80 /*     .. */
  81 /*     .. Local Scalars .. */
  82 /*     .. */
  83 /*     .. Intrinsic Functions .. */
  84 /*     .. */
  85 /*     .. Executable Statements .. */
  86
  87 /*     Compute the eigenvalues */
  88
  89     sm = *a + *c__;
  90     df = *a - *c__;
  91     adf = abs(df);
  92     tb = *b + *b;
  93     ab = abs(tb);
  94     if (abs(*a) > abs(*c__)) {
  95         acmx = *a;
  96         acmn = *c__;
  97     } else {
  98         acmx = *c__;
  99         acmn = *a;
 100     }
 101     if (adf > ab) {
 102 /* Computing 2nd power */
 103         d__1 = ab / adf;
 104         rt = adf * sqrt(d__1 * d__1 + 1.);
 105     } else if (adf < ab) {
 106 /* Computing 2nd power */
 107         d__1 = adf / ab;
 108         rt = ab * sqrt(d__1 * d__1 + 1.);
 109     } else {
 110
 111 /*        Includes case AB=ADF=0 */
 112
 113         rt = ab * sqrt(2.);
 114     }
 115     if (sm < 0.) {
 116         *rt1 = (sm - rt) * .5;
 117         sgn1 = -1;
 118
 119 /*        Order of execution important. */
 120 /*        To get fully accurate smaller eigenvalue, */
 121 /*        next line needs to be executed in higher precision. */
 122
 123         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
 124     } else if (sm > 0.) {
 125         *rt1 = (sm + rt) * .5;
 126         sgn1 = 1;
 127
 128 /*        Order of execution important. */
 129 /*        To get fully accurate smaller eigenvalue, */
 130 /*        next line needs to be executed in higher precision. */
 131
 132         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
 133     } else {
 134
 135 /*        Includes case RT1 = RT2 = 0 */
 136
 137         *rt1 = rt * .5;
 138         *rt2 = rt * -.5;
 139         sgn1 = 1;
 140     }
 141
 142 /*     Compute the eigenvector */
 143
 144     if (df >= 0.) {
 145         cs = df + rt;
 146         sgn2 = 1;
 147     } else {
 148         cs = df - rt;
 149         sgn2 = -1;
 150     }
 151     acs = abs(cs);
 152     if (acs > ab) {
 153         ct = -tb / cs;
 154         *sn1 = 1. / sqrt(ct * ct + 1.);
 155         *cs1 = ct * *sn1;
 156     } else {
 157         if (ab == 0.) {
 158             *cs1 = 1.;
 159             *sn1 = 0.;
 160         } else {
 161             tn = -cs / tb;
 162             *cs1 = 1. / sqrt(tn * tn + 1.);
 163             *sn1 = tn * *cs1;
 164         }
 165     }
 166     if (sgn1 == sgn2) {
 167         tn = *cs1;
 168         *cs1 = -(*sn1);
 169         *sn1 = tn;
 170     }
 171     return 0;
 172
 173 /*     End of DLAEV2 */
 174
 175 } /* dlaev2_ */