1*> \brief \b DLAEDA 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 DLAEDA + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaeda.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaeda.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaeda.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAEDA( 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* DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DLAEDA 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 December 2016 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 DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, 167 $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO ) 168* 169* -- LAPACK computational routine (version 3.7.0) -- 170* -- LAPACK is a software package provided by Univ. of Tennessee, -- 171* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 172* December 2016 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 DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 181* .. 182* 183* ===================================================================== 184* 185* .. Parameters .. 186 DOUBLE PRECISION ZERO, HALF, ONE 187 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 ) 188* .. 189* .. Local Scalars .. 190 INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2, 191 $ PTR, ZPTR1 192* .. 193* .. External Subroutines .. 194 EXTERNAL DCOPY, DGEMV, DROT, XERBLA 195* .. 196* .. Intrinsic Functions .. 197 INTRINSIC DBLE, INT, 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( 'DLAEDA', -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( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 236 BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) ) 237 DO 10 K = 1, MID - BSIZ1 - 1 238 Z( K ) = ZERO 239 10 CONTINUE 240 CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1, 241 $ Z( MID-BSIZ1 ), 1 ) 242 CALL DCOPY( 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 DROT( 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 DROT( 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( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 286 BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+ 287 $ 1 ) ) ) ) 288 IF( BSIZ1.GT.0 ) THEN 289 CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ), 290 $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 ) 291 END IF 292 CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ), 293 $ 1 ) 294 IF( BSIZ2.GT.0 ) THEN 295 CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ), 296 $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 ) 297 END IF 298 CALL DCOPY( 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 DLAEDA 307* 308 END 309