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