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

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