Update to 2.0.0 tree from current Fremantle build

[opencv] / 3rdparty / lapack / slarfg.c

diff --git a/3rdparty/lapack/slarfg.c b/3rdparty/lapack/slarfg.c
new file mode 100644 (file)
index 0000000..bd99e17
--- /dev/null
+++ b/3rdparty/lapack/slarfg.c
@@ -0,0 +1,162 @@
+#include "clapack.h"
+
+/* Subroutine */ int slarfg_(integer *n, real *alpha, real *x, integer *incx, 
+       real *tau)
+{
+    /* System generated locals */
+    integer i__1;
+    real r__1;
+
+    /* Builtin functions */
+    double r_sign(real *, real *);
+
+    /* Local variables */
+    integer j, knt;
+    real beta;
+    extern doublereal snrm2_(integer *, real *, integer *);
+    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
+    real xnorm;
+    extern doublereal slapy2_(real *, real *), slamch_(char *);
+    real safmin, rsafmn;
+
+
+/*  -- LAPACK auxiliary routine (version 3.1) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SLARFG generates a real elementary reflector H of order n, such */
+/*  that */
+
+/*        H * ( alpha ) = ( beta ),   H' * H = I. */
+/*            (   x   )   (   0  ) */
+
+/*  where alpha and beta are scalars, and x is an (n-1)-element real */
+/*  vector. H is represented in the form */
+
+/*        H = I - tau * ( 1 ) * ( 1 v' ) , */
+/*                      ( v ) */
+
+/*  where tau is a real scalar and v is a real (n-1)-element */
+/*  vector. */
+
+/*  If the elements of x are all zero, then tau = 0 and H is taken to be */
+/*  the unit matrix. */
+
+/*  Otherwise  1 <= tau <= 2. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  N       (input) INTEGER */
+/*          The order of the elementary reflector. */
+
+/*  ALPHA   (input/output) REAL */
+/*          On entry, the value alpha. */
+/*          On exit, it is overwritten with the value beta. */
+
+/*  X       (input/output) REAL array, dimension */
+/*                         (1+(N-2)*abs(INCX)) */
+/*          On entry, the vector x. */
+/*          On exit, it is overwritten with the vector v. */
+
+/*  INCX    (input) INTEGER */
+/*          The increment between elements of X. INCX > 0. */
+
+/*  TAU     (output) REAL */
+/*          The value tau. */
+
+/*  ===================================================================== */
+
+/*     .. Parameters .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+    /* Parameter adjustments */
+    --x;
+
+    /* Function Body */
+    if (*n <= 1) {
+       *tau = 0.f;
+       return 0;
+    }
+
+    i__1 = *n - 1;
+    xnorm = snrm2_(&i__1, &x[1], incx);
+
+    if (xnorm == 0.f) {
+
+/*        H  =  I */
+
+       *tau = 0.f;
+    } else {
+
+/*        general case */
+
+       r__1 = slapy2_(alpha, &xnorm);
+       beta = -r_sign(&r__1, alpha);
+       safmin = slamch_("S") / slamch_("E");
+       if (dabs(beta) < safmin) {
+
+/*           XNORM, BETA may be inaccurate; scale X and recompute them */
+
+           rsafmn = 1.f / safmin;
+           knt = 0;
+L10:
+           ++knt;
+           i__1 = *n - 1;
+           sscal_(&i__1, &rsafmn, &x[1], incx);
+           beta *= rsafmn;
+           *alpha *= rsafmn;
+           if (dabs(beta) < safmin) {
+               goto L10;
+           }
+
+/*           New BETA is at most 1, at least SAFMIN */
+
+           i__1 = *n - 1;
+           xnorm = snrm2_(&i__1, &x[1], incx);
+           r__1 = slapy2_(alpha, &xnorm);
+           beta = -r_sign(&r__1, alpha);
+           *tau = (beta - *alpha) / beta;
+           i__1 = *n - 1;
+           r__1 = 1.f / (*alpha - beta);
+           sscal_(&i__1, &r__1, &x[1], incx);
+
+/*           If ALPHA is subnormal, it may lose relative accuracy */
+
+           *alpha = beta;
+           i__1 = knt;
+           for (j = 1; j <= i__1; ++j) {
+               *alpha *= safmin;
+/* L20: */
+           }
+       } else {
+           *tau = (beta - *alpha) / beta;
+           i__1 = *n - 1;
+           r__1 = 1.f / (*alpha - beta);
+           sscal_(&i__1, &r__1, &x[1], incx);
+           *alpha = beta;
+       }
+    }
+
+    return 0;
+
+/*     End of SLARFG */
+
+} /* slarfg_ */