SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
*
* -- LAPACK routine (version 3.2) --
*
* -- Contributed by Osni Marques of the Lawrence Berkeley National --
* -- Laboratory and Beresford Parlett of the Univ. of California at --
* -- Berkeley --
* -- November 2008 --
*
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INFO, N
* ..
* .. Array Arguments ..
DOUBLE PRECISION D( * ), E( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* DLASQ1 computes the singular values of a real N-by-N bidiagonal
* matrix with diagonal D and off-diagonal E. The singular values
* are computed to high relative accuracy, in the absence of
* denormalization, underflow and overflow. The algorithm was first
* presented in
*
* "Accurate singular values and differential qd algorithms" by K. V.
* Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
* 1994,
*
* and the present implementation is described in "An implementation of
* the dqds Algorithm (Positive Case)", LAPACK Working Note.
*
* Arguments
* =========
*
* N (input) INTEGER
* The number of rows and columns in the matrix. N >= 0.
*
* D (input/output) DOUBLE PRECISION array, dimension (N)
* On entry, D contains the diagonal elements of the
* bidiagonal matrix whose SVD is desired. On normal exit,
* D contains the singular values in decreasing order.
*
* E (input/output) DOUBLE PRECISION array, dimension (N)
* On entry, elements E(1:N-1) contain the off-diagonal elements
* of the bidiagonal matrix whose SVD is desired.
* On exit, E is overwritten.
*
* WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: the algorithm failed
* = 1, a split was marked by a positive value in E
* = 2, current block of Z not diagonalized after 30*N
* iterations (in inner while loop)
* = 3, termination criterion of outer while loop not met
* (program created more than N unreduced blocks)
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
* ..
* .. Local Scalars ..
INTEGER I, IINFO
DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT
* ..
* .. Executable Statements ..
*
INFO = 0
IF( N.LT.0 ) THEN
INFO = -2
CALL XERBLA( 'DLASQ1', -INFO )
RETURN
ELSE IF( N.EQ.0 ) THEN
RETURN
ELSE IF( N.EQ.1 ) THEN
D( 1 ) = ABS( D( 1 ) )
RETURN
ELSE IF( N.EQ.2 ) THEN
CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
D( 1 ) = SIGMX
D( 2 ) = SIGMN
RETURN
END IF
*
* Estimate the largest singular value.
*
SIGMX = ZERO
DO 10 I = 1, N - 1
D( I ) = ABS( D( I ) )
SIGMX = MAX( SIGMX, ABS( E( I ) ) )
10 CONTINUE
D( N ) = ABS( D( N ) )
*
* Early return if SIGMX is zero (matrix is already diagonal).
*
IF( SIGMX.EQ.ZERO ) THEN
CALL DLASRT( 'D', N, D, IINFO )
RETURN
END IF
*
DO 20 I = 1, N
SIGMX = MAX( SIGMX, D( I ) )
20 CONTINUE
*
* Copy D and E into WORK (in the Z format) and scale (squaring the
* input data makes scaling by a power of the radix pointless).
*
EPS = DLAMCH( 'Precision' )
SAFMIN = DLAMCH( 'Safe minimum' )
SCALE = SQRT( EPS / SAFMIN )
CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
$ IINFO )
*
* Compute the q's and e's.
*
DO 30 I = 1, 2*N - 1
WORK( I ) = WORK( I )**2
30 CONTINUE
WORK( 2*N ) = ZERO
*
CALL DLASQ2( N, WORK, INFO )
*
IF( INFO.EQ.0 ) THEN
DO 40 I = 1, N
D( I ) = SQRT( WORK( I ) )
40 CONTINUE
CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
END IF
*
RETURN
*
* End of DLASQ1
*
END