![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlasq5.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
2: $ DNM1, DNM2, IEEE )
3: *
4: * -- LAPACK routine (version 3.2) --
5: *
6: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
7: * -- Laboratory and Beresford Parlett of the Univ. of California at --
8: * -- Berkeley --
9: * -- November 2008 --
10: *
11: * -- LAPACK is a software package provided by Univ. of Tennessee, --
12: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13: *
14: * .. Scalar Arguments ..
15: LOGICAL IEEE
16: INTEGER I0, N0, PP
17: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
18: * ..
19: * .. Array Arguments ..
20: DOUBLE PRECISION Z( * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * DLASQ5 computes one dqds transform in ping-pong form, one
27: * version for IEEE machines another for non IEEE machines.
28: *
29: * Arguments
30: * =========
31: *
32: * I0 (input) INTEGER
33: * First index.
34: *
35: * N0 (input) INTEGER
36: * Last index.
37: *
38: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
39: * Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
40: * an extra argument.
41: *
42: * PP (input) INTEGER
43: * PP=0 for ping, PP=1 for pong.
44: *
45: * TAU (input) DOUBLE PRECISION
46: * This is the shift.
47: *
48: * DMIN (output) DOUBLE PRECISION
49: * Minimum value of d.
50: *
51: * DMIN1 (output) DOUBLE PRECISION
52: * Minimum value of d, excluding D( N0 ).
53: *
54: * DMIN2 (output) DOUBLE PRECISION
55: * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
56: *
57: * DN (output) DOUBLE PRECISION
58: * d(N0), the last value of d.
59: *
60: * DNM1 (output) DOUBLE PRECISION
61: * d(N0-1).
62: *
63: * DNM2 (output) DOUBLE PRECISION
64: * d(N0-2).
65: *
66: * IEEE (input) LOGICAL
67: * Flag for IEEE or non IEEE arithmetic.
68: *
69: * =====================================================================
70: *
71: * .. Parameter ..
72: DOUBLE PRECISION ZERO
73: PARAMETER ( ZERO = 0.0D0 )
74: * ..
75: * .. Local Scalars ..
76: INTEGER J4, J4P2
77: DOUBLE PRECISION D, EMIN, TEMP
78: * ..
79: * .. Intrinsic Functions ..
80: INTRINSIC MIN
81: * ..
82: * .. Executable Statements ..
83: *
84: IF( ( N0-I0-1 ).LE.0 )
85: $ RETURN
86: *
87: J4 = 4*I0 + PP - 3
88: EMIN = Z( J4+4 )
89: D = Z( J4 ) - TAU
90: DMIN = D
91: DMIN1 = -Z( J4 )
92: *
93: IF( IEEE ) THEN
94: *
95: * Code for IEEE arithmetic.
96: *
97: IF( PP.EQ.0 ) THEN
98: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
99: Z( J4-2 ) = D + Z( J4-1 )
100: TEMP = Z( J4+1 ) / Z( J4-2 )
101: D = D*TEMP - TAU
102: DMIN = MIN( DMIN, D )
103: Z( J4 ) = Z( J4-1 )*TEMP
104: EMIN = MIN( Z( J4 ), EMIN )
105: 10 CONTINUE
106: ELSE
107: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
108: Z( J4-3 ) = D + Z( J4 )
109: TEMP = Z( J4+2 ) / Z( J4-3 )
110: D = D*TEMP - TAU
111: DMIN = MIN( DMIN, D )
112: Z( J4-1 ) = Z( J4 )*TEMP
113: EMIN = MIN( Z( J4-1 ), EMIN )
114: 20 CONTINUE
115: END IF
116: *
117: * Unroll last two steps.
118: *
119: DNM2 = D
120: DMIN2 = DMIN
121: J4 = 4*( N0-2 ) - PP
122: J4P2 = J4 + 2*PP - 1
123: Z( J4-2 ) = DNM2 + Z( J4P2 )
124: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
125: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
126: DMIN = MIN( DMIN, DNM1 )
127: *
128: DMIN1 = DMIN
129: J4 = J4 + 4
130: J4P2 = J4 + 2*PP - 1
131: Z( J4-2 ) = DNM1 + Z( J4P2 )
132: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
133: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
134: DMIN = MIN( DMIN, DN )
135: *
136: ELSE
137: *
138: * Code for non IEEE arithmetic.
139: *
140: IF( PP.EQ.0 ) THEN
141: DO 30 J4 = 4*I0, 4*( N0-3 ), 4
142: Z( J4-2 ) = D + Z( J4-1 )
143: IF( D.LT.ZERO ) THEN
144: RETURN
145: ELSE
146: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
147: D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
148: END IF
149: DMIN = MIN( DMIN, D )
150: EMIN = MIN( EMIN, Z( J4 ) )
151: 30 CONTINUE
152: ELSE
153: DO 40 J4 = 4*I0, 4*( N0-3 ), 4
154: Z( J4-3 ) = D + Z( J4 )
155: IF( D.LT.ZERO ) THEN
156: RETURN
157: ELSE
158: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
159: D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
160: END IF
161: DMIN = MIN( DMIN, D )
162: EMIN = MIN( EMIN, Z( J4-1 ) )
163: 40 CONTINUE
164: END IF
165: *
166: * Unroll last two steps.
167: *
168: DNM2 = D
169: DMIN2 = DMIN
170: J4 = 4*( N0-2 ) - PP
171: J4P2 = J4 + 2*PP - 1
172: Z( J4-2 ) = DNM2 + Z( J4P2 )
173: IF( DNM2.LT.ZERO ) THEN
174: RETURN
175: ELSE
176: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
177: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
178: END IF
179: DMIN = MIN( DMIN, DNM1 )
180: *
181: DMIN1 = DMIN
182: J4 = J4 + 4
183: J4P2 = J4 + 2*PP - 1
184: Z( J4-2 ) = DNM1 + Z( J4P2 )
185: IF( DNM1.LT.ZERO ) THEN
186: RETURN
187: ELSE
188: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
189: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
190: END IF
191: DMIN = MIN( DMIN, DN )
192: *
193: END IF
194: *
195: Z( J4+2 ) = DN
196: Z( 4*N0-PP ) = EMIN
197: RETURN
198: *
199: * End of DLASQ5
200: *
201: END