Return to dlasq1.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2)                                    --
    4: *
    5: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
    6: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
    7: *  -- Berkeley                                                        --
    8: *  -- November 2008                                                   --
    9: *
   10: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   11: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   12: *
   13: *     .. Scalar Arguments ..
   14:       INTEGER            INFO, N
   15: *     ..
   16: *     .. Array Arguments ..
   17:       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  DLASQ1 computes the singular values of a real N-by-N bidiagonal
   24: *  matrix with diagonal D and off-diagonal E. The singular values
   25: *  are computed to high relative accuracy, in the absence of
   26: *  denormalization, underflow and overflow. The algorithm was first
   27: *  presented in
   28: *
   29: *  "Accurate singular values and differential qd algorithms" by K. V.
   30: *  Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
   31: *  1994,
   32: *
   33: *  and the present implementation is described in "An implementation of
   34: *  the dqds Algorithm (Positive Case)", LAPACK Working Note.
   35: *
   36: *  Arguments
   37: *  =========
   38: *
   39: *  N     (input) INTEGER
   40: *        The number of rows and columns in the matrix. N >= 0.
   41: *
   42: *  D     (input/output) DOUBLE PRECISION array, dimension (N)
   43: *        On entry, D contains the diagonal elements of the
   44: *        bidiagonal matrix whose SVD is desired. On normal exit,
   45: *        D contains the singular values in decreasing order.
   46: *
   47: *  E     (input/output) DOUBLE PRECISION array, dimension (N)
   48: *        On entry, elements E(1:N-1) contain the off-diagonal elements
   49: *        of the bidiagonal matrix whose SVD is desired.
   50: *        On exit, E is overwritten.
   51: *
   52: *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
   53: *
   54: *  INFO  (output) INTEGER
   55: *        = 0: successful exit
   56: *        < 0: if INFO = -i, the i-th argument had an illegal value
   57: *        > 0: the algorithm failed
   58: *             = 1, a split was marked by a positive value in E
   59: *             = 2, current block of Z not diagonalized after 30*N
   60: *                  iterations (in inner while loop)
   61: *             = 3, termination criterion of outer while loop not met 
   62: *                  (program created more than N unreduced blocks)
   63: *
   64: *  =====================================================================
   65: *
   66: *     .. Parameters ..
   67:       DOUBLE PRECISION   ZERO
   68:       PARAMETER          ( ZERO = 0.0D0 )
   69: *     ..
   70: *     .. Local Scalars ..
   71:       INTEGER            I, IINFO
   72:       DOUBLE PRECISION   EPS, SCALE, SAFMIN, SIGMN, SIGMX
   73: *     ..
   74: *     .. External Subroutines ..
   75:       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
   76: *     ..
   77: *     .. External Functions ..
   78:       DOUBLE PRECISION   DLAMCH
   79:       EXTERNAL           DLAMCH
   80: *     ..
   81: *     .. Intrinsic Functions ..
   82:       INTRINSIC          ABS, MAX, SQRT
   83: *     ..
   84: *     .. Executable Statements ..
   85: *
   86:       INFO = 0
   87:       IF( N.LT.0 ) THEN
   88:          INFO = -2
   89:          CALL XERBLA( 'DLASQ1', -INFO )
   90:          RETURN
   91:       ELSE IF( N.EQ.0 ) THEN
   92:          RETURN
   93:       ELSE IF( N.EQ.1 ) THEN
   94:          D( 1 ) = ABS( D( 1 ) )
   95:          RETURN
   96:       ELSE IF( N.EQ.2 ) THEN
   97:          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
   98:          D( 1 ) = SIGMX
   99:          D( 2 ) = SIGMN
  100:          RETURN
  101:       END IF
  102: *
  103: *     Estimate the largest singular value.
  104: *
  105:       SIGMX = ZERO
  106:       DO 10 I = 1, N - 1
  107:          D( I ) = ABS( D( I ) )
  108:          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
  109:    10 CONTINUE
  110:       D( N ) = ABS( D( N ) )
  111: *
  112: *     Early return if SIGMX is zero (matrix is already diagonal).
  113: *
  114:       IF( SIGMX.EQ.ZERO ) THEN
  115:          CALL DLASRT( 'D', N, D, IINFO )
  116:          RETURN
  117:       END IF
  118: *
  119:       DO 20 I = 1, N
  120:          SIGMX = MAX( SIGMX, D( I ) )
  121:    20 CONTINUE
  122: *
  123: *     Copy D and E into WORK (in the Z format) and scale (squaring the
  124: *     input data makes scaling by a power of the radix pointless).
  125: *
  126:       EPS = DLAMCH( 'Precision' )
  127:       SAFMIN = DLAMCH( 'Safe minimum' )
  128:       SCALE = SQRT( EPS / SAFMIN )
  129:       CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
  130:       CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
  131:       CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
  132:      $             IINFO )
  133: *         
  134: *     Compute the q's and e's.
  135: *
  136:       DO 30 I = 1, 2*N - 1
  137:          WORK( I ) = WORK( I )**2
  138:    30 CONTINUE
  139:       WORK( 2*N ) = ZERO
  140: *
  141:       CALL DLASQ2( N, WORK, INFO )
  142: *
  143:       IF( INFO.EQ.0 ) THEN
  144:          DO 40 I = 1, N
  145:             D( I ) = SQRT( WORK( I ) )
  146:    40    CONTINUE
  147:          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
  148:       END IF
  149: *
  150:       RETURN
  151: *
  152: *     End of DLASQ1
  153: *
  154:       END