KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlasq3.f
Go to the documentation of this file.
1  SUBROUTINE dlasq3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
2  $ ITER, NDIV, IEEE )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * May 17, 2000
8 *
9 * .. Scalar Arguments ..
10  LOGICAL IEEE
11  INTEGER I0, ITER, N0, NDIV, NFAIL, PP
12  DOUBLE PRECISION DESIG, DMIN, QMAX, SIGMA
13 * ..
14 * .. Array Arguments ..
15  DOUBLE PRECISION Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
22 * In case of failure it changes shifts, and tries again until output
23 * is positive.
24 *
25 * Arguments
26 * =========
27 *
28 * I0 (input) INTEGER
29 * First index.
30 *
31 * N0 (input) INTEGER
32 * Last index.
33 *
34 * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
35 * Z holds the qd array.
36 *
37 * PP (input) INTEGER
38 * PP=0 for ping, PP=1 for pong.
39 *
40 * DMIN (output) DOUBLE PRECISION
41 * Minimum value of d.
42 *
43 * SIGMA (output) DOUBLE PRECISION
44 * Sum of shifts used in current segment.
45 *
46 * DESIG (input/output) DOUBLE PRECISION
47 * Lower order part of SIGMA
48 *
49 * QMAX (input) DOUBLE PRECISION
50 * Maximum value of q.
51 *
52 * NFAIL (output) INTEGER
53 * Number of times shift was too big.
54 *
55 * ITER (output) INTEGER
56 * Number of iterations.
57 *
58 * NDIV (output) INTEGER
59 * Number of divisions.
60 *
61 * TTYPE (output) INTEGER
62 * Shift type.
63 *
64 * IEEE (input) LOGICAL
65 * Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70  DOUBLE PRECISION CBIAS
71  parameter( cbias = 1.50d0 )
72  DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
73  parameter( zero = 0.0d0, qurtr = 0.250d0, half = 0.5d0,
74  $ one = 1.0d0, two = 2.0d0, hundrd = 100.0d0 )
75 * ..
76 * .. Local Scalars ..
77  INTEGER IPN4, J4, N0IN, NN, TTYPE
78  DOUBLE PRECISION DMIN1, DMIN2, DN, DN1, DN2, EPS, S, SAFMIN, T,
79  $ tau, temp, tol, tol2
80 * ..
81 * .. External Subroutines ..
82  EXTERNAL dlasq4, dlasq5, dlasq6
83 * ..
84 * .. External Function ..
85  DOUBLE PRECISION DLAMCH
86  EXTERNAL dlamch
87 * ..
88 * .. Intrinsic Functions ..
89  INTRINSIC abs, min, sqrt
90 * ..
91 * .. Save statement ..
92  SAVE ttype
93  SAVE dmin1, dmin2, dn, dn1, dn2, tau
94 * ..
95 * .. Data statement ..
96  DATA ttype / 0 /
97  DATA dmin1 / zero /, dmin2 / zero /, dn / zero /,
98  $ dn1 / zero /, dn2 / zero /, tau / zero /
99 * ..
100 * .. Executable Statements ..
101 *
102  n0in = n0
103  eps = dlamch( 'Precision' )
104  safmin = dlamch( 'Safe minimum' )
105  tol = eps*hundrd
106  tol2 = tol**2
107 *
108 * Check for deflation.
109 *
110  10 CONTINUE
111 *
112  IF( n0.LT.i0 )
113  $ RETURN
114  IF( n0.EQ.i0 )
115  $ GO TO 20
116  nn = 4*n0 + pp
117  IF( n0.EQ.( i0+1 ) )
118  $ GO TO 40
119 *
120 * Check whether E(N0-1) is negligible, 1 eigenvalue.
121 *
122  IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
123  $ z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
124  $ GO TO 30
125 *
126  20 CONTINUE
127 *
128  z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
129  n0 = n0 - 1
130  GO TO 10
131 *
132 * Check whether E(N0-2) is negligible, 2 eigenvalues.
133 *
134  30 CONTINUE
135 *
136  IF( z( nn-9 ).GT.tol2*sigma .AND.
137  $ z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
138  $ GO TO 50
139 *
140  40 CONTINUE
141 *
142  IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
143  s = z( nn-3 )
144  z( nn-3 ) = z( nn-7 )
145  z( nn-7 ) = s
146  END IF
147  IF( z( nn-5 ).GT.z( nn-3 )*tol2 ) THEN
148  t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
149  s = z( nn-3 )*( z( nn-5 ) / t )
150  IF( s.LE.t ) THEN
151  s = z( nn-3 )*( z( nn-5 ) /
152  $ ( t*( one+sqrt( one+s / t ) ) ) )
153  ELSE
154  s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
155  END IF
156  t = z( nn-7 ) + ( s+z( nn-5 ) )
157  z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
158  z( nn-7 ) = t
159  END IF
160  z( 4*n0-7 ) = z( nn-7 ) + sigma
161  z( 4*n0-3 ) = z( nn-3 ) + sigma
162  n0 = n0 - 2
163  GO TO 10
164 *
165  50 CONTINUE
166 *
167 * Reverse the qd-array, if warranted.
168 *
169  IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
170  IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
171  ipn4 = 4*( i0+n0 )
172  DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
173  temp = z( j4-3 )
174  z( j4-3 ) = z( ipn4-j4-3 )
175  z( ipn4-j4-3 ) = temp
176  temp = z( j4-2 )
177  z( j4-2 ) = z( ipn4-j4-2 )
178  z( ipn4-j4-2 ) = temp
179  temp = z( j4-1 )
180  z( j4-1 ) = z( ipn4-j4-5 )
181  z( ipn4-j4-5 ) = temp
182  temp = z( j4 )
183  z( j4 ) = z( ipn4-j4-4 )
184  z( ipn4-j4-4 ) = temp
185  60 CONTINUE
186  IF( n0-i0.LE.4 ) THEN
187  z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
188  z( 4*n0-pp ) = z( 4*i0-pp )
189  END IF
190  dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
191  z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
192  $ z( 4*i0+pp+3 ) )
193  z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
194  $ z( 4*i0-pp+4 ) )
195  qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
196  dmin = -zero
197  END IF
198  END IF
199 *
200  70 CONTINUE
201 *
202  IF( dmin.LT.zero .OR. safmin*qmax.LT.min( z( 4*n0+pp-1 ),
203  $ z( 4*n0+pp-9 ), dmin2+z( 4*n0-pp ) ) ) THEN
204 *
205 * Choose a shift.
206 *
207  CALL dlasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
208  $ dn2, tau, ttype )
209 *
210 * Call dqds until DMIN > 0.
211 *
212  80 CONTINUE
213 *
214  CALL dlasq5( i0, n0, z, pp, tau, dmin, dmin1, dmin2, dn,
215  $ dn1, dn2, ieee )
216 *
217  ndiv = ndiv + ( n0-i0+2 )
218  iter = iter + 1
219 *
220 * Check status.
221 *
222  IF( dmin.GE.zero .AND. dmin1.GT.zero ) THEN
223 *
224 * Success.
225 *
226  GO TO 100
227 *
228  ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
229  $ z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
230  $ abs( dn ).LT.tol*sigma ) THEN
231 *
232 * Convergence hidden by negative DN.
233 *
234  z( 4*( n0-1 )-pp+2 ) = zero
235  dmin = zero
236  GO TO 100
237  ELSE IF( dmin.LT.zero ) THEN
238 *
239 * TAU too big. Select new TAU and try again.
240 *
241  nfail = nfail + 1
242  IF( ttype.LT.-22 ) THEN
243 *
244 * Failed twice. Play it safe.
245 *
246  tau = zero
247  ELSE IF( dmin1.GT.zero ) THEN
248 *
249 * Late failure. Gives excellent shift.
250 *
251  tau = ( tau+dmin )*( one-two*eps )
252  ttype = ttype - 11
253  ELSE
254 *
255 * Early failure. Divide by 4.
256 *
257  tau = qurtr*tau
258  ttype = ttype - 12
259  END IF
260  GO TO 80
261  ELSE IF( dmin.NE.dmin ) THEN
262 *
263 * NaN.
264 *
265  tau = zero
266  GO TO 80
267  ELSE
268 *
269 * Possible underflow. Play it safe.
270 *
271  GO TO 90
272  END IF
273  END IF
274 *
275 * Risk of underflow.
276 *
277  90 CONTINUE
278  CALL dlasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
279  ndiv = ndiv + ( n0-i0+2 )
280  iter = iter + 1
281  tau = zero
282 *
283  100 CONTINUE
284  IF( tau.LT.sigma ) THEN
285  desig = desig + tau
286  t = sigma + desig
287  desig = desig - ( t-sigma )
288  ELSE
289  t = sigma + tau
290  desig = sigma - ( t-tau ) + desig
291  END IF
292  sigma = t
293 *
294  RETURN
295 *
296 * End of DLASQ3
297 *
298  END
dlasq3
subroutine dlasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE)
Definition: dlasq3.f:3
dlasq4
subroutine dlasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE)
Definition: dlasq4.f:3
dlasq5
subroutine dlasq5(I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE)
Definition: dlasq5.f:3
dlasq6
subroutine dlasq6(I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)
Definition: dlasq6.f:3