Mercurial > octave-nkf

comparison libcruft/lapack/clarfg.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 CLARFG( 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 COMPLEX ALPHA, TAU
10 * ..
11 * .. Array Arguments ..
12 COMPLEX X( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * CLARFG generates a complex 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, with beta real, and x is an
25 * (n-1)-element complex vector. H is represented in the form
26 *
27 * H = I - tau * ( 1 ) * ( 1 v' ) ,
28 * ( v )
29 *
30 * where tau is a complex scalar and v is a complex (n-1)-element
31 * vector. Note that H is not hermitian.
32 *
33 * If the elements of x are all zero and alpha is real, then tau = 0
34 * and H is taken to be the unit matrix.
35 *
36 * Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 .
37 *
38 * Arguments
39 * =========
40 *
41 * N (input) INTEGER
42 * The order of the elementary reflector.
43 *
44 * ALPHA (input/output) COMPLEX
45 * On entry, the value alpha.
46 * On exit, it is overwritten with the value beta.
47 *
48 * X (input/output) COMPLEX 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) COMPLEX
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 ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
68 * ..
69 * .. External Functions ..
70 REAL SCNRM2, SLAMCH, SLAPY3
71 COMPLEX CLADIV
72 EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV
73 * ..
74 * .. Intrinsic Functions ..
75 INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL CSCAL, CSSCAL
79 * ..
80 * .. Executable Statements ..
81 *
82 IF( N.LE.0 ) THEN
83 TAU = ZERO
84 RETURN
85 END IF
86 *
87 XNORM = SCNRM2( N-1, X, INCX )
88 ALPHR = REAL( ALPHA )
89 ALPHI = AIMAG( ALPHA )
90 *
91 IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
92 *
93 * H = I
94 *
95 TAU = ZERO
96 ELSE
97 *
98 * general case
99 *
100 BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
101 SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
102 RSAFMN = ONE / SAFMIN
103 *
104 IF( ABS( BETA ).LT.SAFMIN ) THEN
105 *
106 * XNORM, BETA may be inaccurate; scale X and recompute them
107 *
108 KNT = 0
109 10 CONTINUE
110 KNT = KNT + 1
111 CALL CSSCAL( N-1, RSAFMN, X, INCX )
112 BETA = BETA*RSAFMN
113 ALPHI = ALPHI*RSAFMN
114 ALPHR = ALPHR*RSAFMN
115 IF( ABS( BETA ).LT.SAFMIN )
116 $ GO TO 10
117 *
118 * New BETA is at most 1, at least SAFMIN
119 *
120 XNORM = SCNRM2( N-1, X, INCX )
121 ALPHA = CMPLX( ALPHR, ALPHI )
122 BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
123 TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
124 ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
125 CALL CSCAL( N-1, ALPHA, X, INCX )
126 *
127 * If ALPHA is subnormal, it may lose relative accuracy
128 *
129 ALPHA = BETA
130 DO 20 J = 1, KNT
131 ALPHA = ALPHA*SAFMIN
132 20 CONTINUE
133 ELSE
134 TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
135 ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
136 CALL CSCAL( N-1, ALPHA, X, INCX )
137 ALPHA = BETA
138 END IF
139 END IF
140 *
141 RETURN
142 *
143 * End of CLARFG
144 *
145 END