1: SUBROUTINE DLARTGS( X, Y, SIGMA, CS, SN )
2: IMPLICIT NONE
3: *
4: * -- LAPACK routine (version 3.3.0) --
5: *
6: * -- Contributed by Brian Sutton of the Randolph-Macon College --
7: * -- November 2010
8: *
9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11: *
12: * .. Scalar Arguments ..
13: DOUBLE PRECISION CS, SIGMA, SN, X, Y
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DLARTGS generates a plane rotation designed to introduce a bulge in
20: * Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD
21: * problem. X and Y are the top-row entries, and SIGMA is the shift.
22: * The computed CS and SN define a plane rotation satisfying
23: *
24: * [ CS SN ] . [ X^2 - SIGMA ] = [ R ],
25: * [ -SN CS ] [ X * Y ] [ 0 ]
26: *
27: * with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the
28: * rotation is by PI/2.
29: *
30: * Arguments
31: * =========
32: *
33: * X (input) DOUBLE PRECISION
34: * The (1,1) entry of an upper bidiagonal matrix.
35: *
36: * Y (input) DOUBLE PRECISION
37: * The (1,2) entry of an upper bidiagonal matrix.
38: *
39: * SIGMA (input) DOUBLE PRECISION
40: * The shift.
41: *
42: * CS (output) DOUBLE PRECISION
43: * The cosine of the rotation.
44: *
45: * SN (output) DOUBLE PRECISION
46: * The sine of the rotation.
47: *
48: * ===================================================================
49: *
50: * .. Parameters ..
51: DOUBLE PRECISION NEGONE, ONE, ZERO
52: PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
53: * ..
54: * .. Local Scalars ..
55: DOUBLE PRECISION R, S, THRESH, W, Z
56: * ..
57: * .. External Functions ..
58: DOUBLE PRECISION DLAMCH
59: EXTERNAL DLAMCH
60: * .. Executable Statements ..
61: *
62: THRESH = DLAMCH('E')
63: *
64: * Compute the first column of B'*B - SIGMA^2*I, up to a scale
65: * factor.
66: *
67: IF( (SIGMA .EQ. ZERO .AND. ABS(X) .LT. THRESH) .OR.
68: $ (ABS(X) .EQ. SIGMA .AND. Y .EQ. ZERO) ) THEN
69: Z = ZERO
70: W = ZERO
71: ELSE IF( SIGMA .EQ. ZERO ) THEN
72: IF( X .GE. ZERO ) THEN
73: Z = X
74: W = Y
75: ELSE
76: Z = -X
77: W = -Y
78: END IF
79: ELSE IF( ABS(X) .LT. THRESH ) THEN
80: Z = -SIGMA*SIGMA
81: W = ZERO
82: ELSE
83: IF( X .GE. ZERO ) THEN
84: S = ONE
85: ELSE
86: S = NEGONE
87: END IF
88: Z = S * (ABS(X)-SIGMA) * (S+SIGMA/X)
89: W = S * Y
90: END IF
91: *
92: * Generate the rotation.
93: * CALL DLARTGP( Z, W, CS, SN, R ) might seem more natural;
94: * reordering the arguments ensures that if Z = 0 then the rotation
95: * is by PI/2.
96: *
97: CALL DLARTGP( W, Z, SN, CS, R )
98: *
99: RETURN
100: *
101: * End DLARTGS
102: *
103: END
104:
CVSweb interface <joel.bertrand@systella.fr>