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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlas2_(doublereal *f, doublereal *g, doublereal *h__,
   4         doublereal *ssmin, doublereal *ssmax)
   5 {
   6     /* System generated locals */
   7     doublereal d__1, d__2;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     doublereal c__, fa, ga, ha, as, at, au, fhmn, fhmx;
  14
  15
  16 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  17 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  18 /*     November 2006 */
  19
  20 /*     .. Scalar Arguments .. */
  21 /*     .. */
  22
  23 /*  Purpose */
  24 /*  ======= */
  25
  26 /*  DLAS2  computes the singular values of the 2-by-2 matrix */
  27 /*     [  F   G  ] */
  28 /*     [  0   H  ]. */
  29 /*  On return, SSMIN is the smaller singular value and SSMAX is the */
  30 /*  larger singular value. */
  31
  32 /*  Arguments */
  33 /*  ========= */
  34
  35 /*  F       (input) DOUBLE PRECISION */
  36 /*          The (1,1) element of the 2-by-2 matrix. */
  37
  38 /*  G       (input) DOUBLE PRECISION */
  39 /*          The (1,2) element of the 2-by-2 matrix. */
  40
  41 /*  H       (input) DOUBLE PRECISION */
  42 /*          The (2,2) element of the 2-by-2 matrix. */
  43
  44 /*  SSMIN   (output) DOUBLE PRECISION */
  45 /*          The smaller singular value. */
  46
  47 /*  SSMAX   (output) DOUBLE PRECISION */
  48 /*          The larger singular value. */
  49
  50 /*  Further Details */
  51 /*  =============== */
  52
  53 /*  Barring over/underflow, all output quantities are correct to within */
  54 /*  a few units in the last place (ulps), even in the absence of a guard */
  55 /*  digit in addition/subtraction. */
  56
  57 /*  In IEEE arithmetic, the code works correctly if one matrix element is */
  58 /*  infinite. */
  59
  60 /*  Overflow will not occur unless the largest singular value itself */
  61 /*  overflows, or is within a few ulps of overflow. (On machines with */
  62 /*  partial overflow, like the Cray, overflow may occur if the largest */
  63 /*  singular value is within a factor of 2 of overflow.) */
  64
  65 /*  Underflow is harmless if underflow is gradual. Otherwise, results */
  66 /*  may correspond to a matrix modified by perturbations of size near */
  67 /*  the underflow threshold. */
  68
  69 /*  ==================================================================== */
  70
  71 /*     .. Parameters .. */
  72 /*     .. */
  73 /*     .. Local Scalars .. */
  74 /*     .. */
  75 /*     .. Intrinsic Functions .. */
  76 /*     .. */
  77 /*     .. Executable Statements .. */
  78
  79     fa = abs(*f);
  80     ga = abs(*g);
  81     ha = abs(*h__);
  82     fhmn = min(fa,ha);
  83     fhmx = max(fa,ha);
  84     if (fhmn == 0.) {
  85         *ssmin = 0.;
  86         if (fhmx == 0.) {
  87             *ssmax = ga;
  88         } else {
  89 /* Computing 2nd power */
  90             d__1 = min(fhmx,ga) / max(fhmx,ga);
  91             *ssmax = max(fhmx,ga) * sqrt(d__1 * d__1 + 1.);
  92         }
  93     } else {
  94         if (ga < fhmx) {
  95             as = fhmn / fhmx + 1.;
  96             at = (fhmx - fhmn) / fhmx;
  97 /* Computing 2nd power */
  98             d__1 = ga / fhmx;
  99             au = d__1 * d__1;
 100             c__ = 2. / (sqrt(as * as + au) + sqrt(at * at + au));
 101             *ssmin = fhmn * c__;
 102             *ssmax = fhmx / c__;
 103         } else {
 104             au = fhmx / ga;
 105             if (au == 0.) {
 106
 107 /*              Avoid possible harmful underflow if exponent range */
 108 /*              asymmetric (true SSMIN may not underflow even if */
 109 /*              AU underflows) */
 110
 111                 *ssmin = fhmn * fhmx / ga;
 112                 *ssmax = ga;
 113             } else {
 114                 as = fhmn / fhmx + 1.;
 115                 at = (fhmx - fhmn) / fhmx;
 116 /* Computing 2nd power */
 117                 d__1 = as * au;
 118 /* Computing 2nd power */
 119                 d__2 = at * au;
 120                 c__ = 1. / (sqrt(d__1 * d__1 + 1.) + sqrt(d__2 * d__2 + 1.));
 121                 *ssmin = fhmn * c__ * au;
 122                 *ssmin += *ssmin;
 123                 *ssmax = ga / (c__ + c__);
 124             }
 125         }
 126     }
 127     return 0;
 128
 129 /*     End of DLAS2 */
 130
 131 } /* dlas2_ */