← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dbdsqr.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dbdsqr.f
 
1
*> \brief \b DBDSQR
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DBDSQR + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
 
22
*                          LDU, C, LDC, WORK, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       CHARACTER          UPLO
 
26
*       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
 
27
*       ..
 
28
*       .. Array Arguments ..
 
29
*       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
 
30
*      $                   VT( LDVT, * ), WORK( * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DBDSQR computes the singular values and, optionally, the right and/or
 
40
*> left singular vectors from the singular value decomposition (SVD) of
 
41
*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
 
42
*> zero-shift QR algorithm.  The SVD of B has the form
 
43
*> 
 
44
*>    B = Q * S * P**T
 
45
*> 
 
46
*> where S is the diagonal matrix of singular values, Q is an orthogonal
 
47
*> matrix of left singular vectors, and P is an orthogonal matrix of
 
48
*> right singular vectors.  If left singular vectors are requested, this
 
49
*> subroutine actually returns U*Q instead of Q, and, if right singular
 
50
*> vectors are requested, this subroutine returns P**T*VT instead of
 
51
*> P**T, for given real input matrices U and VT.  When U and VT are the
 
52
*> orthogonal matrices that reduce a general matrix A to bidiagonal
 
53
*> form:  A = U*B*VT, as computed by DGEBRD, then
 
54
*>
 
55
*>    A = (U*Q) * S * (P**T*VT)
 
56
*>
 
57
*> is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
 
58
*> for a given real input matrix C.
 
59
*>
 
60
*> See "Computing  Small Singular Values of Bidiagonal Matrices With
 
61
*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
 
62
*> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
 
63
*> no. 5, pp. 873-912, Sept 1990) and
 
64
*> "Accurate singular values and differential qd algorithms," by
 
65
*> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
 
66
*> Department, University of California at Berkeley, July 1992
 
67
*> for a detailed description of the algorithm.
 
68
*> \endverbatim
 
69
*
 
70
*  Arguments:
 
71
*  ==========
 
72
*
 
73
*> \param[in] UPLO
 
74
*> \verbatim
 
75
*>          UPLO is CHARACTER*1
 
76
*>          = 'U':  B is upper bidiagonal;
 
77
*>          = 'L':  B is lower bidiagonal.
 
78
*> \endverbatim
 
79
*>
 
80
*> \param[in] N
 
81
*> \verbatim
 
82
*>          N is INTEGER
 
83
*>          The order of the matrix B.  N >= 0.
 
84
*> \endverbatim
 
85
*>
 
86
*> \param[in] NCVT
 
87
*> \verbatim
 
88
*>          NCVT is INTEGER
 
89
*>          The number of columns of the matrix VT. NCVT >= 0.
 
90
*> \endverbatim
 
91
*>
 
92
*> \param[in] NRU
 
93
*> \verbatim
 
94
*>          NRU is INTEGER
 
95
*>          The number of rows of the matrix U. NRU >= 0.
 
96
*> \endverbatim
 
97
*>
 
98
*> \param[in] NCC
 
99
*> \verbatim
 
100
*>          NCC is INTEGER
 
101
*>          The number of columns of the matrix C. NCC >= 0.
 
102
*> \endverbatim
 
103
*>
 
104
*> \param[in,out] D
 
105
*> \verbatim
 
106
*>          D is DOUBLE PRECISION array, dimension (N)
 
107
*>          On entry, the n diagonal elements of the bidiagonal matrix B.
 
108
*>          On exit, if INFO=0, the singular values of B in decreasing
 
109
*>          order.
 
110
*> \endverbatim
 
111
*>
 
112
*> \param[in,out] E
 
113
*> \verbatim
 
114
*>          E is DOUBLE PRECISION array, dimension (N-1)
 
115
*>          On entry, the N-1 offdiagonal elements of the bidiagonal
 
116
*>          matrix B. 
 
117
*>          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
 
118
*>          will contain the diagonal and superdiagonal elements of a
 
119
*>          bidiagonal matrix orthogonally equivalent to the one given
 
120
*>          as input.
 
121
*> \endverbatim
 
