← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/lapack/dlasq3.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/dlasq3.f
 
1
      SUBROUTINE DLASQ3( N, Q, E, QQ, EE, SUP, SIGMA, KEND, OFF, IPHASE,
 
2
     $                   ICONV, EPS, TOL2, SMALL2 )
 
3
*
 
4
*  -- LAPACK routine (version 2.0) --
 
5
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
6
*     Courant Institute, Argonne National Lab, and Rice University
 
7
*     September 30, 1994
 
8
*
 
9
*     .. Scalar Arguments ..
 
10
      INTEGER            ICONV, IPHASE, KEND, N, OFF
 
11
      DOUBLE PRECISION   EPS, SIGMA, SMALL2, SUP, TOL2
 
12
*     ..
 
13
*     .. Array Arguments ..
 
14
      DOUBLE PRECISION   E( * ), EE( * ), Q( * ), QQ( * )
 
15
*     ..
 
16
*
 
17
*     Purpose
 
18
*     =======
 
19
*
 
20
*     DLASQ3 is the workhorse of the whole bidiagonal SVD algorithm.
 
21
*     This can be described as the differential qd with shifts.
 
22
*
 
23
*     Arguments
 
24
*     =========
 
25
*
 
26
*  N       (input/output) INTEGER
 
27
*          On entry, N specifies the number of rows and columns
 
28
*          in the matrix. N must be at least 3.
 
29
*          On exit N is non-negative and less than the input value.
 
30
*
 
31
*  Q       (input/output) DOUBLE PRECISION array, dimension (N)
 
32
*          Q array in ping (see IPHASE below)
 
33
*
 
34
*  E       (input/output) DOUBLE PRECISION array, dimension (N)
 
35
*          E array in ping (see IPHASE below)
 
36
*
 
37
*  QQ      (input/output) DOUBLE PRECISION array, dimension (N)
 
38
*          Q array in pong (see IPHASE below)
 
39
*
 
40
*  EE      (input/output) DOUBLE PRECISION array, dimension (N)
 
41
*          E array in pong (see IPHASE below)
 
42
*
 
43
*  SUP     (input/output) DOUBLE PRECISION
 
44
*          Upper bound for the smallest eigenvalue
 
45
*
 
46
*  SIGMA   (input/output) DOUBLE PRECISION
 
47
*          Accumulated shift for the present submatrix
 
48
*
 
49
*  KEND    (input/output) INTEGER
 
50
*          Index where minimum D(i) occurs in recurrence for
 
51
*          splitting criterion
 
52
*
 
53
*  OFF     (input/output) INTEGER
 
54
*          Offset for arrays
 
55
*
 
56
*  IPHASE  (input/output) INTEGER
 
57
*          If IPHASE = 1 (ping) then data is in Q and E arrays
 
58
*          If IPHASE = 2 (pong) then data is in QQ and EE arrays
 
59
*
 
60
*  ICONV   (input) INTEGER
 
61
*          If ICONV = 0 a bottom part of a matrix (with a split)
 
62
*          If ICONV =-3 a top part of a matrix (with a split)
 
63
*
 
64
*  EPS     (input) DOUBLE PRECISION
 
65
*          Machine epsilon
 
66
*
 
67
*  TOL2    (input) DOUBLE PRECISION
 
68
*          Square of the relative tolerance TOL as defined in DLASQ1
 
69
*
 
70
*  SMALL2  (input) DOUBLE PRECISION
 
71
*          A threshold value as defined in DLASQ1
 
72
*
 
73
*  =====================================================================
 
74
*
 
75
*     .. Parameters ..
 
76
      DOUBLE PRECISION   ONE, ZERO
 
77
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 
78
      INTEGER            NPP
 
79
      PARAMETER          ( NPP = 32 )
 
80
      INTEGER            IPP
 
81
      PARAMETER          ( IPP = 5 )
 
82
      DOUBLE PRECISION   HALF, FOUR
 
83
      PARAMETER          ( HALF = 0.5D+0, FOUR = 4.0D+0 )
 
84
      INTEGER            IFLMAX
 
85
      PARAMETER          ( IFLMAX = 2 )
 
86
*     ..
 
87
*     .. Local Scalars ..
 
88
      LOGICAL            LDEF, LSPLIT
 
89
      INTEGER            I, IC, ICNT, IFL, IP, ISP, K1END, K2END, KE,
 
90
     $                   KS, MAXIT, N1, N2
 
91
      DOUBLE PRECISION   D, DM, QEMAX, T1, TAU, TOLX, TOLY, TOLZ, XX, YY
 
92
*     ..
 
93
*     .. External Subroutines ..
 
