SRC/claed7.f

   1 *> \brief \b CLAED7 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CLAED7 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claed7.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claed7.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claed7.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
  22 *                          LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
  23 *                          GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
  24 *                          INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       INTEGER            CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
  28 *      $                   TLVLS
  29 *       REAL               RHO
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
  33 *      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
  34 *       REAL               D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
  35 *       COMPLEX            Q( LDQ, * ), WORK( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> CLAED7 computes the updated eigensystem of a diagonal
  45 *> matrix after modification by a rank-one symmetric matrix. This
  46 *> routine is used only for the eigenproblem which requires all
  47 *> eigenvalues and optionally eigenvectors of a dense or banded
  48 *> Hermitian matrix that has been reduced to tridiagonal form.
  49 *>
  50 *>   T = Q(in) ( D(in) + RHO * Z*Z**H ) Q**H(in) = Q(out) * D(out) * Q**H(out)
  51 *>
  52 *>   where Z = Q**Hu, u is a vector of length N with ones in the
  53 *>   CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
  54 *>
  55 *>    The eigenvectors of the original matrix are stored in Q, and the
  56 *>    eigenvalues are in D.  The algorithm consists of three stages:
  57 *>
  58 *>       The first stage consists of deflating the size of the problem
  59 *>       when there are multiple eigenvalues or if there is a zero in
  60 *>       the Z vector.  For each such occurrence the dimension of the
  61 *>       secular equation problem is reduced by one.  This stage is
  62 *>       performed by the routine SLAED2.
  63 *>
  64 *>       The second stage consists of calculating the updated
  65 *>       eigenvalues. This is done by finding the roots of the secular
  66 *>       equation via the routine SLAED4 (as called by SLAED3).
  67 *>       This routine also calculates the eigenvectors of the current
  68 *>       problem.
  69 *>
  70 *>       The final stage consists of computing the updated eigenvectors
  71 *>       directly using the updated eigenvalues.  The eigenvectors for
  72 *>       the current problem are multiplied with the eigenvectors from
  73 *>       the overall problem.
  74 *> \endverbatim
  75 *
  76 *  Arguments:
  77 *  ==========
  78 *
  79 *> \param[in] N
  80 *> \verbatim
  81 *>          N is INTEGER
  82 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] CUTPNT
  86 *> \verbatim
  87 *>          CUTPNT is INTEGER
  88 *>         Contains the location of the last eigenvalue in the leading
  89 *>         sub-matrix.  min(1,N) <= CUTPNT <= N.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] QSIZ
  93 *> \verbatim
  94 *>          QSIZ is INTEGER
  95 *>         The dimension of the unitary matrix used to reduce
  96 *>         the full matrix to tridiagonal form.  QSIZ >= N.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] TLVLS
 100 *> \verbatim
 101 *>          TLVLS is INTEGER
 102 *>         The total number of merging levels in the overall divide and
 103 *>         conquer tree.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] CURLVL
 107 *> \verbatim
 108 *>          CURLVL is INTEGER
 109 *>         The current level in the overall merge routine,
 110 *>         0 <= curlvl <= tlvls.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] CURPBM
 114 *> \verbatim
 115 *>          CURPBM is INTEGER
 116 *>         The current problem in the current level in the overall
 117 *>         merge routine (counting from upper left to lower right).
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] D
 121 *> \verbatim
 122 *>          D is REAL array, dimension (N)
 123 *>         On entry, the eigenvalues of the rank-1-perturbed matrix.
 124 *>         On exit, the eigenvalues of the repaired matrix.
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] Q
 128 *> \verbatim
 129 *>          Q is COMPLEX array, dimension (LDQ,N)
 130 *>         On entry, the eigenvectors of the rank-1-perturbed matrix.
 131 *>         On exit, the eigenvectors of the repaired tridiagonal matrix.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LDQ
 135 *> \verbatim
 136 *>          LDQ is INTEGER
 137 *>         The leading dimension of the array Q.  LDQ >= max(1,N).
 138 *> \endverbatim
 139 *>
 140 *> \param[in] RHO
 141 *> \verbatim
 142 *>          RHO is REAL
 143 *>         Contains the subdiagonal element used to create the rank-1
 144 *>         modification.
 145 *> \endverbatim
 146 *>
 147 *> \param[out] INDXQ
 148 *> \verbatim
 149 *>          INDXQ is INTEGER array, dimension (N)
 150 *>         This contains the permutation which will reintegrate the
 151 *>         subproblem just solved back into sorted order,
 152 *>         ie. D( INDXQ( I = 1, N ) ) will be in ascending order.
 153 *> \endverbatim
 154 *>
 155 *> \param[out] IWORK
 156 *> \verbatim
 157 *>          IWORK is INTEGER array, dimension (4*N)
 158 *> \endverbatim
 159 *>
 160 *> \param[out] RWORK
 161 *> \verbatim
 162 *>          RWORK is REAL array,
 163 *>                                 dimension (3*N+2*QSIZ*N)
 164 *> \endverbatim
 165 *>
 166 *> \param[out] WORK
 167 *> \verbatim
 168 *>          WORK is COMPLEX array, dimension (QSIZ*N)
 169 *> \endverbatim
 170 *>
 171 *> \param[in,out] QSTORE
 172 *> \verbatim
 173 *>          QSTORE is REAL array, dimension (N**2+1)
 174 *>         Stores eigenvectors of submatrices encountered during
 175 *>         divide and conquer, packed together. QPTR points to
 176 *>         beginning of the submatrices.
 177 *> \endverbatim
 178 *>
 179 *> \param[in,out] QPTR
 180 *> \verbatim
 181 *>          QPTR is INTEGER array, dimension (N+2)
 182 *>         List of indices pointing to beginning of submatrices stored
 183 *>         in QSTORE. The submatrices are numbered starting at the
 184 *>         bottom left of the divide and conquer tree, from left to
 185 *>         right and bottom to top.
 186 *> \endverbatim
 187 *>
 188 *> \param[in] PRMPTR
 189 *> \verbatim
 190 *>          PRMPTR is INTEGER array, dimension (N lg N)
 191 *>         Contains a list of pointers which indicate where in PERM a
 192 *>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
 193 *>         indicates the size of the permutation and also the size of
 194 *>         the full, non-deflated problem.
 195 *> \endverbatim
 196 *>
 197 *> \param[in] PERM
 198 *> \verbatim
 199 *>          PERM is INTEGER array, dimension (N lg N)
 200 *>         Contains the permutations (from deflation and sorting) to be
 201 *>         applied to each eigenblock.
 202 *> \endverbatim
 203 *>
 204 *> \param[in] GIVPTR
 205 *> \verbatim
 206 *>          GIVPTR is INTEGER array, dimension (N lg N)
 207 *>         Contains a list of pointers which indicate where in GIVCOL a
 208 *>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
 209 *>         indicates the number of Givens rotations.
 210 *> \endverbatim
 211 *>
 212 *> \param[in] GIVCOL
 213 *> \verbatim
 214 *>          GIVCOL is INTEGER array, dimension (2, N lg N)
 215 *>         Each pair of numbers indicates a pair of columns to take place
 216 *>         in a Givens rotation.
 217 *> \endverbatim
 218 *>
 219 *> \param[in] GIVNUM
 220 *> \verbatim
 221 *>          GIVNUM is REAL array, dimension (2, N lg N)
 222 *>         Each number indicates the S value to be used in the
 223 *>         corresponding Givens rotation.
 224 *> \endverbatim
 225 *>
 226 *> \param[out] INFO
 227 *> \verbatim
 228 *>          INFO is INTEGER
 229 *>          = 0:  successful exit.
 230 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 231 *>          > 0:  if INFO = 1, an eigenvalue did not converge
 232 *> \endverbatim
 233 *
 234 *  Authors:
 235 *  ========
 236 *
 237 *> \author Univ. of Tennessee
 238 *> \author Univ. of California Berkeley
 239 *> \author Univ. of Colorado Denver
 240 *> \author NAG Ltd.
 241 *
 242 *> \date June 2016
 243 *
 244 *> \ingroup complexOTHERcomputational
 245 *
 246 *  =====================================================================
 247       SUBROUTINE CLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
 248      $                   LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
 249      $                   GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
 250      $                   INFO )
 251 *
 252 *  -- LAPACK computational routine (version 3.6.1) --
 253 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 254 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 255 *     June 2016
 256 *
 257 *     .. Scalar Arguments ..
 258       INTEGER            CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
 259      $                   TLVLS
 260       REAL               RHO
 261 *     ..
 262 *     .. Array Arguments ..
 263       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
 264      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
 265       REAL               D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
 266       COMPLEX            Q( LDQ, * ), WORK( * )
 267 *     ..
 268 *
 269 *  =====================================================================
 270 *
 271 *     .. Local Scalars ..
 272       INTEGER            COLTYP, CURR, I, IDLMDA, INDX,
 273      $                   INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
 274 *     ..
 275 *     .. External Subroutines ..
 276       EXTERNAL           CLACRM, CLAED8, SLAED9, SLAEDA, SLAMRG, XERBLA
 277 *     ..
 278 *     .. Intrinsic Functions ..
 279       INTRINSIC          MAX, MIN
 280 *     ..
 281 *     .. Executable Statements ..
 282 *
 283 *     Test the input parameters.
 284 *
 285       INFO = 0
 286 *
 287 *     IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
 288 *        INFO = -1
 289 *     ELSE IF( N.LT.0 ) THEN
 290       IF( N.LT.0 ) THEN
 291          INFO = -1
 292       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
 293          INFO = -2
 294       ELSE IF( QSIZ.LT.N ) THEN
 295          INFO = -3
 296       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
 297          INFO = -9
 298       END IF
 299       IF( INFO.NE.0 ) THEN
 300          CALL XERBLA( 'CLAED7', -INFO )
 301          RETURN
 302       END IF
 303 *
 304 *     Quick return if possible
 305 *
 306       IF( N.EQ.0 )
 307      $   RETURN
 308 *
 309 *     The following values are for bookkeeping purposes only.  They are
 310 *     integer pointers which indicate the portion of the workspace
 311 *     used by a particular array in SLAED2 and SLAED3.
 312 *
 313       IZ = 1
 314       IDLMDA = IZ + N
 315       IW = IDLMDA + N
 316       IQ = IW + N
 317 *
 318       INDX = 1
 319       INDXC = INDX + N
 320       COLTYP = INDXC + N
 321       INDXP = COLTYP + N
 322 *
 323 *     Form the z-vector which consists of the last row of Q_1 and the
 324 *     first row of Q_2.
 325 *
 326       PTR = 1 + 2**TLVLS
 327       DO 10 I = 1, CURLVL - 1
 328          PTR = PTR + 2**( TLVLS-I )
 329    10 CONTINUE
 330       CURR = PTR + CURPBM
 331       CALL SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
 332      $             GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
 333      $             RWORK( IZ+N ), INFO )
 334 *
 335 *     When solving the final problem, we no longer need the stored data,
 336 *     so we will overwrite the data from this level onto the previously
 337 *     used storage space.
 338 *
 339       IF( CURLVL.EQ.TLVLS ) THEN
 340          QPTR( CURR ) = 1
 341          PRMPTR( CURR ) = 1
 342          GIVPTR( CURR ) = 1
 343       END IF
 344 *
 345 *     Sort and Deflate eigenvalues.
 346 *
 347       CALL CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
 348      $             RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
 349      $             IWORK( INDXP ), IWORK( INDX ), INDXQ,
 350      $             PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
 351      $             GIVCOL( 1, GIVPTR( CURR ) ),
 352      $             GIVNUM( 1, GIVPTR( CURR ) ), INFO )
 353       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
 354       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
 355 *
 356 *     Solve Secular Equation.
 357 *
 358       IF( K.NE.0 ) THEN
 359          CALL SLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
 360      $                RWORK( IDLMDA ), RWORK( IW ),
 361      $                QSTORE( QPTR( CURR ) ), K, INFO )
 362          CALL CLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
 363      $                LDQ, RWORK( IQ ) )
 364          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
 365          IF( INFO.NE.0 ) THEN
 366             RETURN
 367          END IF
 368 *
 369 *     Prepare the INDXQ sorting premutation.
 370 *
 371          N1 = K
 372          N2 = N - K
 373          CALL SLAMRG( N1, N2, D, 1, -1, INDXQ )
 374       ELSE
 375          QPTR( CURR+1 ) = QPTR( CURR )
 376          DO 20 I = 1, N
 377             INDXQ( I ) = I
 378    20    CONTINUE
 379       END IF
 380 *
 381       RETURN
 382 *
 383 *     End of CLAED7
 384 *
 385       END