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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
* .. Scalar Arguments ..
DOUBLE PRECISION DD1,DD2,DX1,DY1
* ..
* .. Array Arguments ..
DOUBLE PRECISION DPARAM(5)
* ..
*
* Purpose
* =======
*
* CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
* THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*
* DY2)**T.
* WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*
* DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
*
* (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
* H=( ) ( ) ( ) ( )
* (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
* LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
* RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
* VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
*
* THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
* INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
* OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
*
*
* Arguments
* =========
*
* DD1 (input/output) DOUBLE PRECISION
*
* DD2 (input/output) DOUBLE PRECISION
*
* DX1 (input/output) DOUBLE PRECISION
*
* DY1 (input) DOUBLE PRECISION
*
* DPARAM (input/output) DOUBLE PRECISION array, dimension 5
* DPARAM(1)=DFLAG
* DPARAM(2)=DH11
* DPARAM(3)=DH21
* DPARAM(4)=DH12
* DPARAM(5)=DH22
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
+ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
INTEGER IGO
* ..
* .. Intrinsic Functions ..
INTRINSIC DABS
* ..
* .. Data statements ..
*
DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/
DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/
* ..
IF (.NOT.DD1.LT.ZERO) GO TO 10
* GO ZERO-H-D-AND-DX1..
GO TO 60
10 CONTINUE
* CASE-DD1-NONNEGATIVE
DP2 = DD2*DY1
IF (.NOT.DP2.EQ.ZERO) GO TO 20
DFLAG = -TWO
GO TO 260
* REGULAR-CASE..
20 CONTINUE
DP1 = DD1*DX1
DQ2 = DP2*DY1
DQ1 = DP1*DX1
*
IF (.NOT.DABS(DQ1).GT.DABS(DQ2)) GO TO 40
DH21 = -DY1/DX1
DH12 = DP2/DP1
*
DU = ONE - DH12*DH21
*
IF (.NOT.DU.LE.ZERO) GO TO 30
* GO ZERO-H-D-AND-DX1..
GO TO 60
30 CONTINUE
DFLAG = ZERO
DD1 = DD1/DU
DD2 = DD2/DU
DX1 = DX1*DU
* GO SCALE-CHECK..
GO TO 100
40 CONTINUE
IF (.NOT.DQ2.LT.ZERO) GO TO 50
* GO ZERO-H-D-AND-DX1..
GO TO 60
50 CONTINUE
DFLAG = ONE
DH11 = DP1/DP2
DH22 = DX1/DY1
DU = ONE + DH11*DH22
DTEMP = DD2/DU
DD2 = DD1/DU
DD1 = DTEMP
DX1 = DY1*DU
* GO SCALE-CHECK
GO TO 100
* PROCEDURE..ZERO-H-D-AND-DX1..
60 CONTINUE
DFLAG = -ONE
DH11 = ZERO
DH12 = ZERO
DH21 = ZERO
DH22 = ZERO
*
DD1 = ZERO
DD2 = ZERO
DX1 = ZERO
* RETURN..
GO TO 220
* PROCEDURE..FIX-H..
70 CONTINUE
IF (.NOT.DFLAG.GE.ZERO) GO TO 90
*
IF (.NOT.DFLAG.EQ.ZERO) GO TO 80
DH11 = ONE
DH22 = ONE
DFLAG = -ONE
GO TO 90
80 CONTINUE
DH21 = -ONE
DH12 = ONE
DFLAG = -ONE
90 CONTINUE
GO TO IGO(120,150,180,210)
* PROCEDURE..SCALE-CHECK
100 CONTINUE
110 CONTINUE
IF (.NOT.DD1.LE.RGAMSQ) GO TO 130
IF (DD1.EQ.ZERO) GO TO 160
ASSIGN 120 TO IGO
* FIX-H..
GO TO 70
120 CONTINUE
DD1 = DD1*GAM**2
DX1 = DX1/GAM
DH11 = DH11/GAM
DH12 = DH12/GAM
GO TO 110
130 CONTINUE
140 CONTINUE
IF (.NOT.DD1.GE.GAMSQ) GO TO 160
ASSIGN 150 TO IGO
* FIX-H..
GO TO 70
150 CONTINUE
DD1 = DD1/GAM**2
DX1 = DX1*GAM
DH11 = DH11*GAM
DH12 = DH12*GAM
GO TO 140
160 CONTINUE
170 CONTINUE
IF (.NOT.DABS(DD2).LE.RGAMSQ) GO TO 190
IF (DD2.EQ.ZERO) GO TO 220
ASSIGN 180 TO IGO
* FIX-H..
GO TO 70
180 CONTINUE
DD2 = DD2*GAM**2
DH21 = DH21/GAM
DH22 = DH22/GAM
GO TO 170
190 CONTINUE
200 CONTINUE
IF (.NOT.DABS(DD2).GE.GAMSQ) GO TO 220
ASSIGN 210 TO IGO
* FIX-H..
GO TO 70
210 CONTINUE
DD2 = DD2/GAM**2
DH21 = DH21*GAM
DH22 = DH22*GAM
GO TO 200
220 CONTINUE
IF (DFLAG) 250,230,240
230 CONTINUE
DPARAM(3) = DH21
DPARAM(4) = DH12
GO TO 260
240 CONTINUE
DPARAM(2) = DH11
DPARAM(5) = DH22
GO TO 260
250 CONTINUE
DPARAM(2) = DH11
DPARAM(3) = DH21
DPARAM(4) = DH12
DPARAM(5) = DH22
260 CONTINUE
DPARAM(1) = DFLAG
RETURN
END
|