122
*>
 
123
*> \param[in,out] VT
 
124
*> \verbatim
 
125
*>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
 
126
*>          On entry, an N-by-NCVT matrix VT.
 
127
*>          On exit, VT is overwritten by P**T * VT.
 
128
*>          Not referenced if NCVT = 0.
 
129
*> \endverbatim
 
130
*>
 
131
*> \param[in] LDVT
 
132
*> \verbatim
 
133
*>          LDVT is INTEGER
 
134
*>          The leading dimension of the array VT.
 
135
*>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
 
136
*> \endverbatim
 
137
*>
 
138
*> \param[in,out] U
 
139
*> \verbatim
 
140
*>          U is DOUBLE PRECISION array, dimension (LDU, N)
 
141
*>          On entry, an NRU-by-N matrix U.
 
142
*>          On exit, U is overwritten by U * Q.
 
143
*>          Not referenced if NRU = 0.
 
144
*> \endverbatim
 
145
*>
 
146
*> \param[in] LDU
 
147
*> \verbatim
 
148
*>          LDU is INTEGER
 
149
*>          The leading dimension of the array U.  LDU >= max(1,NRU).
 
150
*> \endverbatim
 
151
*>
 
152
*> \param[in,out] C
 
153
*> \verbatim
 
154
*>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
 
155
*>          On entry, an N-by-NCC matrix C.
 
156
*>          On exit, C is overwritten by Q**T * C.
 
157
*>          Not referenced if NCC = 0.
 
158
*> \endverbatim
 
159
*>
 
160
*> \param[in] LDC
 
161
*> \verbatim
 
162
*>          LDC is INTEGER
 
163
*>          The leading dimension of the array C.
 
164
*>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
 
165
*> \endverbatim
 
166
*>
 
167
*> \param[out] WORK
 
168
*> \verbatim
 
169
*>          WORK is DOUBLE PRECISION array, dimension (4*N)
 
170
*> \endverbatim
 
171
*>
 
172
*> \param[out] INFO
 
173
*> \verbatim
 
174
*>          INFO is INTEGER
 
175
*>          = 0:  successful exit
 
176
*>          < 0:  If INFO = -i, the i-th argument had an illegal value
 
177
*>          > 0:
 
178
*>             if NCVT = NRU = NCC = 0,
 
179
*>                = 1, a split was marked by a positive value in E
 
180
*>                = 2, current block of Z not diagonalized after 30*N
 
181
*>                     iterations (in inner while loop)
 
182
*>                = 3, termination criterion of outer while loop not met 
 
183
*>                     (program created more than N unreduced blocks)
 
184
*>             else NCVT = NRU = NCC = 0,
 
185
*>                   the algorithm did not converge; D and E contain the
 
186
*>                   elements of a bidiagonal matrix which is orthogonally
 
187
*>                   similar to the input matrix B;  if INFO = i, i
 
188
*>                   elements of E have not converged to zero.
 
189
*> \endverbatim
 
190
*
 
191
*> \par Internal Parameters:
 
192
*  =========================
 
193
*>
 
194
*> \verbatim
 
195
*>  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
 
196
*>          TOLMUL controls the convergence criterion of the QR loop.
 
197
*>          If it is positive, TOLMUL*EPS is the desired relative
 
198
*>             precision in the computed singular values.
 
199
*>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
 
200
*>             desired absolute accuracy in the computed singular
 
201
*>             values (corresponds to relative accuracy
 
202
*>             abs(TOLMUL*EPS) in the largest singular value.
 
203
*>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
 
204
*>             between 10 (for fast convergence) and .1/EPS
 
205
*>             (for there to be some accuracy in the results).
 
206
*>          Default is to lose at either one eighth or 2 of the
 
207
*>             available decimal digits in each computed singular value
 
208
*>             (whichever is smaller).
 
209
*>
 
210
*>  MAXITR  INTEGER, default = 6
 
211
*>          MAXITR controls the maximum number of passes of the
 
212
*>          algorithm through its inner loop. The algorithms stops
 
213
*>          (and so fails to converge) if the number of passes
 
214
*>          through the inner loop exceeds MAXITR*N**2.
 
215
*> \endverbatim
 
216
*
 
217
*  Authors:
 
218
*  ========
 
219
*
 
220
*> \author Univ. of Tennessee 
 
221
*> \author Univ. of California Berkeley 
 
222
*> \author Univ. of Colorado Denver 
 
223
*> \author NAG Ltd. 
 
224
*
 
225
*> \date November 2011
 
226
*
 
227
*> \ingroup auxOTHERcomputational
 
228
*
 
229
*  =====================================================================
 
230
      SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
 
231
     $                   LDU, C, LDC, WORK, INFO )
 
