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
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