1 *> \brief \b DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DLASQ6 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq6.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq6.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq6.f">
21 * SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
24 * .. Scalar Arguments ..
26 * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
28 * .. Array Arguments ..
29 * DOUBLE PRECISION Z( * )
38 *> DLASQ6 computes one dqd (shift equal to zero) transform in
39 *> ping-pong form, with protection against underflow and overflow.
59 *> Z is DOUBLE PRECISION array, dimension ( 4*N )
60 *> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
67 *> PP=0 for ping, PP=1 for pong.
72 *> DMIN is DOUBLE PRECISION
73 *> Minimum value of d.
78 *> DMIN1 is DOUBLE PRECISION
79 *> Minimum value of d, excluding D( N0 ).
84 *> DMIN2 is DOUBLE PRECISION
85 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
90 *> DN is DOUBLE PRECISION
91 *> d(N0), the last value of d.
96 *> DNM1 is DOUBLE PRECISION
102 *> DNM2 is DOUBLE PRECISION
109 *> \author Univ. of Tennessee
110 *> \author Univ. of California Berkeley
111 *> \author Univ. of Colorado Denver
114 *> \date September 2012
116 *> \ingroup auxOTHERcomputational
118 * =====================================================================
119 SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
122 * -- LAPACK computational routine (version 3.4.2) --
123 * -- LAPACK is a software package provided by Univ. of Tennessee, --
124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127 * .. Scalar Arguments ..
129 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
131 * .. Array Arguments ..
132 DOUBLE PRECISION Z( * )
135 * =====================================================================
138 DOUBLE PRECISION ZERO
139 PARAMETER ( ZERO = 0.0D0 )
141 * .. Local Scalars ..
143 DOUBLE PRECISION D, EMIN, SAFMIN, TEMP
145 * .. External Function ..
146 DOUBLE PRECISION DLAMCH
149 * .. Intrinsic Functions ..
152 * .. Executable Statements ..
154 IF( ( N0-I0-1 ).LE.0 )
161 SAFMIN = DLAMCH( 'Safe minimum' )
168 DO 10 J4 = 4*I0, 4*( N0-3 ), 4
169 Z( J4-2 ) = D + Z( J4-1 )
170 IF( Z( J4-2 ).EQ.ZERO ) THEN
175 ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
176 $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
177 TEMP = Z( J4+1 ) / Z( J4-2 )
178 Z( J4 ) = Z( J4-1 )*TEMP
181 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
182 D = Z( J4+1 )*( D / Z( J4-2 ) )
184 DMIN = MIN( DMIN, D )
185 EMIN = MIN( EMIN, Z( J4 ) )
188 DO 20 J4 = 4*I0, 4*( N0-3 ), 4
189 Z( J4-3 ) = D + Z( J4 )
190 IF( Z( J4-3 ).EQ.ZERO ) THEN
195 ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
196 $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
197 TEMP = Z( J4+2 ) / Z( J4-3 )
198 Z( J4-1 ) = Z( J4 )*TEMP
201 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
202 D = Z( J4+2 )*( D / Z( J4-3 ) )
204 DMIN = MIN( DMIN, D )
205 EMIN = MIN( EMIN, Z( J4-1 ) )
209 * Unroll last two steps.
215 Z( J4-2 ) = DNM2 + Z( J4P2 )
216 IF( Z( J4-2 ).EQ.ZERO ) THEN
221 ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
222 $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
223 TEMP = Z( J4P2+2 ) / Z( J4-2 )
224 Z( J4 ) = Z( J4P2 )*TEMP
227 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
228 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
230 DMIN = MIN( DMIN, DNM1 )
235 Z( J4-2 ) = DNM1 + Z( J4P2 )
236 IF( Z( J4-2 ).EQ.ZERO ) THEN
241 ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
242 $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
243 TEMP = Z( J4P2+2 ) / Z( J4-2 )
244 Z( J4 ) = Z( J4P2 )*TEMP
247 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
248 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
250 DMIN = MIN( DMIN, DN )