94
      EXTERNAL           DCOPY, DLASQ4
 
95
*     ..
 
96
*     .. Intrinsic Functions ..
 
97
      INTRINSIC          ABS, MAX, MIN, SQRT
 
98
*     ..
 
99
*     .. Executable Statements ..
 
100
      ICNT = 0
 
101
      TAU = ZERO
 
102
      DM = SUP
 
103
      TOLX = SIGMA*TOL2
 
104
      TOLZ = MAX( SMALL2, SIGMA )*TOL2
 
105
*
 
106
*     Set maximum number of iterations
 
107
*
 
108
      MAXIT = 100*N
 
109
*
 
110
*     Flipping
 
111
*
 
112
      IC = 2
 
113
      IF( N.GT.3 ) THEN
 
114
         IF( IPHASE.EQ.1 ) THEN
 
115
            DO 10 I = 1, N - 2
 
116
               IF( Q( I ).GT.Q( I+1 ) )
 
117
     $            IC = IC + 1
 
118
               IF( E( I ).GT.E( I+1 ) )
 
119
     $            IC = IC + 1
 
120
   10       CONTINUE
 
121
            IF( Q( N-1 ).GT.Q( N ) )
 
122
     $         IC = IC + 1
 
123
            IF( IC.LT.N ) THEN
 
124
               CALL DCOPY( N, Q, 1, QQ, -1 )
 
125
               CALL DCOPY( N-1, E, 1, EE, -1 )
 
126
               IF( KEND.NE.0 )
 
127
     $            KEND = N - KEND + 1
 
128
               IPHASE = 2
 
129
            END IF
 
130
         ELSE
 
131
            DO 20 I = 1, N - 2
 
132
               IF( QQ( I ).GT.QQ( I+1 ) )
 
133
     $            IC = IC + 1
 
134
               IF( EE( I ).GT.EE( I+1 ) )
 
135
     $            IC = IC + 1
 
136
   20       CONTINUE
 
137
            IF( QQ( N-1 ).GT.QQ( N ) )
 
138
     $         IC = IC + 1
 
139
            IF( IC.LT.N ) THEN
 
140
               CALL DCOPY( N, QQ, 1, Q, -1 )
 
141
               CALL DCOPY( N-1, EE, 1, E, -1 )
 
142
               IF( KEND.NE.0 )
 
143
     $            KEND = N - KEND + 1
 
144
               IPHASE = 1
 
145
            END IF
 
146
         END IF
 
147
      END IF
 
148
      IF( ICONV.EQ.-3 ) THEN
 
149
         IF( IPHASE.EQ.1 ) THEN
 
150
            GO TO 180
 
151
         ELSE
 
152
            GO TO 80
 
153
         END IF
 
154
      END IF
 
155
      IF( IPHASE.EQ.2 )
 
156
     $   GO TO 130
 
157
*
 
158
*     The ping section of the code
 
159
*
 
160
   30 CONTINUE
 
161
      IFL = 0
 
162
*
 
163
*     Compute the shift
 
164
*
 
165
      IF( KEND.EQ.0 .OR. SUP.EQ.ZERO ) THEN
 
166
         TAU = ZERO
 
167
      ELSE IF( ICNT.GT.0 .AND. DM.LE.TOLZ ) THEN
 
168
         TAU = ZERO
 
169
      ELSE
 
170
         IP = MAX( IPP, N / NPP )
 
171
         N2 = 2*IP + 1
 
172
         IF( N2.GE.N ) THEN
 
173
            N1 = 1
 
174
            N2 = N
 
175
         ELSE IF( KEND+IP.GT.N ) THEN
 
176
            N1 = N - 2*IP
 
177
         ELSE IF( KEND-IP.LT.1 ) THEN
 
178
            N1 = 1
 
179
         ELSE
 
180
            N1 = KEND - IP
 
181
         END IF
 
182
         CALL DLASQ4( N2, Q( N1 ), E( N1 ), TAU, SUP )
 
183
      END IF
 
184
   40 CONTINUE
 
185
      ICNT = ICNT + 1
 
186
      IF( ICNT.GT.MAXIT ) THEN
 
187
         SUP = -ONE
 
188
         RETURN
 
189
      END IF
 
190
      IF( TAU.EQ.ZERO ) THEN
 
191
*
 
192
*     dqd algorithm
 
193
*
 
194
         D = Q( 1 )
 
195
         DM = D
 
196
         KE = 0
 
197
         DO 50 I = 1, N - 3
 
198
            QQ( I ) = D + E( I )
 
199
            D = ( D / QQ( I ) )*Q( I+1 )
 
200
            IF( DM.GT.D ) THEN
 
201
               DM = D
 
