1: SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION DD1,DD2,DX1,DY1
4: * ..
5: * .. Array Arguments ..
6: DOUBLE PRECISION DPARAM(5)
7: * ..
8: *
9: * Purpose
10: * =======
11: *
12: * CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
13: * THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*
14: * DY2)**T.
15: * WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
16: *
17: * DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
18: *
19: * (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
20: * H=( ) ( ) ( ) ( )
21: * (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
22: * LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
23: * RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
24: * VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
25: *
26: * THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
27: * INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
28: * OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
29: *
30: *
31: * Arguments
32: * =========
33: *
34: * DD1 (input/output) DOUBLE PRECISION
35: *
36: * DD2 (input/output) DOUBLE PRECISION
37: *
38: * DX1 (input/output) DOUBLE PRECISION
39: *
40: * DY1 (input) DOUBLE PRECISION
41: *
42: * DPARAM (input/output) DOUBLE PRECISION array, dimension 5
43: * DPARAM(1)=DFLAG
44: * DPARAM(2)=DH11
45: * DPARAM(3)=DH21
46: * DPARAM(4)=DH12
47: * DPARAM(5)=DH22
48: *
49: * =====================================================================
50: *
51: * .. Local Scalars ..
52: DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
53: $ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
54: * ..
55: * .. Intrinsic Functions ..
56: INTRINSIC DABS
57: * ..
58: * .. Data statements ..
59: *
60: DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/
61: DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/
62: * ..
63:
64: IF (DD1.LT.ZERO) THEN
65: * GO ZERO-H-D-AND-DX1..
66: DFLAG = -ONE
67: DH11 = ZERO
68: DH12 = ZERO
69: DH21 = ZERO
70: DH22 = ZERO
71: *
72: DD1 = ZERO
73: DD2 = ZERO
74: DX1 = ZERO
75: ELSE
76: * CASE-DD1-NONNEGATIVE
77: DP2 = DD2*DY1
78: IF (DP2.EQ.ZERO) THEN
79: DFLAG = -TWO
80: DPARAM(1) = DFLAG
81: RETURN
82: END IF
83: * REGULAR-CASE..
84: DP1 = DD1*DX1
85: DQ2 = DP2*DY1
86: DQ1 = DP1*DX1
87: *
88: IF (DABS(DQ1).GT.DABS(DQ2)) THEN
89: DH21 = -DY1/DX1
90: DH12 = DP2/DP1
91: *
92: DU = ONE - DH12*DH21
93: *
94: IF (DU.GT.ZERO) THEN
95: DFLAG = ZERO
96: DD1 = DD1/DU
97: DD2 = DD2/DU
98: DX1 = DX1*DU
99: END IF
100: ELSE
101:
102: IF (DQ2.LT.ZERO) THEN
103: * GO ZERO-H-D-AND-DX1..
104: DFLAG = -ONE
105: DH11 = ZERO
106: DH12 = ZERO
107: DH21 = ZERO
108: DH22 = ZERO
109: *
110: DD1 = ZERO
111: DD2 = ZERO
112: DX1 = ZERO
113: ELSE
114: DFLAG = ONE
115: DH11 = DP1/DP2
116: DH22 = DX1/DY1
117: DU = ONE + DH11*DH22
118: DTEMP = DD2/DU
119: DD2 = DD1/DU
120: DD1 = DTEMP
121: DX1 = DY1*DU
122: END IF
123: END IF
124:
125: * PROCEDURE..SCALE-CHECK
126: IF (DD1.NE.ZERO) THEN
127: DO WHILE ((DD1.LE.RGAMSQ) .OR. (DD1.GE.GAMSQ))
128: IF (DFLAG.EQ.ZERO) THEN
129: DH11 = ONE
130: DH22 = ONE
131: DFLAG = -ONE
132: ELSE
133: DH21 = -ONE
134: DH12 = ONE
135: DFLAG = -ONE
136: END IF
137: IF (DD1.LE.RGAMSQ) THEN
138: DD1 = DD1*GAM**2
139: DX1 = DX1/GAM
140: DH11 = DH11/GAM
141: DH12 = DH12/GAM
142: ELSE
143: DD1 = DD1/GAM**2
144: DX1 = DX1*GAM
145: DH11 = DH11*GAM
146: DH12 = DH12*GAM
147: END IF
148: ENDDO
149: END IF
150:
151: IF (DD2.NE.ZERO) THEN
152: DO WHILE ( (DABS(DD2).LE.RGAMSQ) .OR. (DABS(DD2).GE.GAMSQ) )
153: IF (DFLAG.EQ.ZERO) THEN
154: DH11 = ONE
155: DH22 = ONE
156: DFLAG = -ONE
157: ELSE
158: DH21 = -ONE
159: DH12 = ONE
160: DFLAG = -ONE
161: END IF
162: IF (DABS(DD2).LE.RGAMSQ) THEN
163: DD2 = DD2*GAM**2
164: DH21 = DH21/GAM
165: DH22 = DH22/GAM
166: ELSE
167: DD2 = DD2/GAM**2
168: DH21 = DH21*GAM
169: DH22 = DH22*GAM
170: END IF
171: END DO
172: END IF
173:
174: END IF
175:
176: IF (DFLAG.LT.ZERO) THEN
177: DPARAM(2) = DH11
178: DPARAM(3) = DH21
179: DPARAM(4) = DH12
180: DPARAM(5) = DH22
181: ELSE IF (DFLAG.EQ.ZERO) THEN
182: DPARAM(3) = DH21
183: DPARAM(4) = DH12
184: ELSE
185: DPARAM(2) = DH11
186: DPARAM(5) = DH22
187: END IF
188:
189: 260 CONTINUE
190: DPARAM(1) = DFLAG
191: RETURN
192: END
193:
194:
195:
196:
CVSweb interface <joel.bertrand@systella.fr>