LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slasq3.f
Go to the documentation of this file.
1*> \brief \b SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SLASQ3 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq3.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq3.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq3.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
20* ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
21* DN2, G, TAU )
22*
23* .. Scalar Arguments ..
24* LOGICAL IEEE
25* INTEGER I0, ITER, N0, NDIV, NFAIL, PP
26* REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
27* $ QMAX, SIGMA, TAU
28* ..
29* .. Array Arguments ..
30* REAL Z( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
40*> In case of failure it changes shifts, and tries again until output
41*> is positive.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] I0
48*> \verbatim
49*> I0 is INTEGER
50*> First index.
51*> \endverbatim
52*>
53*> \param[in,out] N0
54*> \verbatim
55*> N0 is INTEGER
56*> Last index.
57*> \endverbatim
58*>
59*> \param[in,out] Z
60*> \verbatim
61*> Z is REAL array, dimension ( 4*N0 )
62*> Z holds the qd array.
63*> \endverbatim
64*>
65*> \param[in,out] PP
66*> \verbatim
67*> PP is INTEGER
68*> PP=0 for ping, PP=1 for pong.
69*> PP=2 indicates that flipping was applied to the Z array
70*> and that the initial tests for deflation should not be
71*> performed.
72*> \endverbatim
73*>
74*> \param[out] DMIN
75*> \verbatim
76*> DMIN is REAL
77*> Minimum value of d.
78*> \endverbatim
79*>
80*> \param[out] SIGMA
81*> \verbatim
82*> SIGMA is REAL
83*> Sum of shifts used in current segment.
84*> \endverbatim
85*>
86*> \param[in,out] DESIG
87*> \verbatim
88*> DESIG is REAL
89*> Lower order part of SIGMA
90*> \endverbatim
91*>
92*> \param[in] QMAX
93*> \verbatim
94*> QMAX is REAL
95*> Maximum value of q.
96*> \endverbatim
97*>
98*> \param[in,out] NFAIL
99*> \verbatim
100*> NFAIL is INTEGER
101*> Increment NFAIL by 1 each time the shift was too big.
102*> \endverbatim
103*>
104*> \param[in,out] ITER
105*> \verbatim
106*> ITER is INTEGER
107*> Increment ITER by 1 for each iteration.
108*> \endverbatim
109*>
110*> \param[in,out] NDIV
111*> \verbatim
112*> NDIV is INTEGER
113*> Increment NDIV by 1 for each division.
114*> \endverbatim
115*>
116*> \param[in] IEEE
117*> \verbatim
118*> IEEE is LOGICAL
119*> Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
120*> \endverbatim
121*>
122*> \param[in,out] TTYPE
123*> \verbatim
124*> TTYPE is INTEGER
125*> Shift type.
126*> \endverbatim
127*>
128*> \param[in,out] DMIN1
129*> \verbatim
130*> DMIN1 is REAL
131*> \endverbatim
132*>
133*> \param[in,out] DMIN2
134*> \verbatim
135*> DMIN2 is REAL
136*> \endverbatim
137*>
138*> \param[in,out] DN
139*> \verbatim
140*> DN is REAL
141*> \endverbatim
142*>
143*> \param[in,out] DN1
144*> \verbatim
145*> DN1 is REAL
146*> \endverbatim
147*>
148*> \param[in,out] DN2
149*> \verbatim
150*> DN2 is REAL
151*> \endverbatim
152*>
153*> \param[in,out] G
154*> \verbatim
155*> G is REAL
156*> \endverbatim
157*>
158*> \param[in,out] TAU
159*> \verbatim
160*> TAU is REAL
161*>
162*> These are passed as arguments in order to save their values
163*> between calls to SLASQ3.
164*> \endverbatim
165*
166* Authors:
167* ========
168*
169*> \author Univ. of Tennessee
170*> \author Univ. of California Berkeley
171*> \author Univ. of Colorado Denver
172*> \author NAG Ltd.
173*
174*> \ingroup lasq3
175*
176* =====================================================================
177 SUBROUTINE slasq3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX,
178 $ NFAIL,
179 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
180 $ DN2, G, TAU )
181*
182* -- LAPACK computational routine --
183* -- LAPACK is a software package provided by Univ. of Tennessee, --
184* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185*
186* .. Scalar Arguments ..
187 LOGICAL IEEE
188 INTEGER I0, ITER, N0, NDIV, NFAIL, PP
189 REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
190 $ qmax, sigma, tau
191* ..
192* .. Array Arguments ..
193 REAL Z( * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL CBIAS
200 PARAMETER ( CBIAS = 1.50e0 )
201 REAL ZERO, QURTR, HALF, ONE, TWO, HUNDRD
202 parameter( zero = 0.0e0, qurtr = 0.250e0, half = 0.5e0,
203 $ one = 1.0e0, two = 2.0e0, hundrd = 100.0e0 )
204* ..
205* .. Local Scalars ..
206 INTEGER IPN4, J4, N0IN, NN, TTYPE
207 REAL EPS, S, T, TEMP, TOL, TOL2
208* ..
209* .. External Subroutines ..
210 EXTERNAL slasq4, slasq5, slasq6
211* ..
212* .. External Function ..
213 REAL SLAMCH
214 LOGICAL SISNAN
215 EXTERNAL SISNAN, SLAMCH
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC abs, max, min, sqrt
219* ..
220* .. Executable Statements ..
221*
222 n0in = n0
223 eps = slamch( 'Precision' )
224 tol = eps*hundrd
225 tol2 = tol**2
226*
227* Check for deflation.
228*
229 10 CONTINUE
230*
231 IF( n0.LT.i0 )
232 $ RETURN
233 IF( n0.EQ.i0 )
234 $ GO TO 20
235 nn = 4*n0 + pp
236 IF( n0.EQ.( i0+1 ) )
237 $ GO TO 40
238*
239* Check whether E(N0-1) is negligible, 1 eigenvalue.
240*
241 IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
242 $ z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
243 $ GO TO 30
244*
245 20 CONTINUE
246*
247 z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
248 n0 = n0 - 1
249 GO TO 10
250*
251* Check whether E(N0-2) is negligible, 2 eigenvalues.
252*
253 30 CONTINUE
254*
255 IF( z( nn-9 ).GT.tol2*sigma .AND.
256 $ z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
257 $ GO TO 50
258*
259 40 CONTINUE
260*
261 IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
262 s = z( nn-3 )
263 z( nn-3 ) = z( nn-7 )
264 z( nn-7 ) = s
265 END IF
266 t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
267 IF( z( nn-5 ).GT.z( nn-3 )*tol2.AND.t.NE.zero ) THEN
268 s = z( nn-3 )*( z( nn-5 ) / t )
269 IF( s.LE.t ) THEN
270 s = z( nn-3 )*( z( nn-5 ) /
271 $ ( t*( one+sqrt( one+s / t ) ) ) )
272 ELSE
273 s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
274 END IF
275 t = z( nn-7 ) + ( s+z( nn-5 ) )
276 z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
277 z( nn-7 ) = t
278 END IF
279 z( 4*n0-7 ) = z( nn-7 ) + sigma
280 z( 4*n0-3 ) = z( nn-3 ) + sigma
281 n0 = n0 - 2
282 GO TO 10
283*
284 50 CONTINUE
285 IF( pp.EQ.2 )
286 $ pp = 0
287*
288* Reverse the qd-array, if warranted.
289*
290 IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
291 IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
292 ipn4 = 4*( i0+n0 )
293 DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
294 temp = z( j4-3 )
295 z( j4-3 ) = z( ipn4-j4-3 )
296 z( ipn4-j4-3 ) = temp
297 temp = z( j4-2 )
298 z( j4-2 ) = z( ipn4-j4-2 )
299 z( ipn4-j4-2 ) = temp
300 temp = z( j4-1 )
301 z( j4-1 ) = z( ipn4-j4-5 )
302 z( ipn4-j4-5 ) = temp
303 temp = z( j4 )
304 z( j4 ) = z( ipn4-j4-4 )
305 z( ipn4-j4-4 ) = temp
306 60 CONTINUE
307 IF( n0-i0.LE.4 ) THEN
308 z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
309 z( 4*n0-pp ) = z( 4*i0-pp )
310 END IF
311 dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
312 z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
313 $ z( 4*i0+pp+3 ) )
314 z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
315 $ z( 4*i0-pp+4 ) )
316 qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
317 dmin = -zero
318 END IF
319 END IF
320*
321* Choose a shift.
322*
323 CALL slasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
324 $ dn2, tau, ttype, g )
325*
326* Call dqds until DMIN > 0.
327*
328 70 CONTINUE
329*
330 CALL slasq5( i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn,
331 $ dn1, dn2, ieee, eps )
332*
333 ndiv = ndiv + ( n0-i0+2 )
334 iter = iter + 1
335*
336* Check status.
337*
338 IF( dmin.GE.zero .AND. dmin1.GE.zero ) THEN
339*
340* Success.
341*
342 GO TO 90
343*
344 ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
345 $ z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
346 $ abs( dn ).LT.tol*sigma ) THEN
347*
348* Convergence hidden by negative DN.
349*
350 z( 4*( n0-1 )-pp+2 ) = zero
351 dmin = zero
352 GO TO 90
353 ELSE IF( dmin.LT.zero ) THEN
354*
355* TAU too big. Select new TAU and try again.
356*
357 nfail = nfail + 1
358 IF( ttype.LT.-22 ) THEN
359*
360* Failed twice. Play it safe.
361*
362 tau = zero
363 ELSE IF( dmin1.GT.zero ) THEN
364*
365* Late failure. Gives excellent shift.
366*
367 tau = ( tau+dmin )*( one-two*eps )
368 ttype = ttype - 11
369 ELSE
370*
371* Early failure. Divide by 4.
372*
373 tau = qurtr*tau
374 ttype = ttype - 12
375 END IF
376 GO TO 70
377 ELSE IF( sisnan( dmin ) ) THEN
378*
379* NaN.
380*
381 IF( tau.EQ.zero ) THEN
382 GO TO 80
383 ELSE
384 tau = zero
385 GO TO 70
386 END IF
387 ELSE
388*
389* Possible underflow. Play it safe.
390*
391 GO TO 80
392 END IF
393*
394* Risk of underflow.
395*
396 80 CONTINUE
397 CALL slasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
398 ndiv = ndiv + ( n0-i0+2 )
399 iter = iter + 1
400 tau = zero
401*
402 90 CONTINUE
403 IF( tau.LT.sigma ) THEN
404 desig = desig + tau
405 t = sigma + desig
406 desig = desig - ( t-sigma )
407 ELSE
408 t = sigma + tau
409 desig = sigma - ( t-tau ) + desig
410 END IF
411 sigma = t
412*
413 RETURN
414*
415* End of SLASQ3
416*
417 END
slasq3
subroutine slasq3(i0, n0, z, pp, dmin, sigma, desig, qmax, nfail, iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau)
SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
Definition slasq3.f:181
slasq4
subroutine slasq4(i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, tau, ttype, g)
SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous trans...
Definition slasq4.f:149
slasq5
subroutine slasq5(i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2, ieee, eps)
SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
Definition slasq5.f:143
slasq6
subroutine slasq6(i0, n0, z, pp, dmin, dmin1, dmin2, dn, dnm1, dnm2)
SLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
Definition slasq6.f:117