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
CVSweb interface <joel.bertrand@systella.fr>