LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasd8.f
Go to the documentation of this file.
00001 *> \brief \b DLASD8
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASD8 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd8.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd8.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd8.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
00022 *                          DSIGMA, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            ICOMPQ, INFO, K, LDDIFR
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( LDDIFR, * ),
00029 *      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
00030 *      $                   Z( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLASD8 finds the square roots of the roots of the secular equation,
00040 *> as defined by the values in DSIGMA and Z. It makes the appropriate
00041 *> calls to DLASD4, and stores, for each  element in D, the distance
00042 *> to its two nearest poles (elements in DSIGMA). It also updates
00043 *> the arrays VF and VL, the first and last components of all the
00044 *> right singular vectors of the original bidiagonal matrix.
00045 *>
00046 *> DLASD8 is called from DLASD6.
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] ICOMPQ
00053 *> \verbatim
00054 *>          ICOMPQ is INTEGER
00055 *>          Specifies whether singular vectors are to be computed in
00056 *>          factored form in the calling routine:
00057 *>          = 0: Compute singular values only.
00058 *>          = 1: Compute singular vectors in factored form as well.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] K
00062 *> \verbatim
00063 *>          K is INTEGER
00064 *>          The number of terms in the rational function to be solved
00065 *>          by DLASD4.  K >= 1.
00066 *> \endverbatim
00067 *>
00068 *> \param[out] D
00069 *> \verbatim
00070 *>          D is DOUBLE PRECISION array, dimension ( K )
00071 *>          On output, D contains the updated singular values.
00072 *> \endverbatim
00073 *>
00074 *> \param[in,out] Z
00075 *> \verbatim
00076 *>          Z is DOUBLE PRECISION array, dimension ( K )
00077 *>          On entry, the first K elements of this array contain the
00078 *>          components of the deflation-adjusted updating row vector.
00079 *>          On exit, Z is updated.
00080 *> \endverbatim
00081 *>
00082 *> \param[in,out] VF
00083 *> \verbatim
00084 *>          VF is DOUBLE PRECISION array, dimension ( K )
00085 *>          On entry, VF contains  information passed through DBEDE8.
00086 *>          On exit, VF contains the first K components of the first
00087 *>          components of all right singular vectors of the bidiagonal
00088 *>          matrix.
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] VL
00092 *> \verbatim
00093 *>          VL is DOUBLE PRECISION array, dimension ( K )
00094 *>          On entry, VL contains  information passed through DBEDE8.
00095 *>          On exit, VL contains the first K components of the last
00096 *>          components of all right singular vectors of the bidiagonal
00097 *>          matrix.
00098 *> \endverbatim
00099 *>
00100 *> \param[out] DIFL
00101 *> \verbatim
00102 *>          DIFL is DOUBLE PRECISION array, dimension ( K )
00103 *>          On exit, DIFL(I) = D(I) - DSIGMA(I).
00104 *> \endverbatim
00105 *>
00106 *> \param[out] DIFR
00107 *> \verbatim
00108 *>          DIFR is DOUBLE PRECISION array,
00109 *>                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
00110 *>                   dimension ( K ) if ICOMPQ = 0.
00111 *>          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
00112 *>          defined and will not be referenced.
00113 *>
00114 *>          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
00115 *>          normalizing factors for the right singular vector matrix.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] LDDIFR
00119 *> \verbatim
00120 *>          LDDIFR is INTEGER
00121 *>          The leading dimension of DIFR, must be at least K.
00122 *> \endverbatim
00123 *>
00124 *> \param[in,out] DSIGMA
00125 *> \verbatim
00126 *>          DSIGMA is DOUBLE PRECISION array, dimension ( K )
00127 *>          On entry, the first K elements of this array contain the old
00128 *>          roots of the deflated updating problem.  These are the poles
00129 *>          of the secular equation.
00130 *>          On exit, the elements of DSIGMA may be very slightly altered
00131 *>          in value.
00132 *> \endverbatim
00133 *>
00134 *> \param[out] WORK
00135 *> \verbatim
00136 *>          WORK is DOUBLE PRECISION array, dimension at least 3 * K
00137 *> \endverbatim
00138 *>
00139 *> \param[out] INFO
00140 *> \verbatim
00141 *>          INFO is INTEGER
00142 *>          = 0:  successful exit.
00143 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00144 *>          > 0:  if INFO = 1, a singular value did not converge
00145 *> \endverbatim
00146 *
00147 *  Authors:
00148 *  ========
00149 *
00150 *> \author Univ. of Tennessee 
00151 *> \author Univ. of California Berkeley 
00152 *> \author Univ. of Colorado Denver 
00153 *> \author NAG Ltd. 
00154 *
00155 *> \date November 2011
00156 *
00157 *> \ingroup auxOTHERauxiliary
00158 *
00159 *> \par Contributors:
00160 *  ==================
00161 *>
00162 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00163 *>     California at Berkeley, USA
00164 *>
00165 *  =====================================================================
00166       SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
00167      $                   DSIGMA, WORK, INFO )
00168 *
00169 *  -- LAPACK auxiliary routine (version 3.4.0) --
00170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00172 *     November 2011
00173 *
00174 *     .. Scalar Arguments ..
00175       INTEGER            ICOMPQ, INFO, K, LDDIFR
00176 *     ..
00177 *     .. Array Arguments ..
00178       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( LDDIFR, * ),
00179      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
00180      $                   Z( * )
00181 *     ..
00182 *
00183 *  =====================================================================
00184 *
00185 *     .. Parameters ..
00186       DOUBLE PRECISION   ONE
00187       PARAMETER          ( ONE = 1.0D+0 )
00188 *     ..
00189 *     .. Local Scalars ..
00190       INTEGER            I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
00191       DOUBLE PRECISION   DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
00192 *     ..
00193 *     .. External Subroutines ..
00194       EXTERNAL           DCOPY, DLASCL, DLASD4, DLASET, XERBLA
00195 *     ..
00196 *     .. External Functions ..
00197       DOUBLE PRECISION   DDOT, DLAMC3, DNRM2
00198       EXTERNAL           DDOT, DLAMC3, DNRM2
00199 *     ..
00200 *     .. Intrinsic Functions ..
00201       INTRINSIC          ABS, SIGN, SQRT
00202 *     ..
00203 *     .. Executable Statements ..
00204 *
00205 *     Test the input parameters.
00206 *
00207       INFO = 0
00208 *
00209       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00210          INFO = -1
00211       ELSE IF( K.LT.1 ) THEN
00212          INFO = -2
00213       ELSE IF( LDDIFR.LT.K ) THEN
00214          INFO = -9
00215       END IF
00216       IF( INFO.NE.0 ) THEN
00217          CALL XERBLA( 'DLASD8', -INFO )
00218          RETURN
00219       END IF
00220 *
00221 *     Quick return if possible
00222 *
00223       IF( K.EQ.1 ) THEN
00224          D( 1 ) = ABS( Z( 1 ) )
00225          DIFL( 1 ) = D( 1 )
00226          IF( ICOMPQ.EQ.1 ) THEN
00227             DIFL( 2 ) = ONE
00228             DIFR( 1, 2 ) = ONE
00229          END IF
00230          RETURN
00231       END IF
00232 *
00233 *     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
00234 *     be computed with high relative accuracy (barring over/underflow).
00235 *     This is a problem on machines without a guard digit in
00236 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00237 *     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
00238 *     which on any of these machines zeros out the bottommost
00239 *     bit of DSIGMA(I) if it is 1; this makes the subsequent
00240 *     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
00241 *     occurs. On binary machines with a guard digit (almost all
00242 *     machines) it does not change DSIGMA(I) at all. On hexadecimal
00243 *     and decimal machines with a guard digit, it slightly
00244 *     changes the bottommost bits of DSIGMA(I). It does not account
00245 *     for hexadecimal or decimal machines without guard digits
00246 *     (we know of none). We use a subroutine call to compute
00247 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00248 *     this code.
00249 *
00250       DO 10 I = 1, K
00251          DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
00252    10 CONTINUE
00253 *
00254 *     Book keeping.
00255 *
00256       IWK1 = 1
00257       IWK2 = IWK1 + K
00258       IWK3 = IWK2 + K
00259       IWK2I = IWK2 - 1
00260       IWK3I = IWK3 - 1
00261 *
00262 *     Normalize Z.
00263 *
00264       RHO = DNRM2( K, Z, 1 )
00265       CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
00266       RHO = RHO*RHO
00267 *
00268 *     Initialize WORK(IWK3).
00269 *
00270       CALL DLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
00271 *
00272 *     Compute the updated singular values, the arrays DIFL, DIFR,
00273 *     and the updated Z.
00274 *
00275       DO 40 J = 1, K
00276          CALL DLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
00277      $                WORK( IWK2 ), INFO )
00278 *
00279 *        If the root finder fails, the computation is terminated.
00280 *
00281          IF( INFO.NE.0 ) THEN
00282             CALL XERBLA( 'DLASD4', -INFO )
00283             RETURN
00284          END IF
00285          WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
00286          DIFL( J ) = -WORK( J )
00287          DIFR( J, 1 ) = -WORK( J+1 )
00288          DO 20 I = 1, J - 1
00289             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
00290      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
00291      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
00292      $                        DSIGMA( J ) )
00293    20    CONTINUE
00294          DO 30 I = J + 1, K
00295             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
00296      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
00297      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
00298      $                        DSIGMA( J ) )
00299    30    CONTINUE
00300    40 CONTINUE
00301 *
00302 *     Compute updated Z.
00303 *
00304       DO 50 I = 1, K
00305          Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
00306    50 CONTINUE
00307 *
00308 *     Update VF and VL.
00309 *
00310       DO 80 J = 1, K
00311          DIFLJ = DIFL( J )
00312          DJ = D( J )
00313          DSIGJ = -DSIGMA( J )
00314          IF( J.LT.K ) THEN
00315             DIFRJ = -DIFR( J, 1 )
00316             DSIGJP = -DSIGMA( J+1 )
00317          END IF
00318          WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
00319          DO 60 I = 1, J - 1
00320             WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
00321      $                   / ( DSIGMA( I )+DJ )
00322    60    CONTINUE
00323          DO 70 I = J + 1, K
00324             WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
00325      $                   / ( DSIGMA( I )+DJ )
00326    70    CONTINUE
00327          TEMP = DNRM2( K, WORK, 1 )
00328          WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP
00329          WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP
00330          IF( ICOMPQ.EQ.1 ) THEN
00331             DIFR( J, 2 ) = TEMP
00332          END IF
00333    80 CONTINUE
00334 *
00335       CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 )
00336       CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 )
00337 *
00338       RETURN
00339 *
00340 *     End of DLASD8
00341 *
00342       END
00343 
 All Files Functions