diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/dlasd7.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/dlasd7.f')
-rw-r--r-- | SRC/dlasd7.f | 444 |
1 files changed, 444 insertions, 0 deletions
diff --git a/SRC/dlasd7.f b/SRC/dlasd7.f new file mode 100644 index 00000000..27547aaa --- /dev/null +++ b/SRC/dlasd7.f @@ -0,0 +1,444 @@ + SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, + $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, + $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, + $ C, S, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, + $ NR, SQRE + DOUBLE PRECISION ALPHA, BETA, C, S +* .. +* .. Array Arguments .. + INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ), + $ IDXQ( * ), PERM( * ) + DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ), + $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ), + $ ZW( * ) +* .. +* +* Purpose +* ======= +* +* DLASD7 merges the two sets of singular values together into a single +* sorted set. Then it tries to deflate the size of the problem. There +* are two ways in which deflation can occur: when two or more singular +* values are close together or if there is a tiny entry in the Z +* vector. For each such occurrence the order of the related +* secular equation problem is reduced by one. +* +* DLASD7 is called from DLASD6. +* +* Arguments +* ========= +* +* ICOMPQ (input) INTEGER +* Specifies whether singular vectors are to be computed +* in compact form, as follows: +* = 0: Compute singular values only. +* = 1: Compute singular vectors of upper +* bidiagonal matrix in compact form. +* +* NL (input) INTEGER +* The row dimension of the upper block. NL >= 1. +* +* NR (input) INTEGER +* The row dimension of the lower block. NR >= 1. +* +* SQRE (input) INTEGER +* = 0: the lower block is an NR-by-NR square matrix. +* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. +* +* The bidiagonal matrix has +* N = NL + NR + 1 rows and +* M = N + SQRE >= N columns. +* +* K (output) INTEGER +* Contains the dimension of the non-deflated matrix, this is +* the order of the related secular equation. 1 <= K <=N. +* +* D (input/output) DOUBLE PRECISION array, dimension ( N ) +* On entry D contains the singular values of the two submatrices +* to be combined. On exit D contains the trailing (N-K) updated +* singular values (those which were deflated) sorted into +* increasing order. +* +* Z (output) DOUBLE PRECISION array, dimension ( M ) +* On exit Z contains the updating row vector in the secular +* equation. +* +* ZW (workspace) DOUBLE PRECISION array, dimension ( M ) +* Workspace for Z. +* +* VF (input/output) DOUBLE PRECISION array, dimension ( M ) +* On entry, VF(1:NL+1) contains the first components of all +* right singular vectors of the upper block; and VF(NL+2:M) +* contains the first components of all right singular vectors +* of the lower block. On exit, VF contains the first components +* of all right singular vectors of the bidiagonal matrix. +* +* VFW (workspace) DOUBLE PRECISION array, dimension ( M ) +* Workspace for VF. +* +* VL (input/output) DOUBLE PRECISION array, dimension ( M ) +* On entry, VL(1:NL+1) contains the last components of all +* right singular vectors of the upper block; and VL(NL+2:M) +* contains the last components of all right singular vectors +* of the lower block. On exit, VL contains the last components +* of all right singular vectors of the bidiagonal matrix. +* +* VLW (workspace) DOUBLE PRECISION array, dimension ( M ) +* Workspace for VL. +* +* ALPHA (input) DOUBLE PRECISION +* Contains the diagonal element associated with the added row. +* +* BETA (input) DOUBLE PRECISION +* Contains the off-diagonal element associated with the added +* row. +* +* DSIGMA (output) DOUBLE PRECISION array, dimension ( N ) +* Contains a copy of the diagonal elements (K-1 singular values +* and one zero) in the secular equation. +* +* IDX (workspace) INTEGER array, dimension ( N ) +* This will contain the permutation used to sort the contents of +* D into ascending order. +* +* IDXP (workspace) INTEGER array, dimension ( N ) +* This will contain the permutation used to place deflated +* values of D at the end of the array. On output IDXP(2:K) +* points to the nondeflated D-values and IDXP(K+1:N) +* points to the deflated singular values. +* +* IDXQ (input) INTEGER array, dimension ( N ) +* This contains the permutation which separately sorts the two +* sub-problems in D into ascending order. Note that entries in +* the first half of this permutation must first be moved one +* position backward; and entries in the second half +* must first have NL+1 added to their values. +* +* PERM (output) INTEGER array, dimension ( N ) +* The permutations (from deflation and sorting) to be applied +* to each singular block. Not referenced if ICOMPQ = 0. +* +* GIVPTR (output) INTEGER +* The number of Givens rotations which took place in this +* subproblem. Not referenced if ICOMPQ = 0. +* +* GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) +* Each pair of numbers indicates a pair of columns to take place +* in a Givens rotation. Not referenced if ICOMPQ = 0. +* +* LDGCOL (input) INTEGER +* The leading dimension of GIVCOL, must be at least N. +* +* GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) +* Each number indicates the C or S value to be used in the +* corresponding Givens rotation. Not referenced if ICOMPQ = 0. +* +* LDGNUM (input) INTEGER +* The leading dimension of GIVNUM, must be at least N. +* +* C (output) DOUBLE PRECISION +* C contains garbage if SQRE =0 and the C-value of a Givens +* rotation related to the right null space if SQRE = 1. +* +* S (output) DOUBLE PRECISION +* S contains garbage if SQRE =0 and the S-value of a Givens +* rotation related to the right null space if SQRE = 1. +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* Based on contributions by +* Ming Gu and Huan Ren, Computer Science Division, University of +* California at Berkeley, USA +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO, EIGHT + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, + $ EIGHT = 8.0D+0 ) +* .. +* .. Local Scalars .. +* + INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N, + $ NLP1, NLP2 + DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1 +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DLAMRG, DROT, XERBLA +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2 + EXTERNAL DLAMCH, DLAPY2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + N = NL + NR + 1 + M = N + SQRE +* + IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN + INFO = -1 + ELSE IF( NL.LT.1 ) THEN + INFO = -2 + ELSE IF( NR.LT.1 ) THEN + INFO = -3 + ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN + INFO = -4 + ELSE IF( LDGCOL.LT.N ) THEN + INFO = -22 + ELSE IF( LDGNUM.LT.N ) THEN + INFO = -24 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLASD7', -INFO ) + RETURN + END IF +* + NLP1 = NL + 1 + NLP2 = NL + 2 + IF( ICOMPQ.EQ.1 ) THEN + GIVPTR = 0 + END IF +* +* Generate the first part of the vector Z and move the singular +* values in the first part of D one position backward. +* + Z1 = ALPHA*VL( NLP1 ) + VL( NLP1 ) = ZERO + TAU = VF( NLP1 ) + DO 10 I = NL, 1, -1 + Z( I+1 ) = ALPHA*VL( I ) + VL( I ) = ZERO + VF( I+1 ) = VF( I ) + D( I+1 ) = D( I ) + IDXQ( I+1 ) = IDXQ( I ) + 1 + 10 CONTINUE + VF( 1 ) = TAU +* +* Generate the second part of the vector Z. +* + DO 20 I = NLP2, M + Z( I ) = BETA*VF( I ) + VF( I ) = ZERO + 20 CONTINUE +* +* Sort the singular values into increasing order +* + DO 30 I = NLP2, N + IDXQ( I ) = IDXQ( I ) + NLP1 + 30 CONTINUE +* +* DSIGMA, IDXC, IDXC, and ZW are used as storage space. +* + DO 40 I = 2, N + DSIGMA( I ) = D( IDXQ( I ) ) + ZW( I ) = Z( IDXQ( I ) ) + VFW( I ) = VF( IDXQ( I ) ) + VLW( I ) = VL( IDXQ( I ) ) + 40 CONTINUE +* + CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) ) +* + DO 50 I = 2, N + IDXI = 1 + IDX( I ) + D( I ) = DSIGMA( IDXI ) + Z( I ) = ZW( IDXI ) + VF( I ) = VFW( IDXI ) + VL( I ) = VLW( IDXI ) + 50 CONTINUE +* +* Calculate the allowable deflation tolerence +* + EPS = DLAMCH( 'Epsilon' ) + TOL = MAX( ABS( ALPHA ), ABS( BETA ) ) + TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL ) +* +* There are 2 kinds of deflation -- first a value in the z-vector +* is small, second two (or more) singular values are very close +* together (their difference is small). +* +* If the value in the z-vector is small, we simply permute the +* array so that the corresponding singular value is moved to the +* end. +* +* If two values in the D-vector are close, we perform a two-sided +* rotation designed to make one of the corresponding z-vector +* entries zero, and then permute the array so that the deflated +* singular value is moved to the end. +* +* If there are multiple singular values then the problem deflates. +* Here the number of equal singular values are found. As each equal +* singular value is found, an elementary reflector is computed to +* rotate the corresponding singular subspace so that the +* corresponding components of Z are zero in this new basis. +* + K = 1 + K2 = N + 1 + DO 60 J = 2, N + IF( ABS( Z( J ) ).LE.TOL ) THEN +* +* Deflate due to small z component. +* + K2 = K2 - 1 + IDXP( K2 ) = J + IF( J.EQ.N ) + $ GO TO 100 + ELSE + JPREV = J + GO TO 70 + END IF + 60 CONTINUE + 70 CONTINUE + J = JPREV + 80 CONTINUE + J = J + 1 + IF( J.GT.N ) + $ GO TO 90 + IF( ABS( Z( J ) ).LE.TOL ) THEN +* +* Deflate due to small z component. +* + K2 = K2 - 1 + IDXP( K2 ) = J + ELSE +* +* Check if singular values are close enough to allow deflation. +* + IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN +* +* Deflation is possible. +* + S = Z( JPREV ) + C = Z( J ) +* +* Find sqrt(a**2+b**2) without overflow or +* destructive underflow. +* + TAU = DLAPY2( C, S ) + Z( J ) = TAU + Z( JPREV ) = ZERO + C = C / TAU + S = -S / TAU +* +* Record the appropriate Givens rotation +* + IF( ICOMPQ.EQ.1 ) THEN + GIVPTR = GIVPTR + 1 + IDXJP = IDXQ( IDX( JPREV )+1 ) + IDXJ = IDXQ( IDX( J )+1 ) + IF( IDXJP.LE.NLP1 ) THEN + IDXJP = IDXJP - 1 + END IF + IF( IDXJ.LE.NLP1 ) THEN + IDXJ = IDXJ - 1 + END IF + GIVCOL( GIVPTR, 2 ) = IDXJP + GIVCOL( GIVPTR, 1 ) = IDXJ + GIVNUM( GIVPTR, 2 ) = C + GIVNUM( GIVPTR, 1 ) = S + END IF + CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S ) + CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S ) + K2 = K2 - 1 + IDXP( K2 ) = JPREV + JPREV = J + ELSE + K = K + 1 + ZW( K ) = Z( JPREV ) + DSIGMA( K ) = D( JPREV ) + IDXP( K ) = JPREV + JPREV = J + END IF + END IF + GO TO 80 + 90 CONTINUE +* +* Record the last singular value. +* + K = K + 1 + ZW( K ) = Z( JPREV ) + DSIGMA( K ) = D( JPREV ) + IDXP( K ) = JPREV +* + 100 CONTINUE +* +* Sort the singular values into DSIGMA. The singular values which +* were not deflated go into the first K slots of DSIGMA, except +* that DSIGMA(1) is treated separately. +* + DO 110 J = 2, N + JP = IDXP( J ) + DSIGMA( J ) = D( JP ) + VFW( J ) = VF( JP ) + VLW( J ) = VL( JP ) + 110 CONTINUE + IF( ICOMPQ.EQ.1 ) THEN + DO 120 J = 2, N + JP = IDXP( J ) + PERM( J ) = IDXQ( IDX( JP )+1 ) + IF( PERM( J ).LE.NLP1 ) THEN + PERM( J ) = PERM( J ) - 1 + END IF + 120 CONTINUE + END IF +* +* The deflated singular values go back into the last N - K slots of +* D. +* + CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 ) +* +* Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and +* VL(M). +* + DSIGMA( 1 ) = ZERO + HLFTOL = TOL / TWO + IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL ) + $ DSIGMA( 2 ) = HLFTOL + IF( M.GT.N ) THEN + Z( 1 ) = DLAPY2( Z1, Z( M ) ) + IF( Z( 1 ).LE.TOL ) THEN + C = ONE + S = ZERO + Z( 1 ) = TOL + ELSE + C = Z1 / Z( 1 ) + S = -Z( M ) / Z( 1 ) + END IF + CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S ) + CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S ) + ELSE + IF( ABS( Z1 ).LE.TOL ) THEN + Z( 1 ) = TOL + ELSE + Z( 1 ) = Z1 + END IF + END IF +* +* Restore Z, VF, and VL. +* + CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 ) + CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 ) + CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 ) +* + RETURN +* +* End of DLASD7 +* + END |