LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \htmlonly
9 *> Download SLASQ3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASQ3( 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 * REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
29 * $ QMAX, SIGMA, TAU
30 * ..
31 * .. Array Arguments ..
32 * REAL Z( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> SLASQ3 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] Z
62 *> \verbatim
63 *> Z is REAL array, dimension ( 4*N )
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 REAL
79 *> Minimum value of d.
80 *> \endverbatim
81 *>
82 *> \param[out] SIGMA
83 *> \verbatim
84 *> SIGMA is REAL
85 *> Sum of shifts used in current segment.
86 *> \endverbatim
87 *>
88 *> \param[in,out] DESIG
89 *> \verbatim
90 *> DESIG is REAL
91 *> Lower order part of SIGMA
92 *> \endverbatim
93 *>
94 *> \param[in] QMAX
95 *> \verbatim
96 *> QMAX is REAL
97 *> Maximum value of q.
98 *> \endverbatim
99 *>
100 *> \param[out] NFAIL
101 *> \verbatim
102 *> NFAIL is INTEGER
103 *> Number of times shift was too big.
104 *> \endverbatim
105 *>
106 *> \param[out] ITER
107 *> \verbatim
108 *> ITER is INTEGER
109 *> Number of iterations.
110 *> \endverbatim
111 *>
112 *> \param[out] NDIV
113 *> \verbatim
114 *> NDIV is INTEGER
115 *> Number of divisions.
116 *> \endverbatim
117 *>
118 *> \param[in] IEEE
119 *> \verbatim
120 *> IEEE is LOGICAL
121 *> Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
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 REAL
133 *> \endverbatim
134 *>
135 *> \param[in,out] DMIN2
136 *> \verbatim
137 *> DMIN2 is REAL
138 *> \endverbatim
139 *>
140 *> \param[in,out] DN
141 *> \verbatim
142 *> DN is REAL
143 *> \endverbatim
144 *>
145 *> \param[in,out] DN1
146 *> \verbatim
147 *> DN1 is REAL
148 *> \endverbatim
149 *>
150 *> \param[in,out] DN2
151 *> \verbatim
152 *> DN2 is REAL
153 *> \endverbatim
154 *>
155 *> \param[in,out] G
156 *> \verbatim
157 *> G is REAL
158 *> \endverbatim
159 *>
160 *> \param[in,out] TAU
161 *> \verbatim
162 *> TAU is REAL
163 *>
164 *> These are passed as arguments in order to save their values
165 *> between calls to SLASQ3.
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 *> \date September 2012
177 *
178 *> \ingroup auxOTHERcomputational
179 *
180 * =====================================================================
181  SUBROUTINE slasq3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
182  $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
183  $ dn2, g, tau )
184 *
185 * -- LAPACK computational routine (version 3.4.2) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * September 2012
189 *
190 * .. Scalar Arguments ..
191  LOGICAL ieee
192  INTEGER i0, iter, n0, ndiv, nfail, pp
193  REAL desig, dmin, dmin1, dmin2, dn, dn1, dn2, g,
194  $ qmax, sigma, tau
195 * ..
196 * .. Array Arguments ..
197  REAL z( * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  REAL cbias
204  parameter( cbias = 1.50e0 )
205  REAL zero, qurtr, half, one, two, hundrd
206  parameter( zero = 0.0e0, qurtr = 0.250e0, half = 0.5e0,
207  $ one = 1.0e0, two = 2.0e0, hundrd = 100.0e0 )
208 * ..
209 * .. Local Scalars ..
210  INTEGER ipn4, j4, n0in, nn, ttype
211  REAL eps, s, t, temp, tol, tol2
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL slasq4, slasq5, slasq6
215 * ..
216 * .. External Function ..
217  REAL slamch
218  LOGICAL sisnan
219  EXTERNAL sisnan, slamch
220 * ..
221 * .. Intrinsic Functions ..
222  INTRINSIC abs, max, min, sqrt
223 * ..
224 * .. Executable Statements ..
225 *
226  n0in = n0
227  eps = slamch( 'Precision' )
228  tol = eps*hundrd
229  tol2 = tol**2
230 *
231 * Check for deflation.
232 *
233  10 continue
234 *
235  IF( n0.LT.i0 )
236  $ return
237  IF( n0.EQ.i0 )
238  $ go to 20
239  nn = 4*n0 + pp
240  IF( n0.EQ.( i0+1 ) )
241  $ go to 40
242 *
243 * Check whether E(N0-1) is negligible, 1 eigenvalue.
244 *
245  IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
246  $ z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
247  $ go to 30
248 *
249  20 continue
250 *
251  z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
252  n0 = n0 - 1
253  go to 10
254 *
255 * Check whether E(N0-2) is negligible, 2 eigenvalues.
256 *
257  30 continue
258 *
259  IF( z( nn-9 ).GT.tol2*sigma .AND.
260  $ z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
261  $ go to 50
262 *
263  40 continue
264 *
265  IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
266  s = z( nn-3 )
267  z( nn-3 ) = z( nn-7 )
268  z( nn-7 ) = s
269  END IF
270  t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
271  IF( z( nn-5 ).GT.z( nn-3 )*tol2.AND.t.NE.zero ) THEN
272  s = z( nn-3 )*( z( nn-5 ) / t )
273  IF( s.LE.t ) THEN
274  s = z( nn-3 )*( z( nn-5 ) /
275  $ ( t*( one+sqrt( one+s / t ) ) ) )
276  ELSE
277  s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
278  END IF
279  t = z( nn-7 ) + ( s+z( nn-5 ) )
280  z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
281  z( nn-7 ) = t
282  END IF
283  z( 4*n0-7 ) = z( nn-7 ) + sigma
284  z( 4*n0-3 ) = z( nn-3 ) + sigma
285  n0 = n0 - 2
286  go to 10
287 *
288  50 continue
289  IF( pp.EQ.2 )
290  $ pp = 0
291 *
292 * Reverse the qd-array, if warranted.
293 *
294  IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
295  IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
296  ipn4 = 4*( i0+n0 )
297  DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
298  temp = z( j4-3 )
299  z( j4-3 ) = z( ipn4-j4-3 )
300  z( ipn4-j4-3 ) = temp
301  temp = z( j4-2 )
302  z( j4-2 ) = z( ipn4-j4-2 )
303  z( ipn4-j4-2 ) = temp
304  temp = z( j4-1 )
305  z( j4-1 ) = z( ipn4-j4-5 )
306  z( ipn4-j4-5 ) = temp
307  temp = z( j4 )
308  z( j4 ) = z( ipn4-j4-4 )
309  z( ipn4-j4-4 ) = temp
310  60 continue
311  IF( n0-i0.LE.4 ) THEN
312  z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
313  z( 4*n0-pp ) = z( 4*i0-pp )
314  END IF
315  dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
316  z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
317  $ z( 4*i0+pp+3 ) )
318  z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
319  $ z( 4*i0-pp+4 ) )
320  qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
321  dmin = -zero
322  END IF
323  END IF
324 *
325 * Choose a shift.
326 *
327  CALL slasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
328  $ dn2, tau, ttype, g )
329 *
330 * Call dqds until DMIN > 0.
331 *
332  70 continue
333 *
334  CALL slasq5( i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn,
335  $ dn1, dn2, ieee, eps )
336 *
337  ndiv = ndiv + ( n0-i0+2 )
338  iter = iter + 1
339 *
340 * Check status.
341 *
342  IF( dmin.GE.zero .AND. dmin1.GE.zero ) THEN
343 *
344 * Success.
345 *
346  go to 90
347 *
348  ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
349  $ z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
350  $ abs( dn ).LT.tol*sigma ) THEN
351 *
352 * Convergence hidden by negative DN.
353 *
354  z( 4*( n0-1 )-pp+2 ) = zero
355  dmin = zero
356  go to 90
357  ELSE IF( dmin.LT.zero ) THEN
358 *
359 * TAU too big. Select new TAU and try again.
360 *
361  nfail = nfail + 1
362  IF( ttype.LT.-22 ) THEN
363 *
364 * Failed twice. Play it safe.
365 *
366  tau = zero
367  ELSE IF( dmin1.GT.zero ) THEN
368 *
369 * Late failure. Gives excellent shift.
370 *
371  tau = ( tau+dmin )*( one-two*eps )
372  ttype = ttype - 11
373  ELSE
374 *
375 * Early failure. Divide by 4.
376 *
377  tau = qurtr*tau
378  ttype = ttype - 12
379  END IF
380  go to 70
381  ELSE IF( sisnan( dmin ) ) THEN
382 *
383 * NaN.
384 *
385  IF( tau.EQ.zero ) THEN
386  go to 80
387  ELSE
388  tau = zero
389  go to 70
390  END IF
391  ELSE
392 *
393 * Possible underflow. Play it safe.
394 *
395  go to 80
396  END IF
397 *
398 * Risk of underflow.
399 *
400  80 continue
401  CALL slasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
402  ndiv = ndiv + ( n0-i0+2 )
403  iter = iter + 1
404  tau = zero
405 *
406  90 continue
407  IF( tau.LT.sigma ) THEN
408  desig = desig + tau
409  t = sigma + desig
410  desig = desig - ( t-sigma )
411  ELSE
412  t = sigma + tau
413  desig = sigma - ( t-tau ) + desig
414  END IF
415  sigma = t
416 *
417  return
418 *
419 * End of SLASQ3
420 *
421  END