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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlartg_(doublereal *f, doublereal *g, doublereal *cs,
   4         doublereal *sn, doublereal *r__)
   5 {
   6     /* System generated locals */
   7     integer i__1;
   8     doublereal d__1, d__2;
   9
  10     /* Builtin functions */
  11     double log(doublereal), pow_di(doublereal *, integer *), sqrt(doublereal);
  12
  13     /* Local variables */
  14     integer i__;
  15     doublereal f1, g1, scale;
  16     integer count;
  17     static doublereal safmn2, safmx2;
  18     extern doublereal dlamch_(char *);
  19     static doublereal safmin, eps;
  20     static volatile logical first = TRUE_;
  21
  22 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  23 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  24 /*     November 2006 */
  25
  26 /*     .. Scalar Arguments .. */
  27 /*     .. */
  28
  29 /*  Purpose */
  30 /*  ======= */
  31
  32 /*  DLARTG generate a plane rotation so that */
  33
  34 /*     [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1. */
  35 /*     [ -SN  CS  ]     [ G ]     [ 0 ] */
  36
  37 /*  This is a slower, more accurate version of the BLAS1 routine DROTG, */
  38 /*  with the following other differences: */
  39 /*     F and G are unchanged on return. */
  40 /*     If G=0, then CS=1 and SN=0. */
  41 /*     If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any */
  42 /*        floating point operations (saves work in DBDSQR when */
  43 /*        there are zeros on the diagonal). */
  44
  45 /*  If F exceeds G in magnitude, CS will be positive. */
  46
  47 /*  Arguments */
  48 /*  ========= */
  49
  50 /*  F       (input) DOUBLE PRECISION */
  51 /*          The first component of vector to be rotated. */
  52
  53 /*  G       (input) DOUBLE PRECISION */
  54 /*          The second component of vector to be rotated. */
  55
  56 /*  CS      (output) DOUBLE PRECISION */
  57 /*          The cosine of the rotation. */
  58
  59 /*  SN      (output) DOUBLE PRECISION */
  60 /*          The sine of the rotation. */
  61
  62 /*  R       (output) DOUBLE PRECISION */
  63 /*          The nonzero component of the rotated vector. */
  64
  65 /*  This version has a few statements commented out for thread safety */
  66 /*  (machine parameters are computed on each entry). 10 feb 03, SJH. */
  67
  68 /*  ===================================================================== */
  69
  70 /*     .. Parameters .. */
  71 /*     .. */
  72 /*     .. Local Scalars .. */
  73 /*     LOGICAL            FIRST */
  74 /*     .. */
  75 /*     .. External Functions .. */
  76 /*     .. */
  77 /*     .. Intrinsic Functions .. */
  78 /*     .. */
  79 /*     .. Save statement .. */
  80 /*     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2 */
  81 /*     .. */
  82 /*     .. Data statements .. */
  83 /*     DATA               FIRST / .TRUE. / */
  84 /*     .. */
  85 /*     .. Executable Statements .. */
  86
  87 /*     IF( FIRST ) THEN */
  88     if (first) {
  89         safmin = dlamch_("S");
  90         eps = dlamch_("E");
  91         d__1 = dlamch_("B");
  92         i__1 = (integer) (log(safmin / eps) / log(dlamch_("B")) / 2.);
  93         safmn2 = pow_di(&d__1, &i__1);
  94         safmx2 = 1. / safmn2;
  95         first = FALSE_;
  96     }
  97 /*        FIRST = .FALSE. */
  98 /*     END IF */
  99     if (*g == 0.) {
 100         *cs = 1.;
 101         *sn = 0.;
 102         *r__ = *f;
 103     } else if (*f == 0.) {
 104         *cs = 0.;
 105         *sn = 1.;
 106         *r__ = *g;
 107     } else {
 108         f1 = *f;
 109         g1 = *g;
 110 /* Computing MAX */
 111         d__1 = abs(f1), d__2 = abs(g1);
 112         scale = max(d__1,d__2);
 113         if (scale >= safmx2) {
 114             count = 0;
 115 L10:
 116             ++count;
 117             f1 *= safmn2;
 118             g1 *= safmn2;
 119 /* Computing MAX */
 120             d__1 = abs(f1), d__2 = abs(g1);
 121             scale = max(d__1,d__2);
 122             if (scale >= safmx2) {
 123                 goto L10;
 124             }
 125 /* Computing 2nd power */
 126             d__1 = f1;
 127 /* Computing 2nd power */
 128             d__2 = g1;
 129             *r__ = sqrt(d__1 * d__1 + d__2 * d__2);
 130             *cs = f1 / *r__;
 131             *sn = g1 / *r__;
 132             i__1 = count;
 133             for (i__ = 1; i__ <= i__1; ++i__) {
 134                 *r__ *= safmx2;
 135 /* L20: */
 136             }
 137         } else if (scale <= safmn2) {
 138             count = 0;
 139 L30:
 140             ++count;
 141             f1 *= safmx2;
 142             g1 *= safmx2;
 143 /* Computing MAX */
 144             d__1 = abs(f1), d__2 = abs(g1);
 145             scale = max(d__1,d__2);
 146             if (scale <= safmn2) {
 147                 goto L30;
 148             }
 149 /* Computing 2nd power */
 150             d__1 = f1;
 151 /* Computing 2nd power */
 152             d__2 = g1;
 153             *r__ = sqrt(d__1 * d__1 + d__2 * d__2);
 154             *cs = f1 / *r__;
 155             *sn = g1 / *r__;
 156             i__1 = count;
 157             for (i__ = 1; i__ <= i__1; ++i__) {
 158                 *r__ *= safmn2;
 159 /* L40: */
 160             }
 161         } else {
 162 /* Computing 2nd power */
 163             d__1 = f1;
 164 /* Computing 2nd power */
 165             d__2 = g1;
 166             *r__ = sqrt(d__1 * d__1 + d__2 * d__2);
 167             *cs = f1 / *r__;
 168             *sn = g1 / *r__;
 169         }
 170         if (abs(*f) > abs(*g) && *cs < 0.) {
 171             *cs = -(*cs);
 172             *sn = -(*sn);
 173             *r__ = -(*r__);
 174         }
 175     }
 176     return 0;
 177
 178 /*     End of DLARTG */
 179
 180 } /* dlartg_ */