annotate libcruft/lapack/dlasq2.f @ 7034:68db500cb558

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