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