1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
SUBROUTINE DLARTG( F, G, CS, SN, R )
*
* -- LAPACK auxiliary routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
DOUBLE PRECISION CS, F, G, R, SN
* ..
*
* Purpose
* =======
*
* DLARTG generate a plane rotation so that
*
* [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1.
* [ -SN CS ] [ G ] [ 0 ]
*
* This is a slower, more accurate version of the BLAS1 routine DROTG,
* with the following other differences:
* F and G are unchanged on return.
* If G=0, then CS=1 and SN=0.
* If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
* floating point operations (saves work in DBDSQR when
* there are zeros on the diagonal).
*
* If F exceeds G in magnitude, CS will be positive.
*
* Arguments
* =========
*
* F (input) DOUBLE PRECISION
* The first component of vector to be rotated.
*
* G (input) DOUBLE PRECISION
* The second component of vector to be rotated.
*
* CS (output) DOUBLE PRECISION
* The cosine of the rotation.
*
* SN (output) DOUBLE PRECISION
* The sine of the rotation.
*
* R (output) DOUBLE PRECISION
* The nonzero component of the rotated vector.
*
* This version has a few statements commented out for thread safety
* (machine parameters are computed on each entry). 10 feb 03, SJH.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D0 )
DOUBLE PRECISION TWO
PARAMETER ( TWO = 2.0D0 )
* ..
* .. Local Scalars ..
* LOGICAL FIRST
INTEGER COUNT, I
DOUBLE PRECISION EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, INT, LOG, MAX, SQRT
* ..
* .. Save statement ..
* SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
* ..
* .. Data statements ..
* DATA FIRST / .TRUE. /
* ..
* .. Executable Statements ..
*
* IF( FIRST ) THEN
SAFMIN = DLAMCH( 'S' )
EPS = DLAMCH( 'E' )
SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
$ LOG( DLAMCH( 'B' ) ) / TWO )
SAFMX2 = ONE / SAFMN2
* FIRST = .FALSE.
* END IF
IF( G.EQ.ZERO ) THEN
CS = ONE
SN = ZERO
R = F
ELSE IF( F.EQ.ZERO ) THEN
CS = ZERO
SN = ONE
R = G
ELSE
F1 = F
G1 = G
SCALE = MAX( ABS( F1 ), ABS( G1 ) )
IF( SCALE.GE.SAFMX2 ) THEN
COUNT = 0
10 CONTINUE
COUNT = COUNT + 1
F1 = F1*SAFMN2
G1 = G1*SAFMN2
SCALE = MAX( ABS( F1 ), ABS( G1 ) )
IF( SCALE.GE.SAFMX2 )
$ GO TO 10
R = SQRT( F1**2+G1**2 )
CS = F1 / R
SN = G1 / R
DO 20 I = 1, COUNT
R = R*SAFMX2
20 CONTINUE
ELSE IF( SCALE.LE.SAFMN2 ) THEN
COUNT = 0
30 CONTINUE
COUNT = COUNT + 1
F1 = F1*SAFMX2
G1 = G1*SAFMX2
SCALE = MAX( ABS( F1 ), ABS( G1 ) )
IF( SCALE.LE.SAFMN2 )
$ GO TO 30
R = SQRT( F1**2+G1**2 )
CS = F1 / R
SN = G1 / R
DO 40 I = 1, COUNT
R = R*SAFMN2
40 CONTINUE
ELSE
R = SQRT( F1**2+G1**2 )
CS = F1 / R
SN = G1 / R
END IF
IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN
CS = -CS
SN = -SN
R = -R
END IF
END IF
RETURN
*
* End of DLARTG
*
END
|
