3 /* Subroutine */ int dlarfg_(integer *n, doublereal *alpha, doublereal *x,
4 integer *incx, doublereal *tau)
6 /* System generated locals */
10 /* Builtin functions */
11 double d_sign(doublereal *, doublereal *);
16 extern doublereal dnrm2_(integer *, doublereal *, integer *);
17 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
20 extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
21 doublereal safmin, rsafmn;
24 /* -- LAPACK auxiliary routine (version 3.1) -- */
25 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
28 /* .. Scalar Arguments .. */
30 /* .. Array Arguments .. */
36 /* DLARFG generates a real elementary reflector H of order n, such */
39 /* H * ( alpha ) = ( beta ), H' * H = I. */
42 /* where alpha and beta are scalars, and x is an (n-1)-element real */
43 /* vector. H is represented in the form */
45 /* H = I - tau * ( 1 ) * ( 1 v' ) , */
48 /* where tau is a real scalar and v is a real (n-1)-element */
51 /* If the elements of x are all zero, then tau = 0 and H is taken to be */
52 /* the unit matrix. */
54 /* Otherwise 1 <= tau <= 2. */
59 /* N (input) INTEGER */
60 /* The order of the elementary reflector. */
62 /* ALPHA (input/output) DOUBLE PRECISION */
63 /* On entry, the value alpha. */
64 /* On exit, it is overwritten with the value beta. */
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. */
71 /* INCX (input) INTEGER */
72 /* The increment between elements of X. INCX > 0. */
74 /* TAU (output) DOUBLE PRECISION */
77 /* ===================================================================== */
79 /* .. Parameters .. */
81 /* .. Local Scalars .. */
83 /* .. External Functions .. */
85 /* .. Intrinsic Functions .. */
87 /* .. External Subroutines .. */
89 /* .. Executable Statements .. */
91 /* Parameter adjustments */
101 xnorm = dnrm2_(&i__1, &x[1], incx);
112 d__1 = dlapy2_(alpha, &xnorm);
113 beta = -d_sign(&d__1, alpha);
114 safmin = dlamch_("S") / dlamch_("E");
115 if (abs(beta) < safmin) {
117 /* XNORM, BETA may be inaccurate; scale X and recompute them */
119 rsafmn = 1. / safmin;
124 dscal_(&i__1, &rsafmn, &x[1], incx);
127 if (abs(beta) < safmin) {
131 /* New BETA is at most 1, at least SAFMIN */
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;
139 d__1 = 1. / (*alpha - beta);
140 dscal_(&i__1, &d__1, &x[1], incx);
142 /* If ALPHA is subnormal, it may lose relative accuracy */
146 for (j = 1; j <= i__1; ++j) {
151 *tau = (beta - *alpha) / beta;
153 d__1 = 1. / (*alpha - beta);
154 dscal_(&i__1, &d__1, &x[1], incx);