202
               KE = I
 
203
            END IF
 
204
   50    CONTINUE
 
205
         KE = KE + 1
 
206
*
 
207
*     Penultimate dqd step (in ping)
 
208
*
 
209
         K2END = KE
 
210
         QQ( N-2 ) = D + E( N-2 )
 
211
         D = ( D / QQ( N-2 ) )*Q( N-1 )
 
212
         IF( DM.GT.D ) THEN
 
213
            DM = D
 
214
            KE = N - 1
 
215
         END IF
 
216
*
 
217
*     Final dqd step (in ping)
 
218
*
 
219
         K1END = KE
 
220
         QQ( N-1 ) = D + E( N-1 )
 
221
         D = ( D / QQ( N-1 ) )*Q( N )
 
222
         IF( DM.GT.D ) THEN
 
223
            DM = D
 
224
            KE = N
 
225
         END IF
 
226
         QQ( N ) = D
 
227
      ELSE
 
228
*
 
229
*     The dqds algorithm (in ping)
 
230
*
 
231
         D = Q( 1 ) - TAU
 
232
         DM = D
 
233
         KE = 0
 
234
         IF( D.LT.ZERO )
 
235
     $      GO TO 120
 
236
         DO 60 I = 1, N - 3
 
237
            QQ( I ) = D + E( I )
 
238
            D = ( D / QQ( I ) )*Q( I+1 ) - TAU
 
239
            IF( DM.GT.D ) THEN
 
240
               DM = D
 
241
               KE = I
 
242
               IF( D.LT.ZERO )
 
243
     $            GO TO 120
 
244
            END IF
 
245
   60    CONTINUE
 
246
         KE = KE + 1
 
247
*
 
248
*     Penultimate dqds step (in ping)
 
249
*
 
250
         K2END = KE
 
251
         QQ( N-2 ) = D + E( N-2 )
 
252
         D = ( D / QQ( N-2 ) )*Q( N-1 ) - TAU
 
253
         IF( DM.GT.D ) THEN
 
254
            DM = D
 
255
            KE = N - 1
 
256
            IF( D.LT.ZERO )
 
257
     $         GO TO 120
 
258
         END IF
 
259
*
 
260
*     Final dqds step (in ping)
 
261
*
 
262
         K1END = KE
 
263
         QQ( N-1 ) = D + E( N-1 )
 
264
         D = ( D / QQ( N-1 ) )*Q( N ) - TAU
 
265
         IF( DM.GT.D ) THEN
 
266
            DM = D
 
267
            KE = N
 
268
         END IF
 
269
         QQ( N ) = D
 
270
      END IF
 
271
*
 
272
*        Convergence when QQ(N) is small (in ping)
 
273
*
 
274
      IF( ABS( QQ( N ) ).LE.SIGMA*TOL2 ) THEN
 
275
         QQ( N ) = ZERO
 
276
         DM = ZERO
 
277
         KE = N
 
278
      END IF
 
279
      IF( QQ( N ).LT.ZERO )
 
280
     $   GO TO 120
 
281
*
 
282
*     Non-negative qd array: Update the e's
 
283
*
 
284
      DO 70 I = 1, N - 1
 
285
         EE( I ) = ( E( I ) / QQ( I ) )*Q( I+1 )
 
286
   70 CONTINUE
 
287
*
 
288
*     Updating sigma and iphase in ping
 
289
*
 
290
      SIGMA = SIGMA + TAU
 
291
      IPHASE = 2
 
292
   80 CONTINUE
 
293
      TOLX = SIGMA*TOL2
 
294
      TOLY = SIGMA*EPS
 
295
      TOLZ = MAX( SIGMA, SMALL2 )*TOL2
 
296
*
 
297
*     Checking for deflation and convergence (in ping)
 
298
*
 
299
   90 CONTINUE
 
300
      IF( N.LE.2 )
 
301
     $   RETURN
 
302
*
 
303
*        Deflation: bottom 1x1 (in ping)
 
304
*
 
305
      LDEF = .FALSE.
 
306
      IF( EE( N-1 ).LE.TOLZ ) THEN
 
307
         LDEF = .TRUE.
 
308
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
309
         IF( EE( N-1 ).LE.EPS*( SIGMA+QQ( N ) ) ) THEN
 
310
            IF( EE( N-1 )*( QQ( N ) / ( QQ( N )+SIGMA ) ).LE.TOL2*
 
311
     $          ( QQ( N )+SIGMA ) ) THEN
 
312
               LDEF = .TRUE.
 
313
            END IF
 
314
         END IF
 
315
      ELSE
 
316
         IF( EE( N-1 ).LE.QQ( N )*TOL2 ) THEN
 
