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