LAPACK 3.3.0

dlasd7.f

Go to the documentation of this file.
00001       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
00002      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
00003      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
00004      $                   C, S, INFO )
00005 *
00006 *  -- LAPACK auxiliary routine (version 3.2) --
00007 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00008 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00009 *     November 2006
00010 *
00011 *     .. Scalar Arguments ..
00012       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
00013      $                   NR, SQRE
00014       DOUBLE PRECISION   ALPHA, BETA, C, S
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
00018      $                   IDXQ( * ), PERM( * )
00019       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
00020      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
00021      $                   ZW( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 *
00027 *  DLASD7 merges the two sets of singular values together into a single
00028 *  sorted set. Then it tries to deflate the size of the problem. There
00029 *  are two ways in which deflation can occur:  when two or more singular
00030 *  values are close together or if there is a tiny entry in the Z
00031 *  vector. For each such occurrence the order of the related
00032 *  secular equation problem is reduced by one.
00033 *
00034 *  DLASD7 is called from DLASD6.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  ICOMPQ  (input) INTEGER
00040 *          Specifies whether singular vectors are to be computed
00041 *          in compact form, as follows:
00042 *          = 0: Compute singular values only.
00043 *          = 1: Compute singular vectors of upper
00044 *               bidiagonal matrix in compact form.
00045 *
00046 *  NL     (input) INTEGER
00047 *         The row dimension of the upper block. NL >= 1.
00048 *
00049 *  NR     (input) INTEGER
00050 *         The row dimension of the lower block. NR >= 1.
00051 *
00052 *  SQRE   (input) INTEGER
00053 *         = 0: the lower block is an NR-by-NR square matrix.
00054 *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
00055 *
00056 *         The bidiagonal matrix has
00057 *         N = NL + NR + 1 rows and
00058 *         M = N + SQRE >= N columns.
00059 *
00060 *  K      (output) INTEGER
00061 *         Contains the dimension of the non-deflated matrix, this is
00062 *         the order of the related secular equation. 1 <= K <=N.
00063 *
00064 *  D      (input/output) DOUBLE PRECISION array, dimension ( N )
00065 *         On entry D contains the singular values of the two submatrices
00066 *         to be combined. On exit D contains the trailing (N-K) updated
00067 *         singular values (those which were deflated) sorted into
00068 *         increasing order.
00069 *
00070 *  Z      (output) DOUBLE PRECISION array, dimension ( M )
00071 *         On exit Z contains the updating row vector in the secular
00072 *         equation.
00073 *
00074 *  ZW     (workspace) DOUBLE PRECISION array, dimension ( M )
00075 *         Workspace for Z.
00076 *
00077 *  VF     (input/output) DOUBLE PRECISION array, dimension ( M )
00078 *         On entry, VF(1:NL+1) contains the first components of all
00079 *         right singular vectors of the upper block; and VF(NL+2:M)
00080 *         contains the first components of all right singular vectors
00081 *         of the lower block. On exit, VF contains the first components
00082 *         of all right singular vectors of the bidiagonal matrix.
00083 *
00084 *  VFW    (workspace) DOUBLE PRECISION array, dimension ( M )
00085 *         Workspace for VF.
00086 *
00087 *  VL     (input/output) DOUBLE PRECISION array, dimension ( M )
00088 *         On entry, VL(1:NL+1) contains the  last components of all
00089 *         right singular vectors of the upper block; and VL(NL+2:M)
00090 *         contains the last components of all right singular vectors
00091 *         of the lower block. On exit, VL contains the last components
00092 *         of all right singular vectors of the bidiagonal matrix.
00093 *
00094 *  VLW    (workspace) DOUBLE PRECISION array, dimension ( M )
00095 *         Workspace for VL.
00096 *
00097 *  ALPHA  (input) DOUBLE PRECISION
00098 *         Contains the diagonal element associated with the added row.
00099 *
00100 *  BETA   (input) DOUBLE PRECISION
00101 *         Contains the off-diagonal element associated with the added
00102 *         row.
00103 *
00104 *  DSIGMA (output) DOUBLE PRECISION array, dimension ( N )
00105 *         Contains a copy of the diagonal elements (K-1 singular values
00106 *         and one zero) in the secular equation.
00107 *
00108 *  IDX    (workspace) INTEGER array, dimension ( N )
00109 *         This will contain the permutation used to sort the contents of
00110 *         D into ascending order.
00111 *
00112 *  IDXP   (workspace) INTEGER array, dimension ( N )
00113 *         This will contain the permutation used to place deflated
00114 *         values of D at the end of the array. On output IDXP(2:K)
00115 *         points to the nondeflated D-values and IDXP(K+1:N)
00116 *         points to the deflated singular values.
00117 *
00118 *  IDXQ   (input) INTEGER array, dimension ( N )
00119 *         This contains the permutation which separately sorts the two
00120 *         sub-problems in D into ascending order.  Note that entries in
00121 *         the first half of this permutation must first be moved one
00122 *         position backward; and entries in the second half
00123 *         must first have NL+1 added to their values.
00124 *
00125 *  PERM   (output) INTEGER array, dimension ( N )
00126 *         The permutations (from deflation and sorting) to be applied
00127 *         to each singular block. Not referenced if ICOMPQ = 0.
00128 *
00129 *  GIVPTR (output) INTEGER
00130 *         The number of Givens rotations which took place in this
00131 *         subproblem. Not referenced if ICOMPQ = 0.
00132 *
00133 *  GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 )
00134 *         Each pair of numbers indicates a pair of columns to take place
00135 *         in a Givens rotation. Not referenced if ICOMPQ = 0.
00136 *
00137 *  LDGCOL (input) INTEGER
00138 *         The leading dimension of GIVCOL, must be at least N.
00139 *
00140 *  GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00141 *         Each number indicates the C or S value to be used in the
00142 *         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
00143 *
00144 *  LDGNUM (input) INTEGER
00145 *         The leading dimension of GIVNUM, must be at least N.
00146 *
00147 *  C      (output) DOUBLE PRECISION
00148 *         C contains garbage if SQRE =0 and the C-value of a Givens
00149 *         rotation related to the right null space if SQRE = 1.
00150 *
00151 *  S      (output) DOUBLE PRECISION
00152 *         S contains garbage if SQRE =0 and the S-value of a Givens
00153 *         rotation related to the right null space if SQRE = 1.
00154 *
00155 *  INFO   (output) INTEGER
00156 *         = 0:  successful exit.
00157 *         < 0:  if INFO = -i, the i-th argument had an illegal value.
00158 *
00159 *  Further Details
00160 *  ===============
00161 *
00162 *  Based on contributions by
00163 *     Ming Gu and Huan Ren, Computer Science Division, University of
00164 *     California at Berkeley, USA
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
00170       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
00171      $                   EIGHT = 8.0D+0 )
00172 *     ..
00173 *     .. Local Scalars ..
00174 *
00175       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
00176      $                   NLP1, NLP2
00177       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
00178 *     ..
00179 *     .. External Subroutines ..
00180       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
00181 *     ..
00182 *     .. External Functions ..
00183       DOUBLE PRECISION   DLAMCH, DLAPY2
00184       EXTERNAL           DLAMCH, DLAPY2
00185 *     ..
00186 *     .. Intrinsic Functions ..
00187       INTRINSIC          ABS, MAX
00188 *     ..
00189 *     .. Executable Statements ..
00190 *
00191 *     Test the input parameters.
00192 *
00193       INFO = 0
00194       N = NL + NR + 1
00195       M = N + SQRE
00196 *
00197       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00198          INFO = -1
00199       ELSE IF( NL.LT.1 ) THEN
00200          INFO = -2
00201       ELSE IF( NR.LT.1 ) THEN
00202          INFO = -3
00203       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00204          INFO = -4
00205       ELSE IF( LDGCOL.LT.N ) THEN
00206          INFO = -22
00207       ELSE IF( LDGNUM.LT.N ) THEN
00208          INFO = -24
00209       END IF
00210       IF( INFO.NE.0 ) THEN
00211          CALL XERBLA( 'DLASD7', -INFO )
00212          RETURN
00213       END IF
00214 *
00215       NLP1 = NL + 1
00216       NLP2 = NL + 2
00217       IF( ICOMPQ.EQ.1 ) THEN
00218          GIVPTR = 0
00219       END IF
00220 *
00221 *     Generate the first part of the vector Z and move the singular
00222 *     values in the first part of D one position backward.
00223 *
00224       Z1 = ALPHA*VL( NLP1 )
00225       VL( NLP1 ) = ZERO
00226       TAU = VF( NLP1 )
00227       DO 10 I = NL, 1, -1
00228          Z( I+1 ) = ALPHA*VL( I )
00229          VL( I ) = ZERO
00230          VF( I+1 ) = VF( I )
00231          D( I+1 ) = D( I )
00232          IDXQ( I+1 ) = IDXQ( I ) + 1
00233    10 CONTINUE
00234       VF( 1 ) = TAU
00235 *
00236 *     Generate the second part of the vector Z.
00237 *
00238       DO 20 I = NLP2, M
00239          Z( I ) = BETA*VF( I )
00240          VF( I ) = ZERO
00241    20 CONTINUE
00242 *
00243 *     Sort the singular values into increasing order
00244 *
00245       DO 30 I = NLP2, N
00246          IDXQ( I ) = IDXQ( I ) + NLP1
00247    30 CONTINUE
00248 *
00249 *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
00250 *
00251       DO 40 I = 2, N
00252          DSIGMA( I ) = D( IDXQ( I ) )
00253          ZW( I ) = Z( IDXQ( I ) )
00254          VFW( I ) = VF( IDXQ( I ) )
00255          VLW( I ) = VL( IDXQ( I ) )
00256    40 CONTINUE
00257 *
00258       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
00259 *
00260       DO 50 I = 2, N
00261          IDXI = 1 + IDX( I )
00262          D( I ) = DSIGMA( IDXI )
00263          Z( I ) = ZW( IDXI )
00264          VF( I ) = VFW( IDXI )
00265          VL( I ) = VLW( IDXI )
00266    50 CONTINUE
00267 *
00268 *     Calculate the allowable deflation tolerence
00269 *
00270       EPS = DLAMCH( 'Epsilon' )
00271       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
00272       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
00273 *
00274 *     There are 2 kinds of deflation -- first a value in the z-vector
00275 *     is small, second two (or more) singular values are very close
00276 *     together (their difference is small).
00277 *
00278 *     If the value in the z-vector is small, we simply permute the
00279 *     array so that the corresponding singular value is moved to the
00280 *     end.
00281 *
00282 *     If two values in the D-vector are close, we perform a two-sided
00283 *     rotation designed to make one of the corresponding z-vector
00284 *     entries zero, and then permute the array so that the deflated
00285 *     singular value is moved to the end.
00286 *
00287 *     If there are multiple singular values then the problem deflates.
00288 *     Here the number of equal singular values are found.  As each equal
00289 *     singular value is found, an elementary reflector is computed to
00290 *     rotate the corresponding singular subspace so that the
00291 *     corresponding components of Z are zero in this new basis.
00292 *
00293       K = 1
00294       K2 = N + 1
00295       DO 60 J = 2, N
00296          IF( ABS( Z( J ) ).LE.TOL ) THEN
00297 *
00298 *           Deflate due to small z component.
00299 *
00300             K2 = K2 - 1
00301             IDXP( K2 ) = J
00302             IF( J.EQ.N )
00303      $         GO TO 100
00304          ELSE
00305             JPREV = J
00306             GO TO 70
00307          END IF
00308    60 CONTINUE
00309    70 CONTINUE
00310       J = JPREV
00311    80 CONTINUE
00312       J = J + 1
00313       IF( J.GT.N )
00314      $   GO TO 90
00315       IF( ABS( Z( J ) ).LE.TOL ) THEN
00316 *
00317 *        Deflate due to small z component.
00318 *
00319          K2 = K2 - 1
00320          IDXP( K2 ) = J
00321       ELSE
00322 *
00323 *        Check if singular values are close enough to allow deflation.
00324 *
00325          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
00326 *
00327 *           Deflation is possible.
00328 *
00329             S = Z( JPREV )
00330             C = Z( J )
00331 *
00332 *           Find sqrt(a**2+b**2) without overflow or
00333 *           destructive underflow.
00334 *
00335             TAU = DLAPY2( C, S )
00336             Z( J ) = TAU
00337             Z( JPREV ) = ZERO
00338             C = C / TAU
00339             S = -S / TAU
00340 *
00341 *           Record the appropriate Givens rotation
00342 *
00343             IF( ICOMPQ.EQ.1 ) THEN
00344                GIVPTR = GIVPTR + 1
00345                IDXJP = IDXQ( IDX( JPREV )+1 )
00346                IDXJ = IDXQ( IDX( J )+1 )
00347                IF( IDXJP.LE.NLP1 ) THEN
00348                   IDXJP = IDXJP - 1
00349                END IF
00350                IF( IDXJ.LE.NLP1 ) THEN
00351                   IDXJ = IDXJ - 1
00352                END IF
00353                GIVCOL( GIVPTR, 2 ) = IDXJP
00354                GIVCOL( GIVPTR, 1 ) = IDXJ
00355                GIVNUM( GIVPTR, 2 ) = C
00356                GIVNUM( GIVPTR, 1 ) = S
00357             END IF
00358             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
00359             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
00360             K2 = K2 - 1
00361             IDXP( K2 ) = JPREV
00362             JPREV = J
00363          ELSE
00364             K = K + 1
00365             ZW( K ) = Z( JPREV )
00366             DSIGMA( K ) = D( JPREV )
00367             IDXP( K ) = JPREV
00368             JPREV = J
00369          END IF
00370       END IF
00371       GO TO 80
00372    90 CONTINUE
00373 *
00374 *     Record the last singular value.
00375 *
00376       K = K + 1
00377       ZW( K ) = Z( JPREV )
00378       DSIGMA( K ) = D( JPREV )
00379       IDXP( K ) = JPREV
00380 *
00381   100 CONTINUE
00382 *
00383 *     Sort the singular values into DSIGMA. The singular values which
00384 *     were not deflated go into the first K slots of DSIGMA, except
00385 *     that DSIGMA(1) is treated separately.
00386 *
00387       DO 110 J = 2, N
00388          JP = IDXP( J )
00389          DSIGMA( J ) = D( JP )
00390          VFW( J ) = VF( JP )
00391          VLW( J ) = VL( JP )
00392   110 CONTINUE
00393       IF( ICOMPQ.EQ.1 ) THEN
00394          DO 120 J = 2, N
00395             JP = IDXP( J )
00396             PERM( J ) = IDXQ( IDX( JP )+1 )
00397             IF( PERM( J ).LE.NLP1 ) THEN
00398                PERM( J ) = PERM( J ) - 1
00399             END IF
00400   120    CONTINUE
00401       END IF
00402 *
00403 *     The deflated singular values go back into the last N - K slots of
00404 *     D.
00405 *
00406       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
00407 *
00408 *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
00409 *     VL(M).
00410 *
00411       DSIGMA( 1 ) = ZERO
00412       HLFTOL = TOL / TWO
00413       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
00414      $   DSIGMA( 2 ) = HLFTOL
00415       IF( M.GT.N ) THEN
00416          Z( 1 ) = DLAPY2( Z1, Z( M ) )
00417          IF( Z( 1 ).LE.TOL ) THEN
00418             C = ONE
00419             S = ZERO
00420             Z( 1 ) = TOL
00421          ELSE
00422             C = Z1 / Z( 1 )
00423             S = -Z( M ) / Z( 1 )
00424          END IF
00425          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
00426          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
00427       ELSE
00428          IF( ABS( Z1 ).LE.TOL ) THEN
00429             Z( 1 ) = TOL
00430          ELSE
00431             Z( 1 ) = Z1
00432          END IF
00433       END IF
00434 *
00435 *     Restore Z, VF, and VL.
00436 *
00437       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
00438       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
00439       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
00440 *
00441       RETURN
00442 *
00443 *     End of DLASD7
00444 *
00445       END
 All Files Functions