317
            LDEF = .TRUE.
 
318
         END IF
 
319
      END IF
 
320
      IF( LDEF ) THEN
 
321
         Q( N ) = QQ( N ) + SIGMA
 
322
         N = N - 1
 
323
         ICONV = ICONV + 1
 
324
         GO TO 90
 
325
      END IF
 
326
*
 
327
*        Deflation: bottom 2x2 (in ping)
 
328
*
 
329
      LDEF = .FALSE.
 
330
      IF( EE( N-2 ).LE.TOLZ ) THEN
 
331
         LDEF = .TRUE.
 
332
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
333
         T1 = SIGMA + EE( N-1 )*( SIGMA / ( SIGMA+QQ( N ) ) )
 
334
         IF( EE( N-2 )*( T1 / ( QQ( N-1 )+T1 ) ).LE.TOLY ) THEN
 
335
            IF( EE( N-2 )*( QQ( N-1 ) / ( QQ( N-1 )+T1 ) ).LE.TOLX )
 
336
     $           THEN
 
337
               LDEF = .TRUE.
 
338
            END IF
 
339
         END IF
 
340
      ELSE
 
341
         IF( EE( N-2 ).LE.( QQ( N ) / ( QQ( N )+EE( N-1 )+QQ( N-1 ) ) )*
 
342
     $       QQ( N-1 )*TOL2 ) THEN
 
343
            LDEF = .TRUE.
 
344
         END IF
 
345
      END IF
 
346
      IF( LDEF ) THEN
 
347
         QEMAX = MAX( QQ( N ), QQ( N-1 ), EE( N-1 ) )
 
348
         IF( QEMAX.NE.ZERO ) THEN
 
349
            IF( QEMAX.EQ.QQ( N-1 ) ) THEN
 
350
               XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
 
351
     $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
 
352
     $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
 
353
            ELSE IF( QEMAX.EQ.QQ( N ) ) THEN
 
354
               XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
 
355
     $              SQRT( ( ( QQ( N-1 )-QQ( N )+EE( N-1 ) ) /
 
356
     $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
 
357
            ELSE
 
358
               XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
 
359
     $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
 
360
     $              QEMAX )**2+FOUR*QQ( N-1 ) / QEMAX ) )
 
361
            END IF
 
362
            YY = ( MAX( QQ( N ), QQ( N-1 ) ) / XX )*
 
363
     $           MIN( QQ( N ), QQ( N-1 ) )
 
364
         ELSE
 
365
            XX = ZERO
 
366
            YY = ZERO
 
367
         END IF
 
368
         Q( N-1 ) = SIGMA + XX
 
369
         Q( N ) = YY + SIGMA
 
370
         N = N - 2
 
371
         ICONV = ICONV + 2
 
372
         GO TO 90
 
373
      END IF
 
374
*
 
375
*     Updating bounds before going to pong
 
376
*
 
377
      IF( ICONV.EQ.0 ) THEN
 
378
         KEND = KE
 
379
         SUP = MIN( DM, SUP-TAU )
 
380
      ELSE IF( ICONV.GT.0 ) THEN
 
381
         SUP = MIN( QQ( N ), QQ( N-1 ), QQ( N-2 ), QQ( 1 ), QQ( 2 ),
 
382
     $         QQ( 3 ) )
 
383
         IF( ICONV.EQ.1 ) THEN
 
384
            KEND = K1END
 
385
         ELSE IF( ICONV.EQ.2 ) THEN
 
386
            KEND = K2END
 
387
         ELSE
 
388
            KEND = N
 
389
         END IF
 
390
         ICNT = 0
 
391
         MAXIT = 100*N
 
392
      END IF
 
393
*
 
394
*     Checking for splitting in ping
 
395
*
 
396
      LSPLIT = .FALSE.
 
397
      DO 100 KS = N - 3, 3, -1
 
398
         IF( EE( KS ).LE.TOLY ) THEN
 
399
            IF( EE( KS )*( MIN( QQ( KS+1 ),
 
400
     $          QQ( KS ) ) / ( MIN( QQ( KS+1 ), QQ( KS ) )+SIGMA ) ).LE.
 
401
     $          TOLX ) THEN
 
402
               LSPLIT = .TRUE.
 
403
               GO TO 110
 
404
            END IF
 
405
         END IF
 
406
  100 CONTINUE
 
407
*
 
408
      KS = 2
 
409
      IF( EE( 2 ).LE.TOLZ ) THEN
 
410
         LSPLIT = .TRUE.
 
411
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
412
         T1 = SIGMA + EE( 1 )*( SIGMA / ( SIGMA+QQ( 1 ) ) )
 
