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