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

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