LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaed8.f
Go to the documentation of this file.
00001 *> \brief \b DLAED8
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAED8 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed8.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed8.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed8.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
00022 *                          CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
00023 *                          GIVCOL, GIVNUM, INDXP, INDX, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
00027 *      $                   QSIZ
00028 *       DOUBLE PRECISION   RHO
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
00032 *      $                   INDXQ( * ), PERM( * )
00033 *       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
00034 *      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> DLAED8 merges the two sets of eigenvalues together into a single
00044 *> sorted set.  Then it tries to deflate the size of the problem.
00045 *> There are two ways in which deflation can occur:  when two or more
00046 *> eigenvalues are close together or if there is a tiny element in the
00047 *> Z vector.  For each such occurrence the order of the related secular
00048 *> equation problem is reduced by one.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] ICOMPQ
00055 *> \verbatim
00056 *>          ICOMPQ is INTEGER
00057 *>          = 0:  Compute eigenvalues only.
00058 *>          = 1:  Compute eigenvectors of original dense symmetric matrix
00059 *>                also.  On entry, Q contains the orthogonal matrix used
00060 *>                to reduce the original matrix to tridiagonal form.
00061 *> \endverbatim
00062 *>
00063 *> \param[out] K
00064 *> \verbatim
00065 *>          K is INTEGER
00066 *>         The number of non-deflated eigenvalues, and the order of the
00067 *>         related secular equation.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] N
00071 *> \verbatim
00072 *>          N is INTEGER
00073 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] QSIZ
00077 *> \verbatim
00078 *>          QSIZ is INTEGER
00079 *>         The dimension of the orthogonal matrix used to reduce
00080 *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00081 *> \endverbatim
00082 *>
00083 *> \param[in,out] D
00084 *> \verbatim
00085 *>          D is DOUBLE PRECISION array, dimension (N)
00086 *>         On entry, the eigenvalues of the two submatrices to be
00087 *>         combined.  On exit, the trailing (N-K) updated eigenvalues
00088 *>         (those which were deflated) sorted into increasing order.
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] Q
00092 *> \verbatim
00093 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00094 *>         If ICOMPQ = 0, Q is not referenced.  Otherwise,
00095 *>         on entry, Q contains the eigenvectors of the partially solved
00096 *>         system which has been previously updated in matrix
00097 *>         multiplies with other partially solved eigensystems.
00098 *>         On exit, Q contains the trailing (N-K) updated eigenvectors
00099 *>         (those which were deflated) in its last N-K columns.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] LDQ
00103 *> \verbatim
00104 *>          LDQ is INTEGER
00105 *>         The leading dimension of the array Q.  LDQ >= max(1,N).
00106 *> \endverbatim
00107 *>
00108 *> \param[in] INDXQ
00109 *> \verbatim
00110 *>          INDXQ is INTEGER array, dimension (N)
00111 *>         The permutation which separately sorts the two sub-problems
00112 *>         in D into ascending order.  Note that elements in the second
00113 *>         half of this permutation must first have CUTPNT added to
00114 *>         their values in order to be accurate.
00115 *> \endverbatim
00116 *>
00117 *> \param[in,out] RHO
00118 *> \verbatim
00119 *>          RHO is DOUBLE PRECISION
00120 *>         On entry, the off-diagonal element associated with the rank-1
00121 *>         cut which originally split the two submatrices which are now
00122 *>         being recombined.
00123 *>         On exit, RHO has been modified to the value required by
00124 *>         DLAED3.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] CUTPNT
00128 *> \verbatim
00129 *>          CUTPNT is INTEGER
00130 *>         The location of the last eigenvalue in the leading
00131 *>         sub-matrix.  min(1,N) <= CUTPNT <= N.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] Z
00135 *> \verbatim
00136 *>          Z is DOUBLE PRECISION array, dimension (N)
00137 *>         On entry, Z contains the updating vector (the last row of
00138 *>         the first sub-eigenvector matrix and the first row of the
00139 *>         second sub-eigenvector matrix).
00140 *>         On exit, the contents of Z are destroyed by the updating
00141 *>         process.
00142 *> \endverbatim
00143 *>
00144 *> \param[out] DLAMDA
00145 *> \verbatim
00146 *>          DLAMDA is DOUBLE PRECISION array, dimension (N)
00147 *>         A copy of the first K eigenvalues which will be used by
00148 *>         DLAED3 to form the secular equation.
00149 *> \endverbatim
00150 *>
00151 *> \param[out] Q2
00152 *> \verbatim
00153 *>          Q2 is DOUBLE PRECISION array, dimension (LDQ2,N)
00154 *>         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
00155 *>         a copy of the first K eigenvectors which will be used by
00156 *>         DLAED7 in a matrix multiply (DGEMM) to update the new
00157 *>         eigenvectors.
00158 *> \endverbatim
00159 *>
00160 *> \param[in] LDQ2
00161 *> \verbatim
00162 *>          LDQ2 is INTEGER
00163 *>         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
00164 *> \endverbatim
00165 *>
00166 *> \param[out] W
00167 *> \verbatim
00168 *>          W is DOUBLE PRECISION array, dimension (N)
00169 *>         The first k values of the final deflation-altered z-vector and
00170 *>         will be passed to DLAED3.
00171 *> \endverbatim
00172 *>
00173 *> \param[out] PERM
00174 *> \verbatim
00175 *>          PERM is INTEGER array, dimension (N)
00176 *>         The permutations (from deflation and sorting) to be applied
00177 *>         to each eigenblock.
00178 *> \endverbatim
00179 *>
00180 *> \param[out] GIVPTR
00181 *> \verbatim
00182 *>          GIVPTR is INTEGER
00183 *>         The number of Givens rotations which took place in this
00184 *>         subproblem.
00185 *> \endverbatim
00186 *>
00187 *> \param[out] GIVCOL
00188 *> \verbatim
00189 *>          GIVCOL is INTEGER array, dimension (2, N)
00190 *>         Each pair of numbers indicates a pair of columns to take place
00191 *>         in a Givens rotation.
00192 *> \endverbatim
00193 *>
00194 *> \param[out] GIVNUM
00195 *> \verbatim
00196 *>          GIVNUM is DOUBLE PRECISION array, dimension (2, N)
00197 *>         Each number indicates the S value to be used in the
00198 *>         corresponding Givens rotation.
00199 *> \endverbatim
00200 *>
00201 *> \param[out] INDXP
00202 *> \verbatim
00203 *>          INDXP is INTEGER array, dimension (N)
00204 *>         The permutation used to place deflated values of D at the end
00205 *>         of the array.  INDXP(1:K) points to the nondeflated D-values
00206 *>         and INDXP(K+1:N) points to the deflated eigenvalues.
00207 *> \endverbatim
00208 *>
00209 *> \param[out] INDX
00210 *> \verbatim
00211 *>          INDX is INTEGER array, dimension (N)
00212 *>         The permutation used to sort the contents of D into ascending
00213 *>         order.
00214 *> \endverbatim
00215 *>
00216 *> \param[out] INFO
00217 *> \verbatim
00218 *>          INFO is INTEGER
00219 *>          = 0:  successful exit.
00220 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00221 *> \endverbatim
00222 *
00223 *  Authors:
00224 *  ========
00225 *
00226 *> \author Univ. of Tennessee 
00227 *> \author Univ. of California Berkeley 
00228 *> \author Univ. of Colorado Denver 
00229 *> \author NAG Ltd. 
00230 *
00231 *> \date November 2011
00232 *
00233 *> \ingroup auxOTHERcomputational
00234 *
00235 *> \par Contributors:
00236 *  ==================
00237 *>
00238 *> Jeff Rutter, Computer Science Division, University of California
00239 *> at Berkeley, USA
00240 *
00241 *  =====================================================================
00242       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
00243      $                   CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
00244      $                   GIVCOL, GIVNUM, INDXP, INDX, INFO )
00245 *
00246 *  -- LAPACK computational routine (version 3.4.0) --
00247 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00248 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00249 *     November 2011
00250 *
00251 *     .. Scalar Arguments ..
00252       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
00253      $                   QSIZ
00254       DOUBLE PRECISION   RHO
00255 *     ..
00256 *     .. Array Arguments ..
00257       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
00258      $                   INDXQ( * ), PERM( * )
00259       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
00260      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
00261 *     ..
00262 *
00263 *  =====================================================================
00264 *
00265 *     .. Parameters ..
00266       DOUBLE PRECISION   MONE, ZERO, ONE, TWO, EIGHT
00267       PARAMETER          ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
00268      $                   TWO = 2.0D0, EIGHT = 8.0D0 )
00269 *     ..
00270 *     .. Local Scalars ..
00271 *
00272       INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
00273       DOUBLE PRECISION   C, EPS, S, T, TAU, TOL
00274 *     ..
00275 *     .. External Functions ..
00276       INTEGER            IDAMAX
00277       DOUBLE PRECISION   DLAMCH, DLAPY2
00278       EXTERNAL           IDAMAX, DLAMCH, DLAPY2
00279 *     ..
00280 *     .. External Subroutines ..
00281       EXTERNAL           DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
00282 *     ..
00283 *     .. Intrinsic Functions ..
00284       INTRINSIC          ABS, MAX, MIN, SQRT
00285 *     ..
00286 *     .. Executable Statements ..
00287 *
00288 *     Test the input parameters.
00289 *
00290       INFO = 0
00291 *
00292       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
00293          INFO = -1
00294       ELSE IF( N.LT.0 ) THEN
00295          INFO = -3
00296       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
00297          INFO = -4
00298       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00299          INFO = -7
00300       ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
00301          INFO = -10
00302       ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
00303          INFO = -14
00304       END IF
00305       IF( INFO.NE.0 ) THEN
00306          CALL XERBLA( 'DLAED8', -INFO )
00307          RETURN
00308       END IF
00309 *
00310 *     Need to initialize GIVPTR to O here in case of quick exit
00311 *     to prevent an unspecified code behavior (usually sigfault) 
00312 *     when IWORK array on entry to *stedc is not zeroed 
00313 *     (or at least some IWORK entries which used in *laed7 for GIVPTR).
00314 *
00315       GIVPTR = 0
00316 *
00317 *     Quick return if possible
00318 *
00319       IF( N.EQ.0 )
00320      $   RETURN
00321 *
00322       N1 = CUTPNT
00323       N2 = N - N1
00324       N1P1 = N1 + 1
00325 *
00326       IF( RHO.LT.ZERO ) THEN
00327          CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
00328       END IF
00329 *
00330 *     Normalize z so that norm(z) = 1
00331 *
00332       T = ONE / SQRT( TWO )
00333       DO 10 J = 1, N
00334          INDX( J ) = J
00335    10 CONTINUE
00336       CALL DSCAL( N, T, Z, 1 )
00337       RHO = ABS( TWO*RHO )
00338 *
00339 *     Sort the eigenvalues into increasing order
00340 *
00341       DO 20 I = CUTPNT + 1, N
00342          INDXQ( I ) = INDXQ( I ) + CUTPNT
00343    20 CONTINUE
00344       DO 30 I = 1, N
00345          DLAMDA( I ) = D( INDXQ( I ) )
00346          W( I ) = Z( INDXQ( I ) )
00347    30 CONTINUE
00348       I = 1
00349       J = CUTPNT + 1
00350       CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
00351       DO 40 I = 1, N
00352          D( I ) = DLAMDA( INDX( I ) )
00353          Z( I ) = W( INDX( I ) )
00354    40 CONTINUE
00355 *
00356 *     Calculate the allowable deflation tolerence
00357 *
00358       IMAX = IDAMAX( N, Z, 1 )
00359       JMAX = IDAMAX( N, D, 1 )
00360       EPS = DLAMCH( 'Epsilon' )
00361       TOL = EIGHT*EPS*ABS( D( JMAX ) )
00362 *
00363 *     If the rank-1 modifier is small enough, no more needs to be done
00364 *     except to reorganize Q so that its columns correspond with the
00365 *     elements in D.
00366 *
00367       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
00368          K = 0
00369          IF( ICOMPQ.EQ.0 ) THEN
00370             DO 50 J = 1, N
00371                PERM( J ) = INDXQ( INDX( J ) )
00372    50       CONTINUE
00373          ELSE
00374             DO 60 J = 1, N
00375                PERM( J ) = INDXQ( INDX( J ) )
00376                CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
00377    60       CONTINUE
00378             CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
00379      $                   LDQ )
00380          END IF
00381          RETURN
00382       END IF
00383 *
00384 *     If there are multiple eigenvalues then the problem deflates.  Here
00385 *     the number of equal eigenvalues are found.  As each equal
00386 *     eigenvalue is found, an elementary reflector is computed to rotate
00387 *     the corresponding eigensubspace so that the corresponding
00388 *     components of Z are zero in this new basis.
00389 *
00390       K = 0
00391       K2 = N + 1
00392       DO 70 J = 1, N
00393          IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
00394 *
00395 *           Deflate due to small z component.
00396 *
00397             K2 = K2 - 1
00398             INDXP( K2 ) = J
00399             IF( J.EQ.N )
00400      $         GO TO 110
00401          ELSE
00402             JLAM = J
00403             GO TO 80
00404          END IF
00405    70 CONTINUE
00406    80 CONTINUE
00407       J = J + 1
00408       IF( J.GT.N )
00409      $   GO TO 100
00410       IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
00411 *
00412 *        Deflate due to small z component.
00413 *
00414          K2 = K2 - 1
00415          INDXP( K2 ) = J
00416       ELSE
00417 *
00418 *        Check if eigenvalues are close enough to allow deflation.
00419 *
00420          S = Z( JLAM )
00421          C = Z( J )
00422 *
00423 *        Find sqrt(a**2+b**2) without overflow or
00424 *        destructive underflow.
00425 *
00426          TAU = DLAPY2( C, S )
00427          T = D( J ) - D( JLAM )
00428          C = C / TAU
00429          S = -S / TAU
00430          IF( ABS( T*C*S ).LE.TOL ) THEN
00431 *
00432 *           Deflation is possible.
00433 *
00434             Z( J ) = TAU
00435             Z( JLAM ) = ZERO
00436 *
00437 *           Record the appropriate Givens rotation
00438 *
00439             GIVPTR = GIVPTR + 1
00440             GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
00441             GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
00442             GIVNUM( 1, GIVPTR ) = C
00443             GIVNUM( 2, GIVPTR ) = S
00444             IF( ICOMPQ.EQ.1 ) THEN
00445                CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
00446      $                    Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
00447             END IF
00448             T = D( JLAM )*C*C + D( J )*S*S
00449             D( J ) = D( JLAM )*S*S + D( J )*C*C
00450             D( JLAM ) = T
00451             K2 = K2 - 1
00452             I = 1
00453    90       CONTINUE
00454             IF( K2+I.LE.N ) THEN
00455                IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
00456                   INDXP( K2+I-1 ) = INDXP( K2+I )
00457                   INDXP( K2+I ) = JLAM
00458                   I = I + 1
00459                   GO TO 90
00460                ELSE
00461                   INDXP( K2+I-1 ) = JLAM
00462                END IF
00463             ELSE
00464                INDXP( K2+I-1 ) = JLAM
00465             END IF
00466             JLAM = J
00467          ELSE
00468             K = K + 1
00469             W( K ) = Z( JLAM )
00470             DLAMDA( K ) = D( JLAM )
00471             INDXP( K ) = JLAM
00472             JLAM = J
00473          END IF
00474       END IF
00475       GO TO 80
00476   100 CONTINUE
00477 *
00478 *     Record the last eigenvalue.
00479 *
00480       K = K + 1
00481       W( K ) = Z( JLAM )
00482       DLAMDA( K ) = D( JLAM )
00483       INDXP( K ) = JLAM
00484 *
00485   110 CONTINUE
00486 *
00487 *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
00488 *     and Q2 respectively.  The eigenvalues/vectors which were not
00489 *     deflated go into the first K slots of DLAMDA and Q2 respectively,
00490 *     while those which were deflated go into the last N - K slots.
00491 *
00492       IF( ICOMPQ.EQ.0 ) THEN
00493          DO 120 J = 1, N
00494             JP = INDXP( J )
00495             DLAMDA( J ) = D( JP )
00496             PERM( J ) = INDXQ( INDX( JP ) )
00497   120    CONTINUE
00498       ELSE
00499          DO 130 J = 1, N
00500             JP = INDXP( J )
00501             DLAMDA( J ) = D( JP )
00502             PERM( J ) = INDXQ( INDX( JP ) )
00503             CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
00504   130    CONTINUE
00505       END IF
00506 *
00507 *     The deflated eigenvalues and their corresponding vectors go back
00508 *     into the last N - K slots of D and Q respectively.
00509 *
00510       IF( K.LT.N ) THEN
00511          IF( ICOMPQ.EQ.0 ) THEN
00512             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
00513          ELSE
00514             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
00515             CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
00516      $                   Q( 1, K+1 ), LDQ )
00517          END IF
00518       END IF
00519 *
00520       RETURN
00521 *
00522 *     End of DLAED8
00523 *
00524       END
 All Files Functions