232
*
 
233
*  -- LAPACK computational routine (version 3.4.0) --
 
234
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
235
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
236
*     November 2011
 
237
*
 
238
*     .. Scalar Arguments ..
 
239
      CHARACTER          UPLO
 
240
      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
 
241
*     ..
 
242
*     .. Array Arguments ..
 
243
      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
 
244
     $                   VT( LDVT, * ), WORK( * )
 
245
*     ..
 
246
*
 
247
*  =====================================================================
 
248
*
 
249
*     .. Parameters ..
 
250
      DOUBLE PRECISION   ZERO
 
251
      PARAMETER          ( ZERO = 0.0D0 )
 
252
      DOUBLE PRECISION   ONE
 
253
      PARAMETER          ( ONE = 1.0D0 )
 
254
      DOUBLE PRECISION   NEGONE
 
255
      PARAMETER          ( NEGONE = -1.0D0 )
 
256
      DOUBLE PRECISION   HNDRTH
 
257
      PARAMETER          ( HNDRTH = 0.01D0 )
 
258
      DOUBLE PRECISION   TEN
 
259
      PARAMETER          ( TEN = 10.0D0 )
 
260
      DOUBLE PRECISION   HNDRD
 
261
      PARAMETER          ( HNDRD = 100.0D0 )
 
262
      DOUBLE PRECISION   MEIGTH
 
263
      PARAMETER          ( MEIGTH = -0.125D0 )
 
264
      INTEGER            MAXITR
 
265
      PARAMETER          ( MAXITR = 6 )
 
266
*     ..
 
267
*     .. Local Scalars ..
 
268
      LOGICAL            LOWER, ROTATE
 
269
      INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
 
270
     $                   NM12, NM13, OLDLL, OLDM
 
271
      DOUBLE PRECISION   ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
 
272
     $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
 
273
     $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
 
274
     $                   SN, THRESH, TOL, TOLMUL, UNFL
 
275
*     ..
 
276
*     .. External Functions ..
 
277
      LOGICAL            LSAME
 
278
      DOUBLE PRECISION   DLAMCH
 
279
      EXTERNAL           LSAME, DLAMCH
 
280
*     ..
 
281
*     .. External Subroutines ..
 
282
      EXTERNAL           DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT,
 
283
     $                   DSCAL, DSWAP, XERBLA
 
284
*     ..
 
285
*     .. Intrinsic Functions ..
 
286
      INTRINSIC          ABS, DBLE, MAX, MIN, SIGN, SQRT
 
287
*     ..
 
288
*     .. Executable Statements ..
 
289
*
 
290
*     Test the input parameters.
 
291
*
 
292
      INFO = 0
 
293
      LOWER = LSAME( UPLO, 'L' )
 
294
      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
 
295
         INFO = -1
 
296
      ELSE IF( N.LT.0 ) THEN
 
297
         INFO = -2
 
298
      ELSE IF( NCVT.LT.0 ) THEN
 
299
         INFO = -3
 
300
      ELSE IF( NRU.LT.0 ) THEN
 
301
         INFO = -4
 
302
      ELSE IF( NCC.LT.0 ) THEN
 
303
         INFO = -5
 
304
      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
 
305
     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
 
306
         INFO = -9
 
307
      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
 
308
         INFO = -11
 
309
      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
 
310
     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
 
311
         INFO = -13
 
312
      END IF
 
313
      IF( INFO.NE.0 ) THEN
 
