001:       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2)                                    --
004: *
005: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
006: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
007: *  -- Berkeley                                                        --
008: *  -- November 2008                                                   --
009: *
010: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
011: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
012: *
013: *     .. Scalar Arguments ..
014:       INTEGER            INFO, N
015: *     ..
016: *     .. Array Arguments ..
017:       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  DLASQ1 computes the singular values of a real N-by-N bidiagonal
024: *  matrix with diagonal D and off-diagonal E. The singular values
025: *  are computed to high relative accuracy, in the absence of
026: *  denormalization, underflow and overflow. The algorithm was first
027: *  presented in
028: *
029: *  "Accurate singular values and differential qd algorithms" by K. V.
030: *  Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
031: *  1994,
032: *
033: *  and the present implementation is described in "An implementation of
034: *  the dqds Algorithm (Positive Case)", LAPACK Working Note.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  N     (input) INTEGER
040: *        The number of rows and columns in the matrix. N >= 0.
041: *
042: *  D     (input/output) DOUBLE PRECISION array, dimension (N)
043: *        On entry, D contains the diagonal elements of the
044: *        bidiagonal matrix whose SVD is desired. On normal exit,
045: *        D contains the singular values in decreasing order.
046: *
047: *  E     (input/output) DOUBLE PRECISION array, dimension (N)
048: *        On entry, elements E(1:N-1) contain the off-diagonal elements
049: *        of the bidiagonal matrix whose SVD is desired.
050: *        On exit, E is overwritten.
051: *
052: *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
053: *
054: *  INFO  (output) INTEGER
055: *        = 0: successful exit
056: *        < 0: if INFO = -i, the i-th argument had an illegal value
057: *        > 0: the algorithm failed
058: *             = 1, a split was marked by a positive value in E
059: *             = 2, current block of Z not diagonalized after 30*N
060: *                  iterations (in inner while loop)
061: *             = 3, termination criterion of outer while loop not met 
062: *                  (program created more than N unreduced blocks)
063: *
064: *  =====================================================================
065: *
066: *     .. Parameters ..
067:       DOUBLE PRECISION   ZERO
068:       PARAMETER          ( ZERO = 0.0D0 )
069: *     ..
070: *     .. Local Scalars ..
071:       INTEGER            I, IINFO
072:       DOUBLE PRECISION   EPS, SCALE, SAFMIN, SIGMN, SIGMX
073: *     ..
074: *     .. External Subroutines ..
075:       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
076: *     ..
077: *     .. External Functions ..
078:       DOUBLE PRECISION   DLAMCH
079:       EXTERNAL           DLAMCH
080: *     ..
081: *     .. Intrinsic Functions ..
082:       INTRINSIC          ABS, MAX, SQRT
083: *     ..
084: *     .. Executable Statements ..
085: *
086:       INFO = 0
087:       IF( N.LT.0 ) THEN
088:          INFO = -2
089:          CALL XERBLA( 'DLASQ1', -INFO )
090:          RETURN
091:       ELSE IF( N.EQ.0 ) THEN
092:          RETURN
093:       ELSE IF( N.EQ.1 ) THEN
094:          D( 1 ) = ABS( D( 1 ) )
095:          RETURN
096:       ELSE IF( N.EQ.2 ) THEN
097:          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
098:          D( 1 ) = SIGMX
099:          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
155: