1 SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, 2 $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO ) 3* 4* -- LAPACK routine (instrumented to count operations, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* September 30, 1994 8* 9* .. Scalar Arguments .. 10 INTEGER CURLVL, CURPBM, INFO, N, TLVLS 11* .. 12* .. Array Arguments .. 13 INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ), 14 $ PRMPTR( * ), QPTR( * ) 15 DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 16* .. 17* Common block to return operation count and iteration count 18* ITCNT is unchanged, OPS is only incremented 19* .. Common blocks .. 20 COMMON / LATIME / OPS, ITCNT 21* .. 22* .. Scalars in Common .. 23 DOUBLE PRECISION ITCNT, OPS 24* .. 25* 26* Purpose 27* ======= 28* 29* DLAEDA computes the Z vector corresponding to the merge step in the 30* CURLVLth step of the merge process with TLVLS steps for the CURPBMth 31* problem. 32* 33* Arguments 34* ========= 35* 36* N (input) INTEGER 37* The dimension of the symmetric tridiagonal matrix. N >= 0. 38* 39* TLVLS (input) INTEGER 40* The total number of merging levels in the overall divide and 41* conquer tree. 42* 43* CURLVL (input) INTEGER 44* The current level in the overall merge routine, 45* 0 <= curlvl <= tlvls. 46* 47* CURPBM (input) INTEGER 48* The current problem in the current level in the overall 49* merge routine (counting from upper left to lower right). 50* 51* PRMPTR (input) INTEGER array, dimension (N lg N) 52* Contains a list of pointers which indicate where in PERM a 53* level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) 54* indicates the size of the permutation and incidentally the 55* size of the full, non-deflated problem. 56* 57* PERM (input) INTEGER array, dimension (N lg N) 58* Contains the permutations (from deflation and sorting) to be 59* applied to each eigenblock. 60* 61* GIVPTR (input) INTEGER array, dimension (N lg N) 62* Contains a list of pointers which indicate where in GIVCOL a 63* level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) 64* indicates the number of Givens rotations. 65* 66* GIVCOL (input) INTEGER array, dimension (2, N lg N) 67* Each pair of numbers indicates a pair of columns to take place 68* in a Givens rotation. 69* 70* GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N) 71* Each number indicates the S value to be used in the 72* corresponding Givens rotation. 73* 74* Q (input) DOUBLE PRECISION array, dimension (N**2) 75* Contains the square eigenblocks from previous levels, the 76* starting positions for blocks are given by QPTR. 77* 78* QPTR (input) INTEGER array, dimension (N+2) 79* Contains a list of pointers which indicate where in Q an 80* eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates 81* the size of the block. 82* 83* Z (output) DOUBLE PRECISION array, dimension (N) 84* On output this vector contains the updating vector (the last 85* row of the first sub-eigenvector matrix and the first row of 86* the second sub-eigenvector matrix). 87* 88* ZTEMP (workspace) DOUBLE PRECISION array, dimension (N) 89* 90* INFO (output) INTEGER 91* = 0: successful exit. 92* < 0: if INFO = -i, the i-th argument had an illegal value. 93* 94* Further Details 95* =============== 96* 97* Based on contributions by 98* Jeff Rutter, Computer Science Division, University of California 99* at Berkeley, USA 100* 101* ===================================================================== 102* 103* .. Parameters .. 104 DOUBLE PRECISION ZERO, HALF, ONE 105 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 ) 106* .. 107* .. Local Scalars .. 108 INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2, 109 $ PTR, ZPTR1 110* .. 111* .. External Subroutines .. 112 EXTERNAL DCOPY, DGEMV, DROT, XERBLA 113* .. 114* .. Intrinsic Functions .. 115 INTRINSIC DBLE, INT, SQRT 116* .. 117* .. Executable Statements .. 118* 119* Test the input parameters. 120* 121 INFO = 0 122* 123 IF( N.LT.0 ) THEN 124 INFO = -1 125 END IF 126 IF( INFO.NE.0 ) THEN 127 CALL XERBLA( 'DLAEDA', -INFO ) 128 RETURN 129 END IF 130* 131* Quick return if possible 132* 133 IF( N.EQ.0 ) 134 $ RETURN 135* 136* Determine location of first number in second half. 137* 138 MID = N / 2 + 1 139* 140* Gather last/first rows of appropriate eigenblocks into center of Z 141* 142 PTR = 1 143* 144* Determine location of lowest level subproblem in the full storage 145* scheme 146* 147 CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1 148* 149* Determine size of these matrices. We add HALF to the value of 150* the SQRT in case the machine underestimates one of these square 151* roots. 152* 153 OPS = OPS + 8 154 BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 155 BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) ) 156 DO 10 K = 1, MID - BSIZ1 - 1 157 Z( K ) = ZERO 158 10 CONTINUE 159 CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1, 160 $ Z( MID-BSIZ1 ), 1 ) 161 CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 ) 162 DO 20 K = MID + BSIZ2, N 163 Z( K ) = ZERO 164 20 CONTINUE 165* 166* Loop thru remaining levels 1 -> CURLVL applying the Givens 167* rotations and permutation and then multiplying the center matrices 168* against the current Z. 169* 170 PTR = 2**TLVLS + 1 171 DO 70 K = 1, CURLVL - 1 172 CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1 173 PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR ) 174 PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 ) 175 ZPTR1 = MID - PSIZ1 176* 177* Apply Givens at CURR and CURR+1 178* 179 OPS = OPS + 6*( GIVPTR( CURR+2 )-GIVPTR( CURR ) ) 180 DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1 181 CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1, 182 $ Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ), 183 $ GIVNUM( 2, I ) ) 184 30 CONTINUE 185 DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1 186 CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1, 187 $ Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ), 188 $ GIVNUM( 2, I ) ) 189 40 CONTINUE 190 PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR ) 191 PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 ) 192 DO 50 I = 0, PSIZ1 - 1 193 ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 ) 194 50 CONTINUE 195 DO 60 I = 0, PSIZ2 - 1 196 ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 ) 197 60 CONTINUE 198* 199* Multiply Blocks at CURR and CURR+1 200* 201* Determine size of these matrices. We add HALF to the value of 202* the SQRT in case the machine underestimates one of these 203* square roots. 204* 205 OPS = OPS + 8 206 BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 207 BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+ 208 $ 1 ) ) ) ) 209 IF( BSIZ1.GT.0 ) THEN 210 OPS = OPS + 2*BSIZ1*BSIZ1 211 CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ), 212 $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 ) 213 END IF 214 CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ), 215 $ 1 ) 216 IF( BSIZ2.GT.0 ) THEN 217 OPS = OPS + 2*BSIZ2*BSIZ2 218 CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ), 219 $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 ) 220 END IF 221 CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1, 222 $ Z( MID+BSIZ2 ), 1 ) 223* 224 PTR = PTR + 2**( TLVLS-K ) 225 70 CONTINUE 226* 227 RETURN 228* 229* End of DLAEDA 230* 231 END 232