314
         CALL XERBLA( 'DBDSQR', -INFO )
 
315
         RETURN
 
316
      END IF
 
317
      IF( N.EQ.0 )
 
318
     $   RETURN
 
319
      IF( N.EQ.1 )
 
320
     $   GO TO 160
 
321
*
 
322
*     ROTATE is true if any singular vectors desired, false otherwise
 
323
*
 
324
      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
 
325
*
 
326
*     If no singular vectors desired, use qd algorithm
 
327
*
 
328
      IF( .NOT.ROTATE ) THEN
 
329
         CALL DLASQ1( N, D, E, WORK, INFO )
 
330
*
 
331
*     If INFO equals 2, dqds didn't finish, try to finish
 
332
*         
 
333
         IF( INFO .NE. 2 ) RETURN
 
334
         INFO = 0
 
335
      END IF
 
336
*
 
337
      NM1 = N - 1
 
338
      NM12 = NM1 + NM1
 
339
      NM13 = NM12 + NM1
 
340
      IDIR = 0
 
341
*
 
342
*     Get machine constants
 
343
*
 
344
      EPS = DLAMCH( 'Epsilon' )
 
345
      UNFL = DLAMCH( 'Safe minimum' )
 
346
*
 
347
*     If matrix lower bidiagonal, rotate to be upper bidiagonal
 
348
*     by applying Givens rotations on the left
 
349
*
 
350
      IF( LOWER ) THEN
 
351
         DO 10 I = 1, N - 1
 
352
            CALL DLARTG( D( I ), E( I ), CS, SN, R )
 
353
            D( I ) = R
 
354
            E( I ) = SN*D( I+1 )
 
355
            D( I+1 ) = CS*D( I+1 )
 
356
            WORK( I ) = CS
 
357
            WORK( NM1+I ) = SN
 
358
   10    CONTINUE
 
359
*
 
360
*        Update singular vectors if desired
 
361
*
 
362
         IF( NRU.GT.0 )
 
363
     $      CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
 
364
     $                  LDU )
 
365
         IF( NCC.GT.0 )
 
366
     $      CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C,
 
367
     $                  LDC )
 
368
      END IF
 
369
*
 
370
*     Compute singular values to relative accuracy TOL
 
