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