LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 auxOTHERcomputational
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
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
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
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
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