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