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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlarfg_(integer *n, doublereal *alpha, doublereal *x,
   4         integer *incx, doublereal *tau)
   5 {
   6     /* System generated locals */
   7     integer i__1;
   8     doublereal d__1;
   9
  10     /* Builtin functions */
  11     double d_sign(doublereal *, doublereal *);
  12
  13     /* Local variables */
  14     integer j, knt;
  15     doublereal beta;
  16     extern doublereal dnrm2_(integer *, doublereal *, integer *);
  17     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  18             integer *);
  19     doublereal xnorm;
  20     extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
  21     doublereal safmin, rsafmn;
  22
  23
  24 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  25 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  26 /*     November 2006 */
  27
  28 /*     .. Scalar Arguments .. */
  29 /*     .. */
  30 /*     .. Array Arguments .. */
  31 /*     .. */
  32
  33 /*  Purpose */
  34 /*  ======= */
  35
  36 /*  DLARFG generates a real elementary reflector H of order n, such */
  37 /*  that */
  38
  39 /*        H * ( alpha ) = ( beta ),   H' * H = I. */
  40 /*            (   x   )   (   0  ) */
  41
  42 /*  where alpha and beta are scalars, and x is an (n-1)-element real */
  43 /*  vector. H is represented in the form */
  44
  45 /*        H = I - tau * ( 1 ) * ( 1 v' ) , */
  46 /*                      ( v ) */
  47
  48 /*  where tau is a real scalar and v is a real (n-1)-element */
  49 /*  vector. */
  50
  51 /*  If the elements of x are all zero, then tau = 0 and H is taken to be */
  52 /*  the unit matrix. */
  53
  54 /*  Otherwise  1 <= tau <= 2. */
  55
  56 /*  Arguments */
  57 /*  ========= */
  58
  59 /*  N       (input) INTEGER */
  60 /*          The order of the elementary reflector. */
  61
  62 /*  ALPHA   (input/output) DOUBLE PRECISION */
  63 /*          On entry, the value alpha. */
  64 /*          On exit, it is overwritten with the value beta. */
  65
  66 /*  X       (input/output) DOUBLE PRECISION array, dimension */
  67 /*                         (1+(N-2)*abs(INCX)) */
  68 /*          On entry, the vector x. */
  69 /*          On exit, it is overwritten with the vector v. */
  70
  71 /*  INCX    (input) INTEGER */
  72 /*          The increment between elements of X. INCX > 0. */
  73
  74 /*  TAU     (output) DOUBLE PRECISION */
  75 /*          The value tau. */
  76
  77 /*  ===================================================================== */
  78
  79 /*     .. Parameters .. */
  80 /*     .. */
  81 /*     .. Local Scalars .. */
  82 /*     .. */
  83 /*     .. External Functions .. */
  84 /*     .. */
  85 /*     .. Intrinsic Functions .. */
  86 /*     .. */
  87 /*     .. External Subroutines .. */
  88 /*     .. */
  89 /*     .. Executable Statements .. */
  90
  91     /* Parameter adjustments */
  92     --x;
  93
  94     /* Function Body */
  95     if (*n <= 1) {
  96         *tau = 0.;
  97         return 0;
  98     }
  99
 100     i__1 = *n - 1;
 101     xnorm = dnrm2_(&i__1, &x[1], incx);
 102
 103     if (xnorm == 0.) {
 104
 105 /*        H  =  I */
 106
 107         *tau = 0.;
 108     } else {
 109
 110 /*        general case */
 111
 112         d__1 = dlapy2_(alpha, &xnorm);
 113         beta = -d_sign(&d__1, alpha);
 114         safmin = dlamch_("S") / dlamch_("E");
 115         if (abs(beta) < safmin) {
 116
 117 /*           XNORM, BETA may be inaccurate; scale X and recompute them */
 118
 119             rsafmn = 1. / safmin;
 120             knt = 0;
 121 L10:
 122             ++knt;
 123             i__1 = *n - 1;
 124             dscal_(&i__1, &rsafmn, &x[1], incx);
 125             beta *= rsafmn;
 126             *alpha *= rsafmn;
 127             if (abs(beta) < safmin) {
 128                 goto L10;
 129             }
 130
 131 /*           New BETA is at most 1, at least SAFMIN */
 132
 133             i__1 = *n - 1;
 134             xnorm = dnrm2_(&i__1, &x[1], incx);
 135             d__1 = dlapy2_(alpha, &xnorm);
 136             beta = -d_sign(&d__1, alpha);
 137             *tau = (beta - *alpha) / beta;
 138             i__1 = *n - 1;
 139             d__1 = 1. / (*alpha - beta);
 140             dscal_(&i__1, &d__1, &x[1], incx);
 141
 142 /*           If ALPHA is subnormal, it may lose relative accuracy */
 143
 144             *alpha = beta;
 145             i__1 = knt;
 146             for (j = 1; j <= i__1; ++j) {
 147                 *alpha *= safmin;
 148 /* L20: */
 149             }
 150         } else {
 151             *tau = (beta - *alpha) / beta;
 152             i__1 = *n - 1;
 153             d__1 = 1. / (*alpha - beta);
 154             dscal_(&i__1, &d__1, &x[1], incx);
 155             *alpha = beta;
 156         }
 157     }
 158
 159     return 0;
 160
 161 /*     End of DLARFG */
 162
 163 } /* dlarfg_ */