LAPACK 3.3.1
Linear Algebra PACKage

slasd8.f

Go to the documentation of this file.
00001       SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
00002      $                   DSIGMA, WORK, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.3.0) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            ICOMPQ, INFO, K, LDDIFR
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               D( * ), DIFL( * ), DIFR( LDDIFR, * ),
00014      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
00015      $                   Z( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SLASD8 finds the square roots of the roots of the secular equation,
00022 *  as defined by the values in DSIGMA and Z. It makes the appropriate
00023 *  calls to SLASD4, and stores, for each  element in D, the distance
00024 *  to its two nearest poles (elements in DSIGMA). It also updates
00025 *  the arrays VF and VL, the first and last components of all the
00026 *  right singular vectors of the original bidiagonal matrix.
00027 *
00028 *  SLASD8 is called from SLASD6.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  ICOMPQ  (input) INTEGER
00034 *          Specifies whether singular vectors are to be computed in
00035 *          factored form in the calling routine:
00036 *          = 0: Compute singular values only.
00037 *          = 1: Compute singular vectors in factored form as well.
00038 *
00039 *  K       (input) INTEGER
00040 *          The number of terms in the rational function to be solved
00041 *          by SLASD4.  K >= 1.
00042 *
00043 *  D       (output) REAL array, dimension ( K )
00044 *          On output, D contains the updated singular values.
00045 *
00046 *  Z       (input/output) REAL array, dimension ( K )
00047 *          On entry, the first K elements of this array contain the
00048 *          components of the deflation-adjusted updating row vector.
00049 *          On exit, Z is updated.
00050 *
00051 *  VF      (input/output) REAL array, dimension ( K )
00052 *          On entry, VF contains  information passed through DBEDE8.
00053 *          On exit, VF contains the first K components of the first
00054 *          components of all right singular vectors of the bidiagonal
00055 *          matrix.
00056 *
00057 *  VL      (input/output) REAL array, dimension ( K )
00058 *          On entry, VL contains  information passed through DBEDE8.
00059 *          On exit, VL contains the first K components of the last
00060 *          components of all right singular vectors of the bidiagonal
00061 *          matrix.
00062 *
00063 *  DIFL    (output) REAL array, dimension ( K )
00064 *          On exit, DIFL(I) = D(I) - DSIGMA(I).
00065 *
00066 *  DIFR    (output) REAL array,
00067 *                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
00068 *                   dimension ( K ) if ICOMPQ = 0.
00069 *          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
00070 *          defined and will not be referenced.
00071 *
00072 *          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
00073 *          normalizing factors for the right singular vector matrix.
00074 *
00075 *  LDDIFR  (input) INTEGER
00076 *          The leading dimension of DIFR, must be at least K.
00077 *
00078 *  DSIGMA  (input/output) REAL array, dimension ( K )
00079 *          On entry, the first K elements of this array contain the old
00080 *          roots of the deflated updating problem.  These are the poles
00081 *          of the secular equation.
00082 *          On exit, the elements of DSIGMA may be very slightly altered
00083 *          in value.
00084 *
00085 *  WORK    (workspace) REAL array, dimension at least 3 * K
00086 *
00087 *  INFO    (output) INTEGER
00088 *          = 0:  successful exit.
00089 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00090 *          > 0:  if INFO = 1, a singular value did not converge
00091 *
00092 *  Further Details
00093 *  ===============
00094 *
00095 *  Based on contributions by
00096 *     Ming Gu and Huan Ren, Computer Science Division, University of
00097 *     California at Berkeley, USA
00098 *
00099 *  =====================================================================
00100 *
00101 *     .. Parameters ..
00102       REAL               ONE
00103       PARAMETER          ( ONE = 1.0E+0 )
00104 *     ..
00105 *     .. Local Scalars ..
00106       INTEGER            I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
00107       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
00108 *     ..
00109 *     .. External Subroutines ..
00110       EXTERNAL           SCOPY, SLASCL, SLASD4, SLASET, XERBLA
00111 *     ..
00112 *     .. External Functions ..
00113       REAL               SDOT, SLAMC3, SNRM2
00114       EXTERNAL           SDOT, SLAMC3, SNRM2
00115 *     ..
00116 *     .. Intrinsic Functions ..
00117       INTRINSIC          ABS, SIGN, SQRT
00118 *     ..
00119 *     .. Executable Statements ..
00120 *
00121 *     Test the input parameters.
00122 *
00123       INFO = 0
00124 *
00125       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00126          INFO = -1
00127       ELSE IF( K.LT.1 ) THEN
00128          INFO = -2
00129       ELSE IF( LDDIFR.LT.K ) THEN
00130          INFO = -9
00131       END IF
00132       IF( INFO.NE.0 ) THEN
00133          CALL XERBLA( 'SLASD8', -INFO )
00134          RETURN
00135       END IF
00136 *
00137 *     Quick return if possible
00138 *
00139       IF( K.EQ.1 ) THEN
00140          D( 1 ) = ABS( Z( 1 ) )
00141          DIFL( 1 ) = D( 1 )
00142          IF( ICOMPQ.EQ.1 ) THEN
00143             DIFL( 2 ) = ONE
00144             DIFR( 1, 2 ) = ONE
00145          END IF
00146          RETURN
00147       END IF
00148 *
00149 *     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
00150 *     be computed with high relative accuracy (barring over/underflow).
00151 *     This is a problem on machines without a guard digit in
00152 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00153 *     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
00154 *     which on any of these machines zeros out the bottommost
00155 *     bit of DSIGMA(I) if it is 1; this makes the subsequent
00156 *     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
00157 *     occurs. On binary machines with a guard digit (almost all
00158 *     machines) it does not change DSIGMA(I) at all. On hexadecimal
00159 *     and decimal machines with a guard digit, it slightly
00160 *     changes the bottommost bits of DSIGMA(I). It does not account
00161 *     for hexadecimal or decimal machines without guard digits
00162 *     (we know of none). We use a subroutine call to compute
00163 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00164 *     this code.
00165 *
00166       DO 10 I = 1, K
00167          DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
00168    10 CONTINUE
00169 *
00170 *     Book keeping.
00171 *
00172       IWK1 = 1
00173       IWK2 = IWK1 + K
00174       IWK3 = IWK2 + K
00175       IWK2I = IWK2 - 1
00176       IWK3I = IWK3 - 1
00177 *
00178 *     Normalize Z.
00179 *
00180       RHO = SNRM2( K, Z, 1 )
00181       CALL SLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
00182       RHO = RHO*RHO
00183 *
00184 *     Initialize WORK(IWK3).
00185 *
00186       CALL SLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
00187 *
00188 *     Compute the updated singular values, the arrays DIFL, DIFR,
00189 *     and the updated Z.
00190 *
00191       DO 40 J = 1, K
00192          CALL SLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
00193      $                WORK( IWK2 ), INFO )
00194 *
00195 *        If the root finder fails, the computation is terminated.
00196 *
00197          IF( INFO.NE.0 ) THEN
00198             CALL XERBLA( 'SLASD4', -INFO )
00199             RETURN
00200          END IF
00201          WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
00202          DIFL( J ) = -WORK( J )
00203          DIFR( J, 1 ) = -WORK( J+1 )
00204          DO 20 I = 1, J - 1
00205             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
00206      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
00207      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
00208      $                        DSIGMA( J ) )
00209    20    CONTINUE
00210          DO 30 I = J + 1, K
00211             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
00212      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
00213      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
00214      $                        DSIGMA( J ) )
00215    30    CONTINUE
00216    40 CONTINUE
00217 *
00218 *     Compute updated Z.
00219 *
00220       DO 50 I = 1, K
00221          Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
00222    50 CONTINUE
00223 *
00224 *     Update VF and VL.
00225 *
00226       DO 80 J = 1, K
00227          DIFLJ = DIFL( J )
00228          DJ = D( J )
00229          DSIGJ = -DSIGMA( J )
00230          IF( J.LT.K ) THEN
00231             DIFRJ = -DIFR( J, 1 )
00232             DSIGJP = -DSIGMA( J+1 )
00233          END IF
00234          WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
00235          DO 60 I = 1, J - 1
00236             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
00237      $                   / ( DSIGMA( I )+DJ )
00238    60    CONTINUE
00239          DO 70 I = J + 1, K
00240             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
00241      $                   / ( DSIGMA( I )+DJ )
00242    70    CONTINUE
00243          TEMP = SNRM2( K, WORK, 1 )
00244          WORK( IWK2I+J ) = SDOT( K, WORK, 1, VF, 1 ) / TEMP
00245          WORK( IWK3I+J ) = SDOT( K, WORK, 1, VL, 1 ) / TEMP
00246          IF( ICOMPQ.EQ.1 ) THEN
00247             DIFR( J, 2 ) = TEMP
00248          END IF
00249    80 CONTINUE
00250 *
00251       CALL SCOPY( K, WORK( IWK2 ), 1, VF, 1 )
00252       CALL SCOPY( K, WORK( IWK3 ), 1, VL, 1 )
00253 *
00254       RETURN
00255 *
00256 *     End of SLASD8
00257 *
00258       END
00259 
 All Files Functions