371
*     (By setting TOL to be negative, algorithm will compute
 
372
*     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
 
373
*
 
374
      TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
 
375
      TOL = TOLMUL*EPS
 
376
*
 
377
*     Compute approximate maximum, minimum singular values
 
378
*
 
379
      SMAX = ZERO
 
380
      DO 20 I = 1, N
 
381
         SMAX = MAX( SMAX, ABS( D( I ) ) )
 
382
   20 CONTINUE
 
383
      DO 30 I = 1, N - 1
 
384
         SMAX = MAX( SMAX, ABS( E( I ) ) )
 
385
   30 CONTINUE
 
386
      SMINL = ZERO
 
387
      IF( TOL.GE.ZERO ) THEN
 
388
*
 
389
*        Relative accuracy desired
 
390
*
 
391
         SMINOA = ABS( D( 1 ) )
 
392
         IF( SMINOA.EQ.ZERO )
 
393
     $      GO TO 50
 
394
         MU = SMINOA
 
395
         DO 40 I = 2, N
 
396
            MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
 
397
            SMINOA = MIN( SMINOA, MU )
 
398
            IF( SMINOA.EQ.ZERO )
 
399
     $         GO TO 50
 
400
   40    CONTINUE
 
401
   50    CONTINUE
 
402
         SMINOA = SMINOA / SQRT( DBLE( N ) )
 
403
         THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
 
404
      ELSE
 
405
*
 
406
*        Absolute accuracy desired
 
407
*
 
408
         THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
 
409
      END IF
 
410
*
 
411
*     Prepare for main iteration loop for the singular values
 
412
*     (MAXIT is the maximum number of passes through the inner
 
413
*     loop permitted before nonconvergence signalled.)
 
414
*
 
415
      MAXIT = MAXITR*N*N
 
416
      ITER = 0
 
417
      OLDLL = -1
 
418
      OLDM = -1
 
419
*
 
420
*     M points to last element of unconverged part of matrix
 
421
*
 
422
      M = N
 
423
*
 
424
*     Begin main iteration loop
 
425
*
 
426
   60 CONTINUE
 
427
*
 
428
*     Check for convergence or exceeding iteration count
 
429
*
 
430
      IF( M.LE.1 )
 
431
     $   GO TO 160
 
432
      IF( ITER.GT.MAXIT )
 
433
     $   GO TO 200
 
434
*
 
435
*     Find diagonal block of matrix to work on
 
436
*
 
437
      IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
 
438
     $   D( M ) = ZERO
 
439
      SMAX = ABS( D( M ) )
 
440
      SMIN = SMAX
 
441
      DO 70 LLL = 1, M - 1
 
442
         LL = M - LLL
 
443
         ABSS = ABS( D( LL ) )
 
444
         ABSE = ABS( E( LL ) )
 
445
         IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
 
446
     $      D( LL ) = ZERO
 
447
         IF( ABSE.LE.THRESH )
 
448
     $      GO TO 80
 
449
         SMIN = MIN( SMIN, ABSS )
 
450
         SMAX = MAX( SMAX, ABSS, ABSE )
 
451
   70 CONTINUE
 
452
      LL = 0
 
453
      GO TO 90
 
454
   80 CONTINUE
 
455
      E( LL ) = ZERO
 
456
*
 
457
*     Matrix splits since E(LL) = 0
 
458
*
 
459
      IF( LL.EQ.M-1 ) THEN
 
460
*
 
461
*        Convergence of bottom singular value, return to top of loop
 
462
*
 
463
         M = M - 1
 
464
         GO TO 60
 
465
      END IF
 
466
   90 CONTINUE
 
467
      LL = LL + 1
 
468
*
 
469
*     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
 
470
*
 
471
      IF( LL.EQ.M-1 ) THEN
 
472
*
 
473
*        2 by 2 block, handle separately
 
474
*
 
475
         CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
 
476
     $                COSR, SINL, COSL )
 
477
         D( M-1 ) = SIGMX
 
478
         E( M-1 ) = ZERO
 
479
         D( M ) = SIGMN
 
480
*
 
481
*        Compute singular vectors, if desired
 
482
*
 
483
         IF( NCVT.GT.0 )
 
484
     $      CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
 
485
     $                 SINR )
 
486
         IF( NRU.GT.0 )
 
487
     $      CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
 
488
         IF( NCC.GT.0 )
 
489
     $      CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
 
490
     $                 SINL )
 
491
         M = M - 2
 
492
         GO TO 60
 
493
      END IF
 
494
*
 
495
*     If working on new submatrix, choose shift direction
 
496
*     (from larger end diagonal element towards smaller)
 
497
*
 
498
      IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
 
499
         IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
 
500
*
 
501
*           Chase bulge from top (big end) to bottom (small end)
 
502
*
 
503
            IDIR = 1
 
504
         ELSE
 
505
*
 
506
*           Chase bulge from bottom (big end) to top (small end)
 
507
*
 
508
            IDIR = 2
 
509
         END IF
 
510
      END IF
 
511
*
 
512
*     Apply convergence tests
 
513
*
 
514
      IF( IDIR.EQ.1 ) THEN
 
515
*
 
516
*        Run convergence test in forward direction
 
517
*        First apply standard test to bottom of matrix
 
518
*
 
519
         IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
 
520
     $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
 
521
            E( M-1 ) = ZERO
 
522
            GO TO 60
 
523
         END IF
 
524
*
 
525
         IF( TOL.GE.ZERO ) THEN
 
526
*
 
527
*           If relative accuracy desired,
 
528
*           apply convergence criterion forward
 
529
*
 
530
            MU = ABS( D( LL ) )
 
531
            SMINL = MU
 
532
            DO 100 LLL = LL, M - 1
 
533
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
 
534
                  E( LLL ) = ZERO
 
535
                  GO TO 60
 
536
               END IF
 
537
               MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
 
