SRC/slaeda.f

   1 *> \brief \b SLAEDA used by sstedc. Computes the Z vector determining the rank-one modification of the diagonal 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 SLAEDA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaeda.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaeda.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaeda.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
  22 *                          GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
  29 *      $                   PRMPTR( * ), QPTR( * )
  30 *       REAL               GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLAEDA computes the Z vector corresponding to the merge step in the
  40 *> CURLVLth step of the merge process with TLVLS steps for the CURPBMth
  41 *> problem.
  42 *> \endverbatim
  43 *
  44 *  Arguments:
  45 *  ==========
  46 *
  47 *> \param[in] N
  48 *> \verbatim
  49 *>          N is INTEGER
  50 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
  51 *> \endverbatim
  52 *>
  53 *> \param[in] TLVLS
  54 *> \verbatim
  55 *>          TLVLS is INTEGER
  56 *>         The total number of merging levels in the overall divide and
  57 *>         conquer tree.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] CURLVL
  61 *> \verbatim
  62 *>          CURLVL is INTEGER
  63 *>         The current level in the overall merge routine,
  64 *>         0 <= curlvl <= tlvls.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] CURPBM
  68 *> \verbatim
  69 *>          CURPBM is INTEGER
  70 *>         The current problem in the current level in the overall
  71 *>         merge routine (counting from upper left to lower right).
  72 *> \endverbatim
  73 *>
  74 *> \param[in] PRMPTR
  75 *> \verbatim
  76 *>          PRMPTR is INTEGER array, dimension (N lg N)
  77 *>         Contains a list of pointers which indicate where in PERM a
  78 *>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
  79 *>         indicates the size of the permutation and incidentally the
  80 *>         size of the full, non-deflated problem.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] PERM
  84 *> \verbatim
  85 *>          PERM is INTEGER array, dimension (N lg N)
  86 *>         Contains the permutations (from deflation and sorting) to be
  87 *>         applied to each eigenblock.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] GIVPTR
  91 *> \verbatim
  92 *>          GIVPTR is INTEGER array, dimension (N lg N)
  93 *>         Contains a list of pointers which indicate where in GIVCOL a
  94 *>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
  95 *>         indicates the number of Givens rotations.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] GIVCOL
  99 *> \verbatim
 100 *>          GIVCOL is INTEGER array, dimension (2, N lg N)
 101 *>         Each pair of numbers indicates a pair of columns to take place
 102 *>         in a Givens rotation.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] GIVNUM
 106 *> \verbatim
 107 *>          GIVNUM is REAL array, dimension (2, N lg N)
 108 *>         Each number indicates the S value to be used in the
 109 *>         corresponding Givens rotation.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] Q
 113 *> \verbatim
 114 *>          Q is REAL array, dimension (N**2)
 115 *>         Contains the square eigenblocks from previous levels, the
 116 *>         starting positions for blocks are given by QPTR.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] QPTR
 120 *> \verbatim
 121 *>          QPTR is INTEGER array, dimension (N+2)
 122 *>         Contains a list of pointers which indicate where in Q an
 123 *>         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
 124 *>         the size of the block.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] Z
 128 *> \verbatim
 129 *>          Z is REAL array, dimension (N)
 130 *>         On output this vector contains the updating vector (the last
 131 *>         row of the first sub-eigenvector matrix and the first row of
 132 *>         the second sub-eigenvector matrix).
 133 *> \endverbatim
 134 *>
 135 *> \param[out] ZTEMP
 136 *> \verbatim
 137 *>          ZTEMP is REAL array, dimension (N)
 138 *> \endverbatim
 139 *>
 140 *> \param[out] INFO
 141 *> \verbatim
 142 *>          INFO is INTEGER
 143 *>          = 0:  successful exit.
 144 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 145 *> \endverbatim
 146 *
 147 *  Authors:
 148 *  ========
 149 *
 150 *> \author Univ. of Tennessee
 151 *> \author Univ. of California Berkeley
 152 *> \author Univ. of Colorado Denver
 153 *> \author NAG Ltd.
 154 *
 155 *> \date September 2012
 156 *
 157 *> \ingroup auxOTHERcomputational
 158 *
 159 *> \par Contributors:
 160 *  ==================
 161 *>
 162 *> Jeff Rutter, Computer Science Division, University of California
 163 *> at Berkeley, USA
 164 *
 165 *  =====================================================================
 166       SUBROUTINE SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
 167      $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
 168 *
 169 *  -- LAPACK computational routine (version 3.4.2) --
 170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 172 *     September 2012
 173 *
 174 *     .. Scalar Arguments ..
 175       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
 176 *     ..
 177 *     .. Array Arguments ..
 178       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
 179      $                   PRMPTR( * ), QPTR( * )
 180       REAL               GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       REAL               ZERO, HALF, ONE
 187       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0 )
 188 *     ..
 189 *     .. Local Scalars ..
 190       INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
 191      $                   PTR, ZPTR1
 192 *     ..
 193 *     .. External Subroutines ..
 194       EXTERNAL           SCOPY, SGEMV, SROT, XERBLA
 195 *     ..
 196 *     .. Intrinsic Functions ..
 197       INTRINSIC          INT, REAL, SQRT
 198 *     ..
 199 *     .. Executable Statements ..
 200 *
 201 *     Test the input parameters.
 202 *
 203       INFO = 0
 204 *
 205       IF( N.LT.0 ) THEN
 206          INFO = -1
 207       END IF
 208       IF( INFO.NE.0 ) THEN
 209          CALL XERBLA( 'SLAEDA', -INFO )
 210          RETURN
 211       END IF
 212 *
 213 *     Quick return if possible
 214 *
 215       IF( N.EQ.0 )
 216      $   RETURN
 217 *
 218 *     Determine location of first number in second half.
 219 *
 220       MID = N / 2 + 1
 221 *
 222 *     Gather last/first rows of appropriate eigenblocks into center of Z
 223 *
 224       PTR = 1
 225 *
 226 *     Determine location of lowest level subproblem in the full storage
 227 *     scheme
 228 *
 229       CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
 230 *
 231 *     Determine size of these matrices.  We add HALF to the value of
 232 *     the SQRT in case the machine underestimates one of these square
 233 *     roots.
 234 *
 235       BSIZ1 = INT( HALF+SQRT( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
 236       BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
 237       DO 10 K = 1, MID - BSIZ1 - 1
 238          Z( K ) = ZERO
 239    10 CONTINUE
 240       CALL SCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
 241      $            Z( MID-BSIZ1 ), 1 )
 242       CALL SCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
 243       DO 20 K = MID + BSIZ2, N
 244          Z( K ) = ZERO
 245    20 CONTINUE
 246 *
 247 *     Loop through remaining levels 1 -> CURLVL applying the Givens
 248 *     rotations and permutation and then multiplying the center matrices
 249 *     against the current Z.
 250 *
 251       PTR = 2**TLVLS + 1
 252       DO 70 K = 1, CURLVL - 1
 253          CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
 254          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
 255          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
 256          ZPTR1 = MID - PSIZ1
 257 *
 258 *       Apply Givens at CURR and CURR+1
 259 *
 260          DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
 261             CALL SROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
 262      $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
 263      $                 GIVNUM( 2, I ) )
 264    30    CONTINUE
 265          DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
 266             CALL SROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
 267      $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
 268      $                 GIVNUM( 2, I ) )
 269    40    CONTINUE
 270          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
 271          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
 272          DO 50 I = 0, PSIZ1 - 1
 273             ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
 274    50    CONTINUE
 275          DO 60 I = 0, PSIZ2 - 1
 276             ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
 277    60    CONTINUE
 278 *
 279 *        Multiply Blocks at CURR and CURR+1
 280 *
 281 *        Determine size of these matrices.  We add HALF to the value of
 282 *        the SQRT in case the machine underestimates one of these
 283 *        square roots.
 284 *
 285          BSIZ1 = INT( HALF+SQRT( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
 286          BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+
 287      $           1 ) ) ) )
 288          IF( BSIZ1.GT.0 ) THEN
 289             CALL SGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
 290      $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
 291          END IF
 292          CALL SCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
 293      $               1 )
 294          IF( BSIZ2.GT.0 ) THEN
 295             CALL SGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
 296      $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
 297          END IF
 298          CALL SCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
 299      $               Z( MID+BSIZ2 ), 1 )
 300 *
 301          PTR = PTR + 2**( TLVLS-K )
 302    70 CONTINUE
 303 *
 304       RETURN
 305 *
 306 *     End of SLAEDA
 307 *
 308       END