413
         IF( EE( 2 )*( T1 / ( QQ( 1 )+T1 ) ).LE.TOLY ) THEN
 
414
            IF( EE( 2 )*( QQ( 1 ) / ( QQ( 1 )+T1 ) ).LE.TOLX ) THEN
 
415
               LSPLIT = .TRUE.
 
416
            END IF
 
417
         END IF
 
418
      ELSE
 
419
         IF( EE( 2 ).LE.( QQ( 1 ) / ( QQ( 1 )+EE( 1 )+QQ( 2 ) ) )*
 
420
     $       QQ( 2 )*TOL2 ) THEN
 
421
            LSPLIT = .TRUE.
 
422
         END IF
 
423
      END IF
 
424
      IF( LSPLIT )
 
425
     $   GO TO 110
 
426
*
 
427
      KS = 1
 
428
      IF( EE( 1 ).LE.TOLZ ) THEN
 
429
         LSPLIT = .TRUE.
 
430
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
431
         IF( EE( 1 ).LE.EPS*( SIGMA+QQ( 1 ) ) ) THEN
 
432
            IF( EE( 1 )*( QQ( 1 ) / ( QQ( 1 )+SIGMA ) ).LE.TOL2*
 
433
     $          ( QQ( 1 )+SIGMA ) ) THEN
 
434
               LSPLIT = .TRUE.
 
435
            END IF
 
436
         END IF
 
437
      ELSE
 
438
         IF( EE( 1 ).LE.QQ( 1 )*TOL2 ) THEN
 
439
            LSPLIT = .TRUE.
 
440
         END IF
 
441
      END IF
 
442
*
 
443
  110 CONTINUE
 
444
      IF( LSPLIT ) THEN
 
445
         SUP = MIN( QQ( N ), QQ( N-1 ), QQ( N-2 ) )
 
446
         ISP = -( OFF+1 )
 
447
         OFF = OFF + KS
 
448
         N = N - KS
 
449
         KEND = MAX( 1, KEND-KS )
 
450
         E( KS ) = SIGMA
 
451
         EE( KS ) = ISP
 
452
         ICONV = 0
 
453
         RETURN
 
454
      END IF
 
455
*
 
456
*     Coincidence
 
457
*
 
458
      IF( TAU.EQ.ZERO .AND. DM.LE.TOLZ .AND. KEND.NE.N .AND. ICONV.EQ.
 
459
     $    0 .AND. ICNT.GT.0 ) THEN
 
460
         CALL DCOPY( N-KE, E( KE ), 1, QQ( KE ), 1 )
 
461
         QQ( N ) = ZERO
 
462
         CALL DCOPY( N-KE, Q( KE+1 ), 1, EE( KE ), 1 )
 
463
         SUP = ZERO
 
464
      END IF
 
465
      ICONV = 0
 
466
      GO TO 130
 
467
*
 
468
*     A new shift when the previous failed (in ping)
 
469
*
 
470
  120 CONTINUE
 
471
      IFL = IFL + 1
 
472
      SUP = TAU
 
473
*
 
474
*     SUP is small or
 
475
*     Too many bad shifts (ping)
 
476
*
 
477
      IF( SUP.LE.TOLZ .OR. IFL.GE.IFLMAX ) THEN
 
478
         TAU = ZERO
 
479
         GO TO 40
 
480
*
 
481
*     The asymptotic shift (in ping)
 
482
*
 
483
      ELSE
 
484
         TAU = MAX( TAU+D, ZERO )
 
485
         IF( TAU.LE.TOLZ )
 
486
     $      TAU = ZERO
 
487
         GO TO 40
 
488
      END IF
 
489
*
 
490
*     the pong section of the code
 
491
*
 
492
  130 CONTINUE
 
493
      IFL = 0
 
494
*
 
495
*     Compute the shift (in pong)
 
496
*
 
497
      IF( KEND.EQ.0 .AND. SUP.EQ.ZERO ) THEN
 
498
         TAU = ZERO
 
499
      ELSE IF( ICNT.GT.0 .AND. DM.LE.TOLZ ) THEN
 
500
         TAU = ZERO
 
501
      ELSE
 
502
         IP = MAX( IPP, N / NPP )
 
503
         N2 = 2*IP + 1
 
504
         IF( N2.GE.N ) THEN
 
505
            N1 = 1
 
506
            N2 = N
 
507
         ELSE IF( KEND+IP.GT.N ) THEN
 
508
            N1 = N - 2*IP
 
509
         ELSE IF( KEND-IP.LT.1 ) THEN
 
510
            N1 = 1
 
511
         ELSE
 
512
            N1 = KEND - IP
 
