LAPACK 3.3.0

claed8.f

Go to the documentation of this file.
00001       SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
00002      $                   Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
00003      $                   GIVCOL, GIVNUM, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     June 2010
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
00012       REAL               RHO
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
00016      $                   INDXQ( * ), PERM( * )
00017       REAL               D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
00018      $                   Z( * )
00019       COMPLEX            Q( LDQ, * ), Q2( LDQ2, * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  CLAED8 merges the two sets of eigenvalues together into a single
00026 *  sorted set.  Then it tries to deflate the size of the problem.
00027 *  There are two ways in which deflation can occur:  when two or more
00028 *  eigenvalues are close together or if there is a tiny element in the
00029 *  Z vector.  For each such occurrence the order of the related secular
00030 *  equation problem is reduced by one.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  K      (output) INTEGER
00036 *         Contains the number of non-deflated eigenvalues.
00037 *         This is the order of the related secular equation.
00038 *
00039 *  N      (input) INTEGER
00040 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00041 *
00042 *  QSIZ   (input) INTEGER
00043 *         The dimension of the unitary matrix used to reduce
00044 *         the dense or band matrix to tridiagonal form.
00045 *         QSIZ >= N if ICOMPQ = 1.
00046 *
00047 *  Q      (input/output) COMPLEX array, dimension (LDQ,N)
00048 *         On entry, Q contains the eigenvectors of the partially solved
00049 *         system which has been previously updated in matrix
00050 *         multiplies with other partially solved eigensystems.
00051 *         On exit, Q contains the trailing (N-K) updated eigenvectors
00052 *         (those which were deflated) in its last N-K columns.
00053 *
00054 *  LDQ    (input) INTEGER
00055 *         The leading dimension of the array Q.  LDQ >= max( 1, N ).
00056 *
00057 *  D      (input/output) REAL array, dimension (N)
00058 *         On entry, D contains the eigenvalues of the two submatrices to
00059 *         be combined.  On exit, D contains the trailing (N-K) updated
00060 *         eigenvalues (those which were deflated) sorted into increasing
00061 *         order.
00062 *
00063 *  RHO    (input/output) REAL
00064 *         Contains the off diagonal element associated with the rank-1
00065 *         cut which originally split the two submatrices which are now
00066 *         being recombined. RHO is modified during the computation to
00067 *         the value required by SLAED3.
00068 *
00069 *  CUTPNT (input) INTEGER
00070 *         Contains the location of the last eigenvalue in the leading
00071 *         sub-matrix.  MIN(1,N) <= CUTPNT <= N.
00072 *
00073 *  Z      (input) REAL array, dimension (N)
00074 *         On input this vector contains the updating vector (the last
00075 *         row of the first sub-eigenvector matrix and the first row of
00076 *         the second sub-eigenvector matrix).  The contents of Z are
00077 *         destroyed during the updating process.
00078 *
00079 *  DLAMDA (output) REAL array, dimension (N)
00080 *         Contains a copy of the first K eigenvalues which will be used
00081 *         by SLAED3 to form the secular equation.
00082 *
00083 *  Q2     (output) COMPLEX array, dimension (LDQ2,N)
00084 *         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
00085 *         Contains a copy of the first K eigenvectors which will be used
00086 *         by SLAED7 in a matrix multiply (SGEMM) to update the new
00087 *         eigenvectors.
00088 *
00089 *  LDQ2   (input) INTEGER
00090 *         The leading dimension of the array Q2.  LDQ2 >= max( 1, N ).
00091 *
00092 *  W      (output) REAL array, dimension (N)
00093 *         This will hold the first k values of the final
00094 *         deflation-altered z-vector and will be passed to SLAED3.
00095 *
00096 *  INDXP  (workspace) INTEGER array, dimension (N)
00097 *         This will contain the permutation used to place deflated
00098 *         values of D at the end of the array. On output INDXP(1:K)
00099 *         points to the nondeflated D-values and INDXP(K+1:N)
00100 *         points to the deflated eigenvalues.
00101 *
00102 *  INDX   (workspace) INTEGER array, dimension (N)
00103 *         This will contain the permutation used to sort the contents of
00104 *         D into ascending order.
00105 *
00106 *  INDXQ  (input) INTEGER array, dimension (N)
00107 *         This contains the permutation which separately sorts the two
00108 *         sub-problems in D into ascending order.  Note that elements in
00109 *         the second half of this permutation must first have CUTPNT
00110 *         added to their values in order to be accurate.
00111 *
00112 *  PERM   (output) INTEGER array, dimension (N)
00113 *         Contains the permutations (from deflation and sorting) to be
00114 *         applied to each eigenblock.
00115 *
00116 *  GIVPTR (output) INTEGER
00117 *         Contains the number of Givens rotations which took place in
00118 *         this subproblem.
00119 *
00120 *  GIVCOL (output) INTEGER array, dimension (2, N)
00121 *         Each pair of numbers indicates a pair of columns to take place
00122 *         in a Givens rotation.
00123 *
00124 *  GIVNUM (output) REAL array, dimension (2, N)
00125 *         Each number indicates the S value to be used in the
00126 *         corresponding Givens rotation.
00127 *
00128 *  INFO   (output) INTEGER
00129 *          = 0:  successful exit.
00130 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00131 *
00132 *  =====================================================================
00133 *
00134 *     .. Parameters ..
00135       REAL               MONE, ZERO, ONE, TWO, EIGHT
00136       PARAMETER          ( MONE = -1.0E0, ZERO = 0.0E0, ONE = 1.0E0,
00137      $                   TWO = 2.0E0, EIGHT = 8.0E0 )
00138 *     ..
00139 *     .. Local Scalars ..
00140       INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
00141       REAL               C, EPS, S, T, TAU, TOL
00142 *     ..
00143 *     .. External Functions ..
00144       INTEGER            ISAMAX
00145       REAL               SLAMCH, SLAPY2
00146       EXTERNAL           ISAMAX, SLAMCH, SLAPY2
00147 *     ..
00148 *     .. External Subroutines ..
00149       EXTERNAL           CCOPY, CLACPY, CSROT, SCOPY, SLAMRG, SSCAL,
00150      $                   XERBLA
00151 *     ..
00152 *     .. Intrinsic Functions ..
00153       INTRINSIC          ABS, MAX, MIN, SQRT
00154 *     ..
00155 *     .. Executable Statements ..
00156 *
00157 *     Test the input parameters.
00158 *
00159       INFO = 0
00160 *
00161       IF( N.LT.0 ) THEN
00162          INFO = -2
00163       ELSE IF( QSIZ.LT.N ) THEN
00164          INFO = -3
00165       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00166          INFO = -5
00167       ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
00168          INFO = -8
00169       ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
00170          INFO = -12
00171       END IF
00172       IF( INFO.NE.0 ) THEN
00173          CALL XERBLA( 'CLAED8', -INFO )
00174          RETURN
00175       END IF
00176 *
00177 *     Need to initialize GIVPTR to O here in case of quick exit
00178 *     to prevent an unspecified code behavior (usually sigfault) 
00179 *     when IWORK array on entry to *stedc is not zeroed 
00180 *     (or at least some IWORK entries which used in *laed7 for GIVPTR).
00181 *
00182       GIVPTR = 0
00183 *
00184 *     Quick return if possible
00185 *
00186       IF( N.EQ.0 )
00187      $   RETURN
00188 *
00189       N1 = CUTPNT
00190       N2 = N - N1
00191       N1P1 = N1 + 1
00192 *
00193       IF( RHO.LT.ZERO ) THEN
00194          CALL SSCAL( N2, MONE, Z( N1P1 ), 1 )
00195       END IF
00196 *
00197 *     Normalize z so that norm(z) = 1
00198 *
00199       T = ONE / SQRT( TWO )
00200       DO 10 J = 1, N
00201          INDX( J ) = J
00202    10 CONTINUE
00203       CALL SSCAL( N, T, Z, 1 )
00204       RHO = ABS( TWO*RHO )
00205 *
00206 *     Sort the eigenvalues into increasing order
00207 *
00208       DO 20 I = CUTPNT + 1, N
00209          INDXQ( I ) = INDXQ( I ) + CUTPNT
00210    20 CONTINUE
00211       DO 30 I = 1, N
00212          DLAMDA( I ) = D( INDXQ( I ) )
00213          W( I ) = Z( INDXQ( I ) )
00214    30 CONTINUE
00215       I = 1
00216       J = CUTPNT + 1
00217       CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
00218       DO 40 I = 1, N
00219          D( I ) = DLAMDA( INDX( I ) )
00220          Z( I ) = W( INDX( I ) )
00221    40 CONTINUE
00222 *
00223 *     Calculate the allowable deflation tolerance
00224 *
00225       IMAX = ISAMAX( N, Z, 1 )
00226       JMAX = ISAMAX( N, D, 1 )
00227       EPS = SLAMCH( 'Epsilon' )
00228       TOL = EIGHT*EPS*ABS( D( JMAX ) )
00229 *
00230 *     If the rank-1 modifier is small enough, no more needs to be done
00231 *     -- except to reorganize Q so that its columns correspond with the
00232 *     elements in D.
00233 *
00234       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
00235          K = 0
00236          DO 50 J = 1, N
00237             PERM( J ) = INDXQ( INDX( J ) )
00238             CALL CCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
00239    50    CONTINUE
00240          CALL CLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), LDQ )
00241          RETURN
00242       END IF
00243 *
00244 *     If there are multiple eigenvalues then the problem deflates.  Here
00245 *     the number of equal eigenvalues are found.  As each equal
00246 *     eigenvalue is found, an elementary reflector is computed to rotate
00247 *     the corresponding eigensubspace so that the corresponding
00248 *     components of Z are zero in this new basis.
00249 *
00250       K = 0
00251       K2 = N + 1
00252       DO 60 J = 1, N
00253          IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
00254 *
00255 *           Deflate due to small z component.
00256 *
00257             K2 = K2 - 1
00258             INDXP( K2 ) = J
00259             IF( J.EQ.N )
00260      $         GO TO 100
00261          ELSE
00262             JLAM = J
00263             GO TO 70
00264          END IF
00265    60 CONTINUE
00266    70 CONTINUE
00267       J = J + 1
00268       IF( J.GT.N )
00269      $   GO TO 90
00270       IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
00271 *
00272 *        Deflate due to small z component.
00273 *
00274          K2 = K2 - 1
00275          INDXP( K2 ) = J
00276       ELSE
00277 *
00278 *        Check if eigenvalues are close enough to allow deflation.
00279 *
00280          S = Z( JLAM )
00281          C = Z( J )
00282 *
00283 *        Find sqrt(a**2+b**2) without overflow or
00284 *        destructive underflow.
00285 *
00286          TAU = SLAPY2( C, S )
00287          T = D( J ) - D( JLAM )
00288          C = C / TAU
00289          S = -S / TAU
00290          IF( ABS( T*C*S ).LE.TOL ) THEN
00291 *
00292 *           Deflation is possible.
00293 *
00294             Z( J ) = TAU
00295             Z( JLAM ) = ZERO
00296 *
00297 *           Record the appropriate Givens rotation
00298 *
00299             GIVPTR = GIVPTR + 1
00300             GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
00301             GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
00302             GIVNUM( 1, GIVPTR ) = C
00303             GIVNUM( 2, GIVPTR ) = S
00304             CALL CSROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
00305      $                  Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
00306             T = D( JLAM )*C*C + D( J )*S*S
00307             D( J ) = D( JLAM )*S*S + D( J )*C*C
00308             D( JLAM ) = T
00309             K2 = K2 - 1
00310             I = 1
00311    80       CONTINUE
00312             IF( K2+I.LE.N ) THEN
00313                IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
00314                   INDXP( K2+I-1 ) = INDXP( K2+I )
00315                   INDXP( K2+I ) = JLAM
00316                   I = I + 1
00317                   GO TO 80
00318                ELSE
00319                   INDXP( K2+I-1 ) = JLAM
00320                END IF
00321             ELSE
00322                INDXP( K2+I-1 ) = JLAM
00323             END IF
00324             JLAM = J
00325          ELSE
00326             K = K + 1
00327             W( K ) = Z( JLAM )
00328             DLAMDA( K ) = D( JLAM )
00329             INDXP( K ) = JLAM
00330             JLAM = J
00331          END IF
00332       END IF
00333       GO TO 70
00334    90 CONTINUE
00335 *
00336 *     Record the last eigenvalue.
00337 *
00338       K = K + 1
00339       W( K ) = Z( JLAM )
00340       DLAMDA( K ) = D( JLAM )
00341       INDXP( K ) = JLAM
00342 *
00343   100 CONTINUE
00344 *
00345 *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
00346 *     and Q2 respectively.  The eigenvalues/vectors which were not
00347 *     deflated go into the first K slots of DLAMDA and Q2 respectively,
00348 *     while those which were deflated go into the last N - K slots.
00349 *
00350       DO 110 J = 1, N
00351          JP = INDXP( J )
00352          DLAMDA( J ) = D( JP )
00353          PERM( J ) = INDXQ( INDX( JP ) )
00354          CALL CCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
00355   110 CONTINUE
00356 *
00357 *     The deflated eigenvalues and their corresponding vectors go back
00358 *     into the last N - K slots of D and Q respectively.
00359 *
00360       IF( K.LT.N ) THEN
00361          CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
00362          CALL CLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ),
00363      $                LDQ )
00364       END IF
00365 *
00366       RETURN
00367 *
00368 *     End of CLAED8
00369 *
00370       END
 All Files Functions