SRC/slaed7.f

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