513
         END IF
 
514
         CALL DLASQ4( N2, QQ( N1 ), EE( N1 ), TAU, SUP )
 
515
      END IF
 
516
  140 CONTINUE
 
517
      ICNT = ICNT + 1
 
518
      IF( ICNT.GT.MAXIT ) THEN
 
519
         SUP = -SUP
 
520
         RETURN
 
521
      END IF
 
522
      IF( TAU.EQ.ZERO ) THEN
 
523
*
 
524
*     The dqd algorithm (in pong)
 
525
*
 
526
         D = QQ( 1 )
 
527
         DM = D
 
528
         KE = 0
 
529
         DO 150 I = 1, N - 3
 
530
            Q( I ) = D + EE( I )
 
531
            D = ( D / Q( I ) )*QQ( I+1 )
 
532
            IF( DM.GT.D ) THEN
 
533
               DM = D
 
534
               KE = I
 
535
            END IF
 
536
  150    CONTINUE
 
537
         KE = KE + 1
 
538
*
 
539
*     Penultimate dqd step (in pong)
 
540
*
 
541
         K2END = KE
 
542
         Q( N-2 ) = D + EE( N-2 )
 
543
         D = ( D / Q( N-2 ) )*QQ( N-1 )
 
544
         IF( DM.GT.D ) THEN
 
545
            DM = D
 
546
            KE = N - 1
 
547
         END IF
 
548
*
 
549
*     Final dqd step (in pong)
 
550
*
 
551
         K1END = KE
 
552
         Q( N-1 ) = D + EE( N-1 )
 
553
         D = ( D / Q( N-1 ) )*QQ( N )
 
554
         IF( DM.GT.D ) THEN
 
555
            DM = D
 
556
            KE = N
 
557
         END IF
 
558
         Q( N ) = D
 
559
      ELSE
 
560
*
 
561
*     The dqds algorithm (in pong)
 
562
*
 
563
         D = QQ( 1 ) - TAU
 
564
         DM = D
 
565
         KE = 0
 
566
         IF( D.LT.ZERO )
 
567
     $      GO TO 220
 
568
         DO 160 I = 1, N - 3
 
569
            Q( I ) = D + EE( I )
 
570
            D = ( D / Q( I ) )*QQ( I+1 ) - TAU
 
571
            IF( DM.GT.D ) THEN
 
572
               DM = D
 
573
               KE = I
 
574
               IF( D.LT.ZERO )
 
575
     $            GO TO 220
 
576
            END IF
 
577
  160    CONTINUE
 
578
         KE = KE + 1
 
579
*
 
580
*     Penultimate dqds step (in pong)
 
581
*
 
582
         K2END = KE
 
583
         Q( N-2 ) = D + EE( N-2 )
 
584
         D = ( D / Q( N-2 ) )*QQ( N-1 ) - TAU
 
585
         IF( DM.GT.D ) THEN
 
586
            DM = D
 
587
            KE = N - 1
 
588
            IF( D.LT.ZERO )
 
589
     $         GO TO 220
 
590
         END IF
 
591
*
 
592
*     Final dqds step (in pong)
 
593
*
 
594
         K1END = KE
 
595
         Q( N-1 ) = D + EE( N-1 )
 
596
         D = ( D / Q( N-1 ) )*QQ( N ) - TAU
 
597
         IF( DM.GT.D ) THEN
 
598
            DM = D
 
599
            KE = N
 
600
         END IF
 
601
         Q( N ) = D
 
602
      END IF
 
603
*
 
604
*        Convergence when is small (in pong)
 
605
*
 
606
      IF( ABS( Q( N ) ).LE.SIGMA*TOL2 ) THEN
 
607
         Q( N ) = ZERO
 
608
         DM = ZERO
 
609
         KE = N
 
610
      END IF
 
611
      IF( Q( N ).LT.ZERO )
 
612
     $   GO TO 220
 
613
*
 
614
*     Non-negative qd array: Update the e's
 
615
*
 
616
      DO 170 I = 1, N - 1
 
617
         E( I ) = ( EE( I ) / Q( I ) )*QQ( I+1 )
 
618
  170 CONTINUE
 
619
*
 
620
*     Updating sigma and iphase in pong
 
621
*
 
622
      SIGMA = SIGMA + TAU
 
623
  180 CONTINUE
 
624
      IPHASE = 1
 
625
      TOLX = SIGMA*TOL2
 
626
      TOLY = SIGMA*EPS
 
627
*
 
628
*     Checking for deflation and convergence (in pong)
 
629
*
 
630
  190 CONTINUE
 
631
      IF( N.LE.2 )
 
632
     $   RETURN
 
