Mercurial > octave-nkf

comparison libcruft/lapack/slarfg.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INCX, N
9 REAL ALPHA, TAU
10 * ..
11 * .. Array Arguments ..
12 REAL X( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SLARFG generates a real elementary reflector H of order n, such
19 * that
20 *
21 * H * ( alpha ) = ( beta ), H' * H = I.
22 * ( x ) ( 0 )
23 *
24 * where alpha and beta are scalars, and x is an (n-1)-element real
25 * vector. H is represented in the form
26 *
27 * H = I - tau * ( 1 ) * ( 1 v' ) ,
28 * ( v )
29 *
30 * where tau is a real scalar and v is a real (n-1)-element
31 * vector.
32 *
33 * If the elements of x are all zero, then tau = 0 and H is taken to be
34 * the unit matrix.
35 *
36 * Otherwise 1 <= tau <= 2.
37 *
38 * Arguments
39 * =========
40 *
41 * N (input) INTEGER
42 * The order of the elementary reflector.
43 *
44 * ALPHA (input/output) REAL
45 * On entry, the value alpha.
46 * On exit, it is overwritten with the value beta.
47 *
48 * X (input/output) REAL array, dimension
49 * (1+(N-2)*abs(INCX))
50 * On entry, the vector x.
51 * On exit, it is overwritten with the vector v.
52 *
53 * INCX (input) INTEGER
54 * The increment between elements of X. INCX > 0.
55 *
56 * TAU (output) REAL
57 * The value tau.
58 *
59 * =====================================================================
60 *
61 * .. Parameters ..
62 REAL ONE, ZERO
63 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER J, KNT
67 REAL BETA, RSAFMN, SAFMIN, XNORM
68 * ..
69 * .. External Functions ..
70 REAL SLAMCH, SLAPY2, SNRM2
71 EXTERNAL SLAMCH, SLAPY2, SNRM2
72 * ..
73 * .. Intrinsic Functions ..
74 INTRINSIC ABS, SIGN
75 * ..
76 * .. External Subroutines ..
77 EXTERNAL SSCAL
78 * ..
79 * .. Executable Statements ..
80 *
81 IF( N.LE.1 ) THEN
82 TAU = ZERO
83 RETURN
84 END IF
85 *
86 XNORM = SNRM2( N-1, X, INCX )
87 *
88 IF( XNORM.EQ.ZERO ) THEN
89 *
90 * H = I
91 *
92 TAU = ZERO
93 ELSE
94 *
95 * general case
96 *
97 BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
98 SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
99 IF( ABS( BETA ).LT.SAFMIN ) THEN
100 *
101 * XNORM, BETA may be inaccurate; scale X and recompute them
102 *
103 RSAFMN = ONE / SAFMIN
104 KNT = 0
105 10 CONTINUE
106 KNT = KNT + 1
107 CALL SSCAL( N-1, RSAFMN, X, INCX )
108 BETA = BETA*RSAFMN
109 ALPHA = ALPHA*RSAFMN
110 IF( ABS( BETA ).LT.SAFMIN )
111 $ GO TO 10
112 *
113 * New BETA is at most 1, at least SAFMIN
114 *
115 XNORM = SNRM2( N-1, X, INCX )
116 BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
117 TAU = ( BETA-ALPHA ) / BETA
118 CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
119 *
120 * If ALPHA is subnormal, it may lose relative accuracy
121 *
122 ALPHA = BETA
123 DO 20 J = 1, KNT
124 ALPHA = ALPHA*SAFMIN
125 20 CONTINUE
126 ELSE
127 TAU = ( BETA-ALPHA ) / BETA
128 CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
129 ALPHA = BETA
130 END IF
131 END IF
132 *
133 RETURN
134 *
135 * End of SLARFG
136 *
137 END