538
               SMINL = MIN( SMINL, MU )
 
539
  100       CONTINUE
 
540
         END IF
 
541
*
 
542
      ELSE
 
543
*
 
544
*        Run convergence test in backward direction
 
545
*        First apply standard test to top of matrix
 
546
*
 
547
         IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
 
548
     $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
 
549
            E( LL ) = ZERO
 
550
            GO TO 60
 
551
         END IF
 
552
*
 
553
         IF( TOL.GE.ZERO ) THEN
 
554
*
 
555
*           If relative accuracy desired,
 
556
*           apply convergence criterion backward
 
557
*
 
558
            MU = ABS( D( M ) )
 
559
            SMINL = MU
 
560
            DO 110 LLL = M - 1, LL, -1
 
561
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
 
562
                  E( LLL ) = ZERO
 
563
                  GO TO 60
 
564
               END IF
 
565
               MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
 
566
               SMINL = MIN( SMINL, MU )
 
567
  110       CONTINUE
 
568
         END IF
 
569
      END IF
 
570
      OLDLL = LL
 
571
      OLDM = M
 
572
*
 
573
*     Compute shift.  First, test if shifting would ruin relative
 
574
*     accuracy, and if so set the shift to zero.
 
575
*
 
576
      IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
 
577
     $    MAX( EPS, HNDRTH*TOL ) ) THEN
 
578
*
 
579
*        Use a zero shift to avoid loss of relative accuracy
 
580
*
 
581
         SHIFT = ZERO
 
582
      ELSE
 
583
*
 
584
*        Compute the shift from 2-by-2 block at end of matrix
 
585
*
 
586
         IF( IDIR.EQ.1 ) THEN
 
587
            SLL = ABS( D( LL ) )
 
588
            CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
 
589
         ELSE
 
590
            SLL = ABS( D( M ) )
 
591
            CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
 
592
         END IF
 
593
*
 
594
*        Test if shift negligible, and if so set to zero
 
595
*
 
596
         IF( SLL.GT.ZERO ) THEN
 
597
            IF( ( SHIFT / SLL )**2.LT.EPS )
 
598
     $         SHIFT = ZERO
 
599
         END IF
 
600
      END IF
 
601
*
 
602
*     Increment iteration count
 
603
*
 
604
      ITER = ITER + M - LL
 
605
*
 
606
*     If SHIFT = 0, do simplified QR iteration
 
607
*
 
608
      IF( SHIFT.EQ.ZERO ) THEN
 
609
         IF( IDIR.EQ.1 ) THEN
 
610
*
 
611
*           Chase bulge from top to bottom
 
612
*           Save cosines and sines for later singular vector updates
 
613
*
 
614
            CS = ONE
 
615
            OLDCS = ONE
 
616
            DO 120 I = LL, M - 1
 
617
               CALL DLARTG( D( I )*CS, E( I ), CS, SN, R )
 
618
               IF( I.GT.LL )
 
619
     $            E( I-1 ) = OLDSN*R
 
620
               CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
 
621
               WORK( I-LL+1 ) = CS
 
622
               WORK( I-LL+1+NM1 ) = SN
 
623
               WORK( I-LL+1+NM12 ) = OLDCS
 
624
               WORK( I-LL+1+NM13 ) = OLDSN
 
625
  120       CONTINUE
 
626
            H = D( M )*CS
 
627
            D( M ) = H*OLDCS
 
628
            E( M-1 ) = H*OLDSN
 
629
*
 
630
*           Update singular vectors
 
631
*
 
632
            IF( NCVT.GT.0 )
 
633
     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
 
634
     $                     WORK( N ), VT( LL, 1 ), LDVT )
 
635
            IF( NRU.GT.0 )
 
636
     $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
 
637
     $                     WORK( NM13+1 ), U( 1, LL ), LDU )
 
638
            IF( NCC.GT.0 )
 
639
     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
 
640
     $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
 
641
*
 
642
*           Test convergence
 
643
*
 
644
            IF( ABS( E( M-1 ) ).LE.THRESH )
 