633
*
 
634
*        Deflation: bottom 1x1 (in pong)
 
635
*
 
636
      LDEF = .FALSE.
 
637
      IF( E( N-1 ).LE.TOLZ ) THEN
 
638
         LDEF = .TRUE.
 
639
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
640
         IF( E( N-1 ).LE.EPS*( SIGMA+Q( N ) ) ) THEN
 
641
            IF( E( N-1 )*( Q( N ) / ( Q( N )+SIGMA ) ).LE.TOL2*
 
642
     $          ( Q( N )+SIGMA ) ) THEN
 
643
               LDEF = .TRUE.
 
644
            END IF
 
645
         END IF
 
646
      ELSE
 
647
         IF( E( N-1 ).LE.Q( N )*TOL2 ) THEN
 
648
            LDEF = .TRUE.
 
649
         END IF
 
650
      END IF
 
651
      IF( LDEF ) THEN
 
652
         Q( N ) = Q( N ) + SIGMA
 
653
         N = N - 1
 
654
         ICONV = ICONV + 1
 
655
         GO TO 190
 
656
      END IF
 
657
*
 
658
*        Deflation: bottom 2x2 (in pong)
 
659
*
 
660
      LDEF = .FALSE.
 
661
      IF( E( N-2 ).LE.TOLZ ) THEN
 
662
         LDEF = .TRUE.
 
663
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
664
         T1 = SIGMA + E( N-1 )*( SIGMA / ( SIGMA+Q( N ) ) )
 
665
         IF( E( N-2 )*( T1 / ( Q( N-1 )+T1 ) ).LE.TOLY ) THEN
 
666
            IF( E( N-2 )*( Q( N-1 ) / ( Q( N-1 )+T1 ) ).LE.TOLX ) THEN
 
667
               LDEF = .TRUE.
 
668
            END IF
 
669
         END IF
 
670
      ELSE
 
671
         IF( E( N-2 ).LE.( Q( N ) / ( Q( N )+EE( N-1 )+Q( N-1 ) )*Q( N-
 
672
     $       1 ) )*TOL2 ) THEN
 
673
            LDEF = .TRUE.
 
674
         END IF
 
675
      END IF
 
676
      IF( LDEF ) THEN
 
677
         QEMAX = MAX( Q( N ), Q( N-1 ), E( N-1 ) )
 
678
         IF( QEMAX.NE.ZERO ) THEN
 
679
            IF( QEMAX.EQ.Q( N-1 ) ) THEN
 
680
               XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
 
681
     $              SQRT( ( ( Q( N )-Q( N-1 )+E( N-1 ) ) / QEMAX )**2+
 
682
     $              FOUR*E( N-1 ) / QEMAX ) )
 
683
            ELSE IF( QEMAX.EQ.Q( N ) ) THEN
 
684
               XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
 
685
     $              SQRT( ( ( Q( N-1 )-Q( N )+E( N-1 ) ) / QEMAX )**2+
 
686
     $              FOUR*E( N-1 ) / QEMAX ) )
 
687
            ELSE
 
688
               XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
 
689
     $              SQRT( ( ( Q( N )-Q( N-1 )+E( N-1 ) ) / QEMAX )**2+
 
690
     $              FOUR*Q( N-1 ) / QEMAX ) )
 
691
            END IF
 
692
            YY = ( MAX( Q( N ), Q( N-1 ) ) / XX )*
 
693
     $           MIN( Q( N ), Q( N-1 ) )
 
694
         ELSE
 
695
            XX = ZERO
 
696
            YY = ZERO
 
697
         END IF
 
698
         Q( N-1 ) = SIGMA + XX
 
699
         Q( N ) = YY + SIGMA
 
700
         N = N - 2
 
701
         ICONV = ICONV + 2
 
702
         GO TO 190
 
703
      END IF
 
704
*
 
705
*     Updating bounds before going to pong
 
706
*
 
707
      IF( ICONV.EQ.0 ) THEN
 
708
         KEND = KE
 
709
         SUP = MIN( DM, SUP-TAU )
 
710
      ELSE IF( ICONV.GT.0 ) THEN
 
711
         SUP = MIN( Q( N ), Q( N-1 ), Q( N-2 ), Q( 1 ), Q( 2 ), Q( 3 ) )
 
712
         IF( ICONV.EQ.1 ) THEN
 
713
            KEND = K1END
 
714
         ELSE IF( ICONV.EQ.2 ) THEN
 
715
            KEND = K2END
 
716
         ELSE
 
717
            KEND = N
 
718
         END IF
 
719
         ICNT = 0
 
720
         MAXIT = 100*N
 
721
      END IF
 
722
*
 
