KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlasq2.f
Go to the documentation of this file.
1  SUBROUTINE dlasq2( N, Z, INFO )
2 *
3 * -- LAPACK routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * October 31, 1999
7 *
8 * .. Scalar Arguments ..
9  INTEGER INFO, N
10 * ..
11 * .. Array Arguments ..
12  DOUBLE PRECISION Z( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DLASQ2 computes all the eigenvalues of the symmetric positive
19 * definite tridiagonal matrix associated with the qd array Z to high
20 * relative accuracy are computed to high relative accuracy, in the
21 * absence of denormalization, underflow and overflow.
22 *
23 * To see the relation of Z to the tridiagonal matrix, let L be a
24 * unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
25 * let U be an upper bidiagonal matrix with 1's above and diagonal
26 * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
27 * symmetric tridiagonal to which it is similar.
28 *
29 * Note : DLASQ2 defines a logical variable, IEEE, which is true
30 * on machines which follow ieee-754 floating-point standard in their
31 * handling of infinities and NaNs, and false otherwise. This variable
32 * is passed to DLASQ3.
33 *
34 * Arguments
35 * =========
36 *
37 * N (input) INTEGER
38 * The number of rows and columns in the matrix. N >= 0.
39 *
40 * Z (workspace) DOUBLE PRECISION array, dimension ( 4*N )
41 * On entry Z holds the qd array. On exit, entries 1 to N hold
42 * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
43 * trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
44 * N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
45 * holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
46 * shifts that failed.
47 *
48 * INFO (output) INTEGER
49 * = 0: successful exit
50 * < 0: if the i-th argument is a scalar and had an illegal
51 * value, then INFO = -i, if the i-th argument is an
52 * array and the j-entry had an illegal value, then
53 * INFO = -(i*100+j)
54 * > 0: the algorithm failed
55 * = 1, a split was marked by a positive value in E
56 * = 2, current block of Z not diagonalized after 30*N
57 * iterations (in inner while loop)
58 * = 3, termination criterion of outer while loop not met
59 * (program created more than N unreduced blocks)
60 *
61 * Further Details
62 * ===============
63 * Local Variables: I0:N0 defines a current unreduced segment of Z.
64 * The shifts are accumulated in SIGMA. Iteration count is in ITER.
65 * Ping-pong is controlled by PP (alternates between 0 and 1).
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70  DOUBLE PRECISION CBIAS
71  parameter( cbias = 1.50d0 )
72  DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD
73  parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0,
74  $ two = 2.0d0, four = 4.0d0, hundrd = 100.0d0 )
75 * ..
76 * .. Local Scalars ..
77  LOGICAL IEEE
78  INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
79  $ N0, NBIG, NDIV, NFAIL, PP, SPLT
80  DOUBLE PRECISION D, DESIG, DMIN, E, EMAX, EMIN, EPS, OLDEMN,
81  $ QMAX, QMIN, S, SAFMIN, SIGMA, T, TEMP, TOL,
82  $ TOL2, TRACE, ZMAX
83 * ..
84 * .. External Subroutines ..
85  EXTERNAL dlasq3, dlasrt, xerbla
86 * ..
87 * .. External Functions ..
88  INTEGER ILAENV
89  DOUBLE PRECISION DLAMCH
90  EXTERNAL dlamch, ilaenv
91 * ..
92 * .. Intrinsic Functions ..
93  INTRINSIC dble, max, min, sqrt
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input arguments.
98 * (in case DLASQ2 is not called by DLASQ1)
99 *
100  info = 0
101  eps = dlamch( 'Precision' )
102  safmin = dlamch( 'Safe minimum' )
103  tol = eps*hundrd
104  tol2 = tol**2
105 *
106  IF( n.LT.0 ) THEN
107  info = -1
108  CALL xerbla( 'DLASQ2', 1 )
109  RETURN
110  ELSE IF( n.EQ.0 ) THEN
111  RETURN
112  ELSE IF( n.EQ.1 ) THEN
113 *
114 * 1-by-1 case.
115 *
116  IF( z( 1 ).LT.zero ) THEN
117  info = -201
118  CALL xerbla( 'DLASQ2', 2 )
119  END IF
120  RETURN
121  ELSE IF( n.EQ.2 ) THEN
122 *
123 * 2-by-2 case.
124 *
125  IF( z( 2 ).LT.zero .OR. z( 3 ).LT.zero ) THEN
126  info = -2
127  CALL xerbla( 'DLASQ2', 2 )
128  RETURN
129  ELSE IF( z( 3 ).GT.z( 1 ) ) THEN
130  d = z( 3 )
131  z( 3 ) = z( 1 )
132  z( 1 ) = d
133  END IF
134  z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
135  IF( z( 2 ).GT.z( 3 )*tol2 ) THEN
136  t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
137  s = z( 3 )*( z( 2 ) / t )
138  IF( s.LE.t ) THEN
139  s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
140  ELSE
141  s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
142  END IF
143  t = z( 1 ) + ( s+z( 2 ) )
144  z( 3 ) = z( 3 )*( z( 1 ) / t )
145  z( 1 ) = t
146  END IF
147  z( 2 ) = z( 3 )
148  z( 6 ) = z( 2 ) + z( 1 )
149  RETURN
150  END IF
151 *
152 * Check for negative data and compute sums of q's and e's.
153 *
154  z( 2*n ) = zero
155  emin = z( 2 )
156  qmax = zero
157  zmax = zero
158  d = zero
159  e = zero
160 *
161  DO 10 k = 1, 2*( n-1 ), 2
162  IF( z( k ).LT.zero ) THEN
163  info = -( 200+k )
164  CALL xerbla( 'DLASQ2', 2 )
165  RETURN
166  ELSE IF( z( k+1 ).LT.zero ) THEN
167  info = -( 200+k+1 )
168  CALL xerbla( 'DLASQ2', 2 )
169  RETURN
170  END IF
171  d = d + z( k )
172  e = e + z( k+1 )
173  qmax = max( qmax, z( k ) )
174  emin = min( emin, z( k+1 ) )
175  zmax = max( qmax, zmax, z( k+1 ) )
176  10 CONTINUE
177  IF( z( 2*n-1 ).LT.zero ) THEN
178  info = -( 200+2*n-1 )
179  CALL xerbla( 'DLASQ2', 2 )
180  RETURN
181  END IF
182  d = d + z( 2*n-1 )
183  qmax = max( qmax, z( 2*n-1 ) )
184  zmax = max( qmax, zmax )
185 *
186 * Check for diagonality.
187 *
188  IF( e.EQ.zero ) THEN
189  DO 20 k = 2, n
190  z( k ) = z( 2*k-1 )
191  20 CONTINUE
192  CALL dlasrt( 'D', n, z, iinfo )
193  z( 2*n-1 ) = d
194  RETURN
195  END IF
196 *
197  trace = d + e
198 *
199 * Check for zero data.
200 *
201  IF( trace.EQ.zero ) THEN
202  z( 2*n-1 ) = zero
203  RETURN
204  END IF
205 *
206 * Check whether the machine is IEEE conformable.
207 *
208  ieee = ilaenv( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
209  $ ilaenv( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
210 *
211 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
212 *
213  DO 30 k = 2*n, 2, -2
214  z( 2*k ) = zero
215  z( 2*k-1 ) = z( k )
216  z( 2*k-2 ) = zero
217  z( 2*k-3 ) = z( k-1 )
218  30 CONTINUE
219 *
220  i0 = 1
221  n0 = n
222 *
223 * Reverse the qd-array, if warranted.
224 *
225  IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
226  ipn4 = 4*( i0+n0 )
227  DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
228  temp = z( i4-3 )
229  z( i4-3 ) = z( ipn4-i4-3 )
230  z( ipn4-i4-3 ) = temp
231  temp = z( i4-1 )
232  z( i4-1 ) = z( ipn4-i4-5 )
233  z( ipn4-i4-5 ) = temp
234  40 CONTINUE
235  END IF
236 *
237 * Initial split checking via dqd and Li's test.
238 *
239  pp = 0
240 *
241  DO 80 k = 1, 2
242 *
243  d = z( 4*n0+pp-3 )
244  DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
245  IF( z( i4-1 ).LE.tol2*d ) THEN
246  z( i4-1 ) = -zero
247  d = z( i4-3 )
248  ELSE
249  d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
250  END IF
251  50 CONTINUE
252 *
253 * dqd maps Z to ZZ plus Li's test.
254 *
255  emin = z( 4*i0+pp+1 )
256  d = z( 4*i0+pp-3 )
257  DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
258  z( i4-2*pp-2 ) = d + z( i4-1 )
259  IF( z( i4-1 ).LE.tol2*d ) THEN
260  z( i4-1 ) = -zero
261  z( i4-2*pp-2 ) = d
262  z( i4-2*pp ) = zero
263  d = z( i4+1 )
264  ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
265  $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
266  temp = z( i4+1 ) / z( i4-2*pp-2 )
267  z( i4-2*pp ) = z( i4-1 )*temp
268  d = d*temp
269  ELSE
270  z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
271  d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
272  END IF
273  emin = min( emin, z( i4-2*pp ) )
274  60 CONTINUE
275  z( 4*n0-pp-2 ) = d
276 *
277 * Now find qmax.
278 *
279  qmax = z( 4*i0-pp-2 )
280  DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
281  qmax = max( qmax, z( i4 ) )
282  70 CONTINUE
283 *
284 * Prepare for the next iteration on K.
285 *
286  pp = 1 - pp
287  80 CONTINUE
288 *
289  iter = 2
290  nfail = 0
291  ndiv = 2*( n0-i0 )
292 *
293  DO 140 iwhila = 1, n + 1
294  IF( n0.LT.1 )
295  $ GO TO 150
296 *
297 * While array unfinished do
298 *
299 * E(N0) holds the value of SIGMA when submatrix in I0:N0
300 * splits from the rest of the array, but is negated.
301 *
302  desig = zero
303  IF( n0.EQ.n ) THEN
304  sigma = zero
305  ELSE
306  sigma = -z( 4*n0-1 )
307  END IF
308  IF( sigma.LT.zero ) THEN
309  info = 1
310  RETURN
311  END IF
312 *
313 * Find last unreduced submatrix's top index I0, find QMAX and
314 * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
315 *
316  emax = zero
317  IF( n0.GT.i0 ) THEN
318  emin = abs( z( 4*n0-5 ) )
319  ELSE
320  emin = zero
321  END IF
322  qmin = z( 4*n0-3 )
323  qmax = qmin
324  DO 90 i4 = 4*n0, 8, -4
325  IF( z( i4-5 ).LE.zero )
326  $ GO TO 100
327  IF( qmin.GE.four*emax ) THEN
328  qmin = min( qmin, z( i4-3 ) )
329  emax = max( emax, z( i4-5 ) )
330  END IF
331  qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
332  emin = min( emin, z( i4-5 ) )
333  90 CONTINUE
334  i4 = 4
335 *
336  100 CONTINUE
337  i0 = i4 / 4
338 *
339 * Store EMIN for passing to DLASQ3.
340 *
341  z( 4*n0-1 ) = emin
342 *
343 * Put -(initial shift) into DMIN.
344 *
345  dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
346 *
347 * Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.
348 *
349  pp = 0
350 *
351  nbig = 30*( n0-i0+1 )
352  DO 120 iwhilb = 1, nbig
353  IF( i0.GT.n0 )
354  $ GO TO 130
355 *
356 * While submatrix unfinished take a good dqds step.
357 *
358  CALL dlasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
359  $ iter, ndiv, ieee )
360 *
361  pp = 1 - pp
362 *
363 * When EMIN is very small check for splits.
364 *
365  IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
366  IF( z( 4*n0 ).LE.tol2*qmax .OR.
367  $ z( 4*n0-1 ).LE.tol2*sigma ) THEN
368  splt = i0 - 1
369  qmax = z( 4*i0-3 )
370  emin = z( 4*i0-1 )
371  oldemn = z( 4*i0 )
372  DO 110 i4 = 4*i0, 4*( n0-3 ), 4
373  IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
374  $ z( i4-1 ).LE.tol2*sigma ) THEN
375  z( i4-1 ) = -sigma
376  splt = i4 / 4
377  qmax = zero
378  emin = z( i4+3 )
379  oldemn = z( i4+4 )
380  ELSE
381  qmax = max( qmax, z( i4+1 ) )
382  emin = min( emin, z( i4-1 ) )
383  oldemn = min( oldemn, z( i4 ) )
384  END IF
385  110 CONTINUE
386  z( 4*n0-1 ) = emin
387  z( 4*n0 ) = oldemn
388  i0 = splt + 1
389  END IF
390  END IF
391 *
392  120 CONTINUE
393 *
394  info = 2
395  RETURN
396 *
397 * end IWHILB
398 *
399  130 CONTINUE
400 *
401  140 CONTINUE
402 *
403  info = 3
404  RETURN
405 *
406 * end IWHILA
407 *
408  150 CONTINUE
409 *
410 * Move q's to the front.
411 *
412  DO 160 k = 2, n
413  z( k ) = z( 4*k-3 )
414  160 CONTINUE
415 *
416 * Sort and compute sum of eigenvalues.
417 *
418  CALL dlasrt( 'D', n, z, iinfo )
419 *
420  e = zero
421  DO 170 k = n, 1, -1
422  e = e + z( k )
423  170 CONTINUE
424 *
425 * Store trace, sum(eigenvalues) and information on performance.
426 *
427  z( 2*n+1 ) = trace
428  z( 2*n+2 ) = e
429  z( 2*n+3 ) = dble( iter )
430  z( 2*n+4 ) = dble( ndiv ) / dble( n**2 )
431  z( 2*n+5 ) = hundrd*nfail / dble( iter )
432  RETURN
433 *
434 * End of DLASQ2
435 *
436  END
dlasq2
subroutine dlasq2(N, Z, INFO)
Definition: dlasq2.f:2
dlasq3
subroutine dlasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE)
Definition: dlasq3.f:3
dlasrt
subroutine dlasrt(ID, N, D, INFO)
Definition: dlasrt.f:2
xerbla
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2