645
     $         E( M-1 ) = ZERO
 
646
*
 
647
         ELSE
 
648
*
 
649
*           Chase bulge from bottom to top
 
650
*           Save cosines and sines for later singular vector updates
 
651
*
 
652
            CS = ONE
 
653
            OLDCS = ONE
 
654
            DO 130 I = M, LL + 1, -1
 
655
               CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
 
656
               IF( I.LT.M )
 
657
     $            E( I ) = OLDSN*R
 
658
               CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
 
659
               WORK( I-LL ) = CS
 
660
               WORK( I-LL+NM1 ) = -SN
 
661
               WORK( I-LL+NM12 ) = OLDCS
 
662
               WORK( I-LL+NM13 ) = -OLDSN
 
663
  130       CONTINUE
 
664
            H = D( LL )*CS
 
665
            D( LL ) = H*OLDCS
 
666
            E( LL ) = H*OLDSN
 
667
*
 
668
*           Update singular vectors
 
669
*
 
670
            IF( NCVT.GT.0 )
 
671
     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
 
672
     $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
 
673
            IF( NRU.GT.0 )
 
674
     $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
 
675
     $                     WORK( N ), U( 1, LL ), LDU )
 
676
            IF( NCC.GT.0 )
 
677
     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
 
678
     $                     WORK( N ), C( LL, 1 ), LDC )
 
679
*
 
680
*           Test convergence
 
681
*
 
682
            IF( ABS( E( LL ) ).LE.THRESH )
 
683
     $         E( LL ) = ZERO
 
684
         END IF
 
685
      ELSE
 
686
*
 
687
*        Use nonzero shift
 
688
*
 
689
         IF( IDIR.EQ.1 ) THEN
 
690
*
 
691
*           Chase bulge from top to bottom
 
692
*           Save cosines and sines for later singular vector updates
 
693
*
 
694
            F = ( ABS( D( LL ) )-SHIFT )*
 
695
     $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
 
696
            G = E( LL )
 
697
            DO 140 I = LL, M - 1
 
698
               CALL DLARTG( F, G, COSR, SINR, R )
 
699
               IF( I.GT.LL )
 
700
     $            E( I-1 ) = R
 
701
               F = COSR*D( I ) + SINR*E( I )
 
702
               E( I ) = COSR*E( I ) - SINR*D( I )
 
703
               G = SINR*D( I+1 )
 
704
               D( I+1 ) = COSR*D( I+1 )
 
705
               CALL DLARTG( F, G, COSL, SINL, R )
 
706
               D( I ) = R
 
707
               F = COSL*E( I ) + SINL*D( I+1 )
 
708
               D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
 
709
               IF( I.LT.M-1 ) THEN
 
710
                  G = SINL*E( I+1 )
 
711
                  E( I+1 ) = COSL*E( I+1 )
 
712
               END IF
 
713
               WORK( I-LL+1 ) = COSR
 
714
               WORK( I-LL+1+NM1 ) = SINR
 
715
               WORK( I-LL+1+NM12 ) = COSL
 
716
               WORK( I-LL+1+NM13 ) = SINL
 
717
  140       CONTINUE
 
718
            E( M-1 ) = F
 
719
*
 
720
*           Update singular vectors
 
721
*
 
722
            IF( NCVT.GT.0 )
 
723
     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
 
724
     $                     WORK( N ), VT( LL, 1 ), LDVT )
 
725
            IF( NRU.GT.0 )
 
726
     $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
 
727
     $                     WORK( NM13+1 ), U( 1, LL ), LDU )
 
728
            IF( NCC.GT.0 )
 
729
     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
 
730
     $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
 
731
*
 
732
*           Test convergence
 
733
*
 
734
            IF( ABS( E( M-1 ) ).LE.THRESH )
 
735
     $         E( M-1 ) = ZERO
 
736
*
 
737
         ELSE
 
738
*
 
739
*           Chase bulge from bottom to top
 
740
*           Save cosines and sines for later singular vector updates
 
741
*
 
742
            F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
 
743
     $          D( M ) )
 
744
            G = E( M-1 )
 
