SUBROUTINE DLASQ1( N, D, E, WORK, INFO ) * * -- LAPACK routine (instrumented to count ops, version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1999 * * .. Scalar Arguments .. INTEGER INFO, N * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ), WORK( * ) * .. * .. Common block to return operation count .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. DOUBLE PRECISION ITCNT, OPS * .. * * 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 DLAS2, DLASQ2, DLASRT, XERBLA * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, 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). * OPS = OPS + DBLE( 1 + 2*N ) 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. * OPS = OPS + DBLE( 2*N-1 ) 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 OPS = OPS + DBLE( 2*N ) 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