Return to dlasq1.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dlasq1.f, revision 1.1

1.1     ! bertrand    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