745
            DO 150 I = M, LL + 1, -1
 
746
               CALL DLARTG( F, G, COSR, SINR, R )
 
747
               IF( I.LT.M )
 
748
     $            E( I ) = R
 
749
               F = COSR*D( I ) + SINR*E( I-1 )
 
750
               E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
 
751
               G = SINR*D( I-1 )
 
752
               D( I-1 ) = COSR*D( I-1 )
 
753
               CALL DLARTG( F, G, COSL, SINL, R )
 
754
               D( I ) = R
 
755
               F = COSL*E( I-1 ) + SINL*D( I-1 )
 
756
               D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
 
757
               IF( I.GT.LL+1 ) THEN
 
758
                  G = SINL*E( I-2 )
 
759
                  E( I-2 ) = COSL*E( I-2 )
 
760
               END IF
 
761
               WORK( I-LL ) = COSR
 
762
               WORK( I-LL+NM1 ) = -SINR
 
763
               WORK( I-LL+NM12 ) = COSL
 
764
               WORK( I-LL+NM13 ) = -SINL
 
765
  150       CONTINUE
 
766
            E( LL ) = F
 
767
*
 
768
*           Test convergence
 
769
*
 
770
            IF( ABS( E( LL ) ).LE.THRESH )
 
771
     $         E( LL ) = ZERO
 
772
*
 
773
*           Update singular vectors if desired
 
774
*
 
775
            IF( NCVT.GT.0 )
 
776
     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
 
777
     $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
 
778
            IF( NRU.GT.0 )
 
779
     $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
 
780
     $                     WORK( N ), U( 1, LL ), LDU )
 
781
            IF( NCC.GT.0 )
 
782
     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
 
783
     $                     WORK( N ), C( LL, 1 ), LDC )
 
784
         END IF
 
785
      END IF
 
786
*
 
787
*     QR iteration finished, go back and check convergence
 
788
*
 
789
      GO TO 60
 
790
*
 
791
*     All singular values converged, so make them positive
 
792
*
 
793
  160 CONTINUE
 
794
      DO 170 I = 1, N
 
795
         IF( D( I ).LT.ZERO ) THEN
 
796
            D( I ) = -D( I )
 
797
*
 
798
*           Change sign of singular vectors, if desired
 
799
*
 
800
            IF( NCVT.GT.0 )
 
801
     $         CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
 
802
         END IF
 
803
  170 CONTINUE
 
804
*
 
805
*     Sort the singular values into decreasing order (insertion sort on
 
806
*     singular values, but only one transposition per singular vector)
 
807
*
 
808
      DO 190 I = 1, N - 1
 
809
*
 
810
*        Scan for smallest D(I)
 
811
*
 
812
         ISUB = 1
 
813
         SMIN = D( 1 )
 
814
         DO 180 J = 2, N + 1 - I
 
815
            IF( D( J ).LE.SMIN ) THEN
 
816
               ISUB = J
 
817
               SMIN = D( J )
 
818
            END IF
 
819
  180    CONTINUE
 
820
         IF( ISUB.NE.N+1-I ) THEN
 
821
*
 
822
*           Swap singular values and vectors
 
823
*
 
824
            D( ISUB ) = D( N+1-I )
 
825
            D( N+1-I ) = SMIN
 
826
            IF( NCVT.GT.0 )
 
827
     $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
 
828
     $                     LDVT )
 
829
            IF( NRU.GT.0 )
 
830
     $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
 
831
            IF( NCC.GT.0 )
 
832
     $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
 
833
         END IF
 
834
  190 CONTINUE
 
835
      GO TO 220
 
836
*
 
837
*     Maximum number of iterations exceeded, failure to converge
 
838
*
 
839
  200 CONTINUE
 
840
      INFO = 0
 
841
      DO 210 I = 1, N - 1
 
842
         IF( E( I ).NE.ZERO )
 
843
     $      INFO = INFO + 1
 
844
  210 CONTINUE
 
845
  220 CONTINUE
 
846
      RETURN
 
847
*
 
848
*     End of DBDSQR
 
849
*
 
850
      END