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 / dlasq1.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  223 lines
  1.       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INFO, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
  13. *     ..
  14. *
  15. *     Purpose
  16. *     =======
  17. *
  18. *     DLASQ1 computes the singular values of a real N-by-N bidiagonal
  19. *     matrix with diagonal D and off-diagonal E. The singular values are
  20. *     computed to high relative accuracy, barring over/underflow or
  21. *     denormalization. The algorithm is described in
  22. *
  23. *     "Accurate singular values and differential qd algorithms," by
  24. *     K. V. Fernando and B. N. Parlett,
  25. *     Numer. Math., Vol-67, No. 2, pp. 191-230,1994.
  26. *
  27. *     See also
  28. *     "Implementation of differential qd algorithms," by
  29. *     K. V. Fernando and B. N. Parlett, Technical Report,
  30. *     Department of Mathematics, University of California at Berkeley,
  31. *     1994 (Under preparation).
  32. *
  33. *     Arguments
  34. *     =========
  35. *
  36. *  N       (input) INTEGER
  37. *          The number of rows and columns in the matrix. N >= 0.
  38. *
  39. *  D       (input/output) DOUBLE PRECISION array, dimension (N)
  40. *          On entry, D contains the diagonal elements of the
  41. *          bidiagonal matrix whose SVD is desired. On normal exit,
  42. *          D contains the singular values in decreasing order.
  43. *
  44. *  E       (input/output) DOUBLE PRECISION array, dimension (N)
  45. *          On entry, elements E(1:N-1) contain the off-diagonal elements
  46. *          of the bidiagonal matrix whose SVD is desired.
  47. *          On exit, E is overwritten.
  48. *
  49. *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
  50. *
  51. *  INFO    (output) INTEGER
  52. *          = 0:  successful exit
  53. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  54. *          > 0:  if INFO = i, the algorithm did not converge;  i
  55. *                specifies how many superdiagonals did not converge.
  56. *
  57. *  =====================================================================
  58. *
  59. *     .. Parameters ..
  60.       DOUBLE PRECISION   MEIGTH
  61.       PARAMETER          ( MEIGTH = -0.125D0 )
  62.       DOUBLE PRECISION   ZERO
  63.       PARAMETER          ( ZERO = 0.0D0 )
  64.       DOUBLE PRECISION   ONE
  65.       PARAMETER          ( ONE = 1.0D0 )
  66.       DOUBLE PRECISION   TEN
  67.       PARAMETER          ( TEN = 10.0D0 )
  68.       DOUBLE PRECISION   HUNDRD
  69.       PARAMETER          ( HUNDRD = 100.0D0 )
  70.       DOUBLE PRECISION   TWO56
  71.       PARAMETER          ( TWO56 = 256.0D0 )
  72. *     ..
  73. *     .. Local Scalars ..
  74.       LOGICAL            RESTRT
  75.       INTEGER            I, IERR, J, KE, KEND, M, NY
  76.       DOUBLE PRECISION   DM, DX, EPS, SCL, SFMIN, SIG1, SIG2, SIGMN,
  77.      $                   SIGMX, SMALL2, THRESH, TOL, TOL2, TOLMUL
  78. *     ..
  79. *     .. External Functions ..
  80.       DOUBLE PRECISION   DLAMCH
  81.       EXTERNAL           DLAMCH
  82. *     ..
  83. *     .. External Subroutines ..
  84.       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
  85. *     ..
  86. *     .. Intrinsic Functions ..
  87.       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
  88. *     ..
  89. *     .. Executable Statements ..
  90.       INFO = 0
  91.       IF( N.LT.0 ) THEN
  92.          INFO = -2
  93.          CALL XERBLA( 'DLASQ1', -INFO )
  94.          RETURN
  95.       ELSE IF( N.EQ.0 ) THEN
  96.          RETURN
  97.       ELSE IF( N.EQ.1 ) THEN
  98.          D( 1 ) = ABS( D( 1 ) )
  99.          RETURN
  100.       ELSE IF( N.EQ.2 ) THEN
  101.          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
  102.          D( 1 ) = SIGMX
  103.          D( 2 ) = SIGMN
  104.          RETURN
  105.       END IF
  106. *
  107. *     Estimate the largest singular value
  108. *
  109.       SIGMX = ZERO
  110.       DO 10 I = 1, N - 1
  111.          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
  112.    10 CONTINUE
  113. *
  114. *     Early return if sigmx is zero (matrix is already diagonal)
  115. *
  116.       IF( SIGMX.EQ.ZERO )
  117.      $   GO TO 70
  118. *
  119.       DO 20 I = 1, N
  120.          D( I ) = ABS( D( I ) )
  121.          SIGMX = MAX( SIGMX, D( I ) )
  122.    20 CONTINUE
  123. *
  124. *     Get machine parameters
  125. *
  126.       EPS = DLAMCH( 'EPSILON' )
  127.       SFMIN = DLAMCH( 'SAFE MINIMUM' )
  128. *
  129. *     Compute singular values to relative accuracy TOL
  130. *     It is assumed that tol**2 does not underflow.
  131. *
  132.       TOLMUL = MAX( TEN, MIN( HUNDRD, EPS**( -MEIGTH ) ) )
  133.       TOL = TOLMUL*EPS
  134.       TOL2 = TOL**2
  135. *
  136.       THRESH = SIGMX*SQRT( SFMIN )*TOL
  137. *
  138. *     Scale matrix so the square of the largest element is
  139. *     1 / ( 256 * SFMIN )
  140. *
  141.       SCL = SQRT( ONE / ( TWO56*SFMIN ) )
  142.       SMALL2 = ONE / ( TWO56*TOLMUL**2 )
  143.       CALL DCOPY( N, D, 1, WORK( 1 ), 1 )
  144.       CALL DCOPY( N-1, E, 1, WORK( N+1 ), 1 )
  145.       CALL DLASCL( 'G', 0, 0, SIGMX, SCL, N, 1, WORK( 1 ), N, IERR )
  146.       CALL DLASCL( 'G', 0, 0, SIGMX, SCL, N-1, 1, WORK( N+1 ), N-1,
  147.      $             IERR )
  148. *
  149. *     Square D and E (the input for the qd algorithm)
  150. *
  151.       DO 30 J = 1, 2*N - 1
  152.          WORK( J ) = WORK( J )**2
  153.    30 CONTINUE
  154. *
  155. *     Apply qd algorithm
  156. *
  157.       M = 0
  158.       E( N ) = ZERO
  159.       DX = WORK( 1 )
  160.       DM = DX
  161.       KE = 0
  162.       RESTRT = .FALSE.
  163.       DO 60 I = 1, N
  164.          IF( ABS( E( I ) ).LE.THRESH .OR. WORK( N+I ).LE.TOL2*
  165.      $       ( DM / DBLE( I-M ) ) ) THEN
  166.             NY = I - M
  167.             IF( NY.EQ.1 ) THEN
  168.                GO TO 50
  169.             ELSE IF( NY.EQ.2 ) THEN
  170.                CALL DLAS2( D( M+1 ), E( M+1 ), D( M+2 ), SIG1, SIG2 )
  171.                D( M+1 ) = SIG1
  172.                D( M+2 ) = SIG2
  173.             ELSE
  174.                KEND = KE + 1 - M
  175.                CALL DLASQ2( NY, D( M+1 ), E( M+1 ), WORK( M+1 ),
  176.      $                      WORK( M+N+1 ), EPS, TOL2, SMALL2, DM, KEND,
  177.      $                      INFO )
  178. *
  179. *                 Return, INFO = number of unconverged superdiagonals
  180. *
  181.                IF( INFO.NE.0 ) THEN
  182.                   INFO = INFO + I
  183.                   RETURN
  184.                END IF
  185. *
  186. *                 Undo scaling
  187. *
  188.                DO 40 J = M + 1, M + NY
  189.                   D( J ) = SQRT( D( J ) )
  190.    40          CONTINUE
  191.                CALL DLASCL( 'G', 0, 0, SCL, SIGMX, NY, 1, D( M+1 ), NY,
  192.      $                      IERR )
  193.             END IF
  194.    50       CONTINUE
  195.             M = I
  196.             IF( I.NE.N ) THEN
  197.                DX = WORK( I+1 )
  198.                DM = DX
  199.                KE = I
  200.                RESTRT = .TRUE.
  201.             END IF
  202.          END IF
  203.          IF( I.NE.N .AND. .NOT.RESTRT ) THEN
  204.             DX = WORK( I+1 )*( DX / ( DX+WORK( N+I ) ) )
  205.             IF( DM.GT.DX ) THEN
  206.                DM = DX
  207.                KE = I
  208.             END IF
  209.          END IF
  210.          RESTRT = .FALSE.
  211.    60 CONTINUE
  212.       KEND = KE + 1
  213. *
  214. *     Sort the singular values into decreasing order
  215. *
  216.    70 CONTINUE
  217.       CALL DLASRT( 'D', N, D, INFO )
  218.       RETURN
  219. *
  220. *     End of DLASQ1
  221. *
  222.       END
  223.