149 lines
4.3 KiB
FortranFixed
149 lines
4.3 KiB
FortranFixed
SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. 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
|