![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq6.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 dlasq6.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
2: $ DNM1, DNM2 )
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: INTEGER I0, N0, PP
16: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION Z( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DLASQ6 computes one dqd (shift equal to zero) transform in
26: * ping-pong form, with protection against underflow and overflow.
27: *
28: * Arguments
29: * =========
30: *
31: * I0 (input) INTEGER
32: * First index.
33: *
34: * N0 (input) INTEGER
35: * Last index.
36: *
37: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
38: * Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
39: * an extra argument.
40: *
41: * PP (input) INTEGER
42: * PP=0 for ping, PP=1 for pong.
43: *
44: * DMIN (output) DOUBLE PRECISION
45: * Minimum value of d.
46: *
47: * DMIN1 (output) DOUBLE PRECISION
48: * Minimum value of d, excluding D( N0 ).
49: *
50: * DMIN2 (output) DOUBLE PRECISION
51: * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
52: *
53: * DN (output) DOUBLE PRECISION
54: * d(N0), the last value of d.
55: *
56: * DNM1 (output) DOUBLE PRECISION
57: * d(N0-1).
58: *
59: * DNM2 (output) DOUBLE PRECISION
60: * d(N0-2).
61: *
62: * =====================================================================
63: *
64: * .. Parameter ..
65: DOUBLE PRECISION ZERO
66: PARAMETER ( ZERO = 0.0D0 )
67: * ..
68: * .. Local Scalars ..
69: INTEGER J4, J4P2
70: DOUBLE PRECISION D, EMIN, SAFMIN, TEMP
71: * ..
72: * .. External Function ..
73: DOUBLE PRECISION DLAMCH
74: EXTERNAL DLAMCH
75: * ..
76: * .. Intrinsic Functions ..
77: INTRINSIC MIN
78: * ..
79: * .. Executable Statements ..
80: *
81: IF( ( N0-I0-1 ).LE.0 )
82: $ RETURN
83: *
84: SAFMIN = DLAMCH( 'Safe minimum' )
85: J4 = 4*I0 + PP - 3
86: EMIN = Z( J4+4 )
87: D = Z( J4 )
88: DMIN = D
89: *
90: IF( PP.EQ.0 ) THEN
91: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
92: Z( J4-2 ) = D + Z( J4-1 )
93: IF( Z( J4-2 ).EQ.ZERO ) THEN
94: Z( J4 ) = ZERO
95: D = Z( J4+1 )
96: DMIN = D
97: EMIN = ZERO
98: ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
99: $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
100: TEMP = Z( J4+1 ) / Z( J4-2 )
101: Z( J4 ) = Z( J4-1 )*TEMP
102: D = D*TEMP
103: ELSE
104: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
105: D = Z( J4+1 )*( D / Z( J4-2 ) )
106: END IF
107: DMIN = MIN( DMIN, D )
108: EMIN = MIN( EMIN, Z( J4 ) )
109: 10 CONTINUE
110: ELSE
111: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
112: Z( J4-3 ) = D + Z( J4 )
113: IF( Z( J4-3 ).EQ.ZERO ) THEN
114: Z( J4-1 ) = ZERO
115: D = Z( J4+2 )
116: DMIN = D
117: EMIN = ZERO
118: ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
119: $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
120: TEMP = Z( J4+2 ) / Z( J4-3 )
121: Z( J4-1 ) = Z( J4 )*TEMP
122: D = D*TEMP
123: ELSE
124: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
125: D = Z( J4+2 )*( D / Z( J4-3 ) )
126: END IF
127: DMIN = MIN( DMIN, D )
128: EMIN = MIN( EMIN, Z( J4-1 ) )
129: 20 CONTINUE
130: END IF
131: *
132: * Unroll last two steps.
133: *
134: DNM2 = D
135: DMIN2 = DMIN
136: J4 = 4*( N0-2 ) - PP
137: J4P2 = J4 + 2*PP - 1
138: Z( J4-2 ) = DNM2 + Z( J4P2 )
139: IF( Z( J4-2 ).EQ.ZERO ) THEN
140: Z( J4 ) = ZERO
141: DNM1 = Z( J4P2+2 )
142: DMIN = DNM1
143: EMIN = ZERO
144: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
145: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
146: TEMP = Z( J4P2+2 ) / Z( J4-2 )
147: Z( J4 ) = Z( J4P2 )*TEMP
148: DNM1 = DNM2*TEMP
149: ELSE
150: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
151: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
152: END IF
153: DMIN = MIN( DMIN, DNM1 )
154: *
155: DMIN1 = DMIN
156: J4 = J4 + 4
157: J4P2 = J4 + 2*PP - 1
158: Z( J4-2 ) = DNM1 + Z( J4P2 )
159: IF( Z( J4-2 ).EQ.ZERO ) THEN
160: Z( J4 ) = ZERO
161: DN = Z( J4P2+2 )
162: DMIN = DN
163: EMIN = ZERO
164: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
165: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
166: TEMP = Z( J4P2+2 ) / Z( J4-2 )
167: Z( J4 ) = Z( J4P2 )*TEMP
168: DN = DNM1*TEMP
169: ELSE
170: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
171: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
172: END IF
173: DMIN = MIN( DMIN, DN )
174: *
175: Z( J4+2 ) = DN
176: Z( 4*N0-PP ) = EMIN
177: RETURN
178: *
179: * End of DLASQ6
180: *
181: END