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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static real c_b3 = 2.f;
   6 static real c_b4 = 1.f;
   7
   8 /* Subroutine */ int slasv2_(real *f, real *g, real *h__, real *ssmin, real *
   9         ssmax, real *snr, real *csr, real *snl, real *csl)
  10 {
  11     /* System generated locals */
  12     real r__1;
  13
  14     /* Builtin functions */
  15     double sqrt(doublereal), r_sign(real *, real *);
  16
  17     /* Local variables */
  18     real a, d__, l, m, r__, s, t, fa, ga, ha, ft, gt, ht, mm, tt, clt, crt,
  19             slt, srt;
  20     integer pmax;
  21     real temp;
  22     logical swap;
  23     real tsign;
  24     logical gasmal;
  25     extern doublereal slamch_(char *);
  26
  27
  28 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  29 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  30 /*     November 2006 */
  31
  32 /*     .. Scalar Arguments .. */
  33 /*     .. */
  34
  35 /*  Purpose */
  36 /*  ======= */
  37
  38 /*  SLASV2 computes the singular value decomposition of a 2-by-2 */
  39 /*  triangular matrix */
  40 /*     [  F   G  ] */
  41 /*     [  0   H  ]. */
  42 /*  On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the */
  43 /*  smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and */
  44 /*  right singular vectors for abs(SSMAX), giving the decomposition */
  45
  46 /*     [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ] */
  47 /*     [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ]. */
  48
  49 /*  Arguments */
  50 /*  ========= */
  51
  52 /*  F       (input) REAL */
  53 /*          The (1,1) element of the 2-by-2 matrix. */
  54
  55 /*  G       (input) REAL */
  56 /*          The (1,2) element of the 2-by-2 matrix. */
  57
  58 /*  H       (input) REAL */
  59 /*          The (2,2) element of the 2-by-2 matrix. */
  60
  61 /*  SSMIN   (output) REAL */
  62 /*          abs(SSMIN) is the smaller singular value. */
  63
  64 /*  SSMAX   (output) REAL */
  65 /*          abs(SSMAX) is the larger singular value. */
  66
  67 /*  SNL     (output) REAL */
  68 /*  CSL     (output) REAL */
  69 /*          The vector (CSL, SNL) is a unit left singular vector for the */
  70 /*          singular value abs(SSMAX). */
  71
  72 /*  SNR     (output) REAL */
  73 /*  CSR     (output) REAL */
  74 /*          The vector (CSR, SNR) is a unit right singular vector for the */
  75 /*          singular value abs(SSMAX). */
  76
  77 /*  Further Details */
  78 /*  =============== */
  79
  80 /*  Any input parameter may be aliased with any output parameter. */
  81
  82 /*  Barring over/underflow and assuming a guard digit in subtraction, all */
  83 /*  output quantities are correct to within a few units in the last */
  84 /*  place (ulps). */
  85
  86 /*  In IEEE arithmetic, the code works correctly if one matrix element is */
  87 /*  infinite. */
  88
  89 /*  Overflow will not occur unless the largest singular value itself */
  90 /*  overflows or is within a few ulps of overflow. (On machines with */
  91 /*  partial overflow, like the Cray, overflow may occur if the largest */
  92 /*  singular value is within a factor of 2 of overflow.) */
  93
  94 /*  Underflow is harmless if underflow is gradual. Otherwise, results */
  95 /*  may correspond to a matrix modified by perturbations of size near */
  96 /*  the underflow threshold. */
  97
  98 /* ===================================================================== */
  99
 100 /*     .. Parameters .. */
 101 /*     .. */
 102 /*     .. Local Scalars .. */
 103 /*     .. */
 104 /*     .. Intrinsic Functions .. */
 105 /*     .. */
 106 /*     .. External Functions .. */
 107 /*     .. */
 108 /*     .. Executable Statements .. */
 109
 110     ft = *f;
 111     fa = dabs(ft);
 112     ht = *h__;
 113     ha = dabs(*h__);
 114
 115 /*     PMAX points to the maximum absolute element of matrix */
 116 /*       PMAX = 1 if F largest in absolute values */
 117 /*       PMAX = 2 if G largest in absolute values */
 118 /*       PMAX = 3 if H largest in absolute values */
 119
 120     pmax = 1;
 121     swap = ha > fa;
 122     if (swap) {
 123         pmax = 3;
 124         temp = ft;
 125         ft = ht;
 126         ht = temp;
 127         temp = fa;
 128         fa = ha;
 129         ha = temp;
 130
 131 /*        Now FA .ge. HA */
 132
 133     }
 134     gt = *g;
 135     ga = dabs(gt);
 136     if (ga == 0.f) {
 137
 138 /*        Diagonal matrix */
 139
 140         *ssmin = ha;
 141         *ssmax = fa;
 142         clt = 1.f;
 143         crt = 1.f;
 144         slt = 0.f;
 145         srt = 0.f;
 146     } else {
 147         gasmal = TRUE_;
 148         if (ga > fa) {
 149             pmax = 2;
 150             if (fa / ga < slamch_("EPS")) {
 151
 152 /*              Case of very large GA */
 153
 154                 gasmal = FALSE_;
 155                 *ssmax = ga;
 156                 if (ha > 1.f) {
 157                     *ssmin = fa / (ga / ha);
 158                 } else {
 159                     *ssmin = fa / ga * ha;
 160                 }
 161                 clt = 1.f;
 162                 slt = ht / gt;
 163                 srt = 1.f;
 164                 crt = ft / gt;
 165             }
 166         }
 167         if (gasmal) {
 168
 169 /*           Normal case */
 170
 171             d__ = fa - ha;
 172             if (d__ == fa) {
 173
 174 /*              Copes with infinite F or H */
 175
 176                 l = 1.f;
 177             } else {
 178                 l = d__ / fa;
 179             }
 180
 181 /*           Note that 0 .le. L .le. 1 */
 182
 183             m = gt / ft;
 184
 185 /*           Note that abs(M) .le. 1/macheps */
 186
 187             t = 2.f - l;
 188
 189 /*           Note that T .ge. 1 */
 190
 191             mm = m * m;
 192             tt = t * t;
 193             s = sqrt(tt + mm);
 194
 195 /*           Note that 1 .le. S .le. 1 + 1/macheps */
 196
 197             if (l == 0.f) {
 198                 r__ = dabs(m);
 199             } else {
 200                 r__ = sqrt(l * l + mm);
 201             }
 202
 203 /*           Note that 0 .le. R .le. 1 + 1/macheps */
 204
 205             a = (s + r__) * .5f;
 206
 207 /*           Note that 1 .le. A .le. 1 + abs(M) */
 208
 209             *ssmin = ha / a;
 210             *ssmax = fa * a;
 211             if (mm == 0.f) {
 212
 213 /*              Note that M is very tiny */
 214
 215                 if (l == 0.f) {
 216                     t = r_sign(&c_b3, &ft) * r_sign(&c_b4, &gt);
 217                 } else {
 218                     t = gt / r_sign(&d__, &ft) + m / t;
 219                 }
 220             } else {
 221                 t = (m / (s + t) + m / (r__ + l)) * (a + 1.f);
 222             }
 223             l = sqrt(t * t + 4.f);
 224             crt = 2.f / l;
 225             srt = t / l;
 226             clt = (crt + srt * m) / a;
 227             slt = ht / ft * srt / a;
 228         }
 229     }
 230     if (swap) {
 231         *csl = srt;
 232         *snl = crt;
 233         *csr = slt;
 234         *snr = clt;
 235     } else {
 236         *csl = clt;
 237         *snl = slt;
 238         *csr = crt;
 239         *snr = srt;
 240     }
 241
 242 /*     Correct signs of SSMAX and SSMIN */
 243
 244     if (pmax == 1) {
 245         tsign = r_sign(&c_b4, csr) * r_sign(&c_b4, csl) * r_sign(&c_b4, f);
 246     }
 247     if (pmax == 2) {
 248         tsign = r_sign(&c_b4, snr) * r_sign(&c_b4, csl) * r_sign(&c_b4, g);
 249     }
 250     if (pmax == 3) {
 251         tsign = r_sign(&c_b4, snr) * r_sign(&c_b4, snl) * r_sign(&c_b4, h__);
 252     }
 253     *ssmax = r_sign(ssmax, &tsign);
 254     r__1 = tsign * r_sign(&c_b4, f) * r_sign(&c_b4, h__);
 255     *ssmin = r_sign(ssmin, &r__1);
 256     return 0;
 257
 258 /*     End of SLASV2 */
 259
 260 } /* slasv2_ */