diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/dlasq1.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/dlasq1.f | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dlasq1.f b/2.3-1/src/fortran/lapack/dlasq1.f new file mode 100644 index 00000000..6f4c3413 --- /dev/null +++ b/2.3-1/src/fortran/lapack/dlasq1.f @@ -0,0 +1,148 @@ + 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 |