723
*     Checking for splitting in pong
 
724
*
 
725
      LSPLIT = .FALSE.
 
726
      DO 200 KS = N - 3, 3, -1
 
727
         IF( E( KS ).LE.TOLY ) THEN
 
728
            IF( E( KS )*( MIN( Q( KS+1 ), Q( KS ) ) / ( MIN( Q( KS+1 ),
 
729
     $          Q( KS ) )+SIGMA ) ).LE.TOLX ) THEN
 
730
               LSPLIT = .TRUE.
 
731
               GO TO 210
 
732
            END IF
 
733
         END IF
 
734
  200 CONTINUE
 
735
*
 
736
      KS = 2
 
737
      IF( E( 2 ).LE.TOLZ ) THEN
 
738
         LSPLIT = .TRUE.
 
739
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
740
         T1 = SIGMA + E( 1 )*( SIGMA / ( SIGMA+Q( 1 ) ) )
 
741
         IF( E( 2 )*( T1 / ( Q( 1 )+T1 ) ).LE.TOLY ) THEN
 
742
            IF( E( 2 )*( Q( 1 ) / ( Q( 1 )+T1 ) ).LE.TOLX ) THEN
 
743
               LSPLIT = .TRUE.
 
744
            END IF
 
745
         END IF
 
746
      ELSE
 
747
         IF( E( 2 ).LE.( Q( 1 ) / ( Q( 1 )+E( 1 )+Q( 2 ) ) )*Q( 2 )*
 
748
     $       TOL2 ) THEN
 
749
            LSPLIT = .TRUE.
 
750
         END IF
 
751
      END IF
 
752
      IF( LSPLIT )
 
753
     $   GO TO 210
 
754
*
 
755
      KS = 1
 
756
      IF( E( 1 ).LE.TOLZ ) THEN
 
757
         LSPLIT = .TRUE.
 
758
      ELSE IF( SIGMA.GT.ZERO ) THEN
 
759
         IF( E( 1 ).LE.EPS*( SIGMA+Q( 1 ) ) ) THEN
 
760
            IF( E( 1 )*( Q( 1 ) / ( Q( 1 )+SIGMA ) ).LE.TOL2*
 
761
     $          ( Q( 1 )+SIGMA ) ) THEN
 
762
               LSPLIT = .TRUE.
 
763
            END IF
 
764
         END IF
 
765
      ELSE
 
766
         IF( E( 1 ).LE.Q( 1 )*TOL2 ) THEN
 
767
            LSPLIT = .TRUE.
 
768
         END IF
 
769
      END IF
 
770
*
 
771
  210 CONTINUE
 
772
      IF( LSPLIT ) THEN
 
773
         SUP = MIN( Q( N ), Q( N-1 ), Q( N-2 ) )
 
774
         ISP = OFF + 1
 
775
         OFF = OFF + KS
 
776
         KEND = MAX( 1, KEND-KS )
 
777
         N = N - KS
 
778
         E( KS ) = SIGMA
 
779
         EE( KS ) = ISP
 
780
         ICONV = 0
 
781
         RETURN
 
782
      END IF
 
783
*
 
784
*     Coincidence
 
785
*
 
786
      IF( TAU.EQ.ZERO .AND. DM.LE.TOLZ .AND. KEND.NE.N .AND. ICONV.EQ.
 
787
     $    0 .AND. ICNT.GT.0 ) THEN
 
788
         CALL DCOPY( N-KE, EE( KE ), 1, Q( KE ), 1 )
 
789
         Q( N ) = ZERO
 
790
         CALL DCOPY( N-KE, QQ( KE+1 ), 1, E( KE ), 1 )
 
791
         SUP = ZERO
 
792
      END IF
 
793
      ICONV = 0
 
794
      GO TO 30
 
795
*
 
796
*     Computation of a new shift when the previous failed (in pong)
 
797
*
 
798
  220 CONTINUE
 
799
      IFL = IFL + 1
 
800
      SUP = TAU
 
801
*
 
802
*     SUP is small or
 
803
*     Too many bad shifts (in pong)
 
804
*
 
805
      IF( SUP.LE.TOLZ .OR. IFL.GE.IFLMAX ) THEN
 
806
         TAU = ZERO
 
807
         GO TO 140
 
808
*
 
809
*     The asymptotic shift (in pong)
 
810
*
 
811
      ELSE
 
812
         TAU = MAX( TAU+D, ZERO )
 
813
         IF( TAU.LE.TOLZ )
 
814
     $      TAU = ZERO
 
815
         GO TO 140
 
816
      END IF
 
817
*
 
818
*     End of DLASQ3
 
819
*
 
820
      END