001:       SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
002:      $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
013:      $                   PRMPTR( * ), QPTR( * )
014:       DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLAEDA computes the Z vector corresponding to the merge step in the
021: *  CURLVLth step of the merge process with TLVLS steps for the CURPBMth
022: *  problem.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  N      (input) INTEGER
028: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
029: *
030: *  TLVLS  (input) INTEGER
031: *         The total number of merging levels in the overall divide and
032: *         conquer tree.
033: *
034: *  CURLVL (input) INTEGER
035: *         The current level in the overall merge routine,
036: *         0 <= curlvl <= tlvls.
037: *
038: *  CURPBM (input) INTEGER
039: *         The current problem in the current level in the overall
040: *         merge routine (counting from upper left to lower right).
041: *
042: *  PRMPTR (input) INTEGER array, dimension (N lg N)
043: *         Contains a list of pointers which indicate where in PERM a
044: *         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
045: *         indicates the size of the permutation and incidentally the
046: *         size of the full, non-deflated problem.
047: *
048: *  PERM   (input) INTEGER array, dimension (N lg N)
049: *         Contains the permutations (from deflation and sorting) to be
050: *         applied to each eigenblock.
051: *
052: *  GIVPTR (input) INTEGER array, dimension (N lg N)
053: *         Contains a list of pointers which indicate where in GIVCOL a
054: *         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
055: *         indicates the number of Givens rotations.
056: *
057: *  GIVCOL (input) INTEGER array, dimension (2, N lg N)
058: *         Each pair of numbers indicates a pair of columns to take place
059: *         in a Givens rotation.
060: *
061: *  GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
062: *         Each number indicates the S value to be used in the
063: *         corresponding Givens rotation.
064: *
065: *  Q      (input) DOUBLE PRECISION array, dimension (N**2)
066: *         Contains the square eigenblocks from previous levels, the
067: *         starting positions for blocks are given by QPTR.
068: *
069: *  QPTR   (input) INTEGER array, dimension (N+2)
070: *         Contains a list of pointers which indicate where in Q an
071: *         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
072: *         the size of the block.
073: *
074: *  Z      (output) DOUBLE PRECISION array, dimension (N)
075: *         On output this vector contains the updating vector (the last
076: *         row of the first sub-eigenvector matrix and the first row of
077: *         the second sub-eigenvector matrix).
078: *
079: *  ZTEMP  (workspace) DOUBLE PRECISION array, dimension (N)
080: *
081: *  INFO   (output) INTEGER
082: *          = 0:  successful exit.
083: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
084: *
085: *  Further Details
086: *  ===============
087: *
088: *  Based on contributions by
089: *     Jeff Rutter, Computer Science Division, University of California
090: *     at Berkeley, USA
091: *
092: *  =====================================================================
093: *
094: *     .. Parameters ..
095:       DOUBLE PRECISION   ZERO, HALF, ONE
096:       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
097: *     ..
098: *     .. Local Scalars ..
099:       INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
100:      $                   PTR, ZPTR1
101: *     ..
102: *     .. External Subroutines ..
103:       EXTERNAL           DCOPY, DGEMV, DROT, XERBLA
104: *     ..
105: *     .. Intrinsic Functions ..
106:       INTRINSIC          DBLE, INT, SQRT
107: *     ..
108: *     .. Executable Statements ..
109: *
110: *     Test the input parameters.
111: *
112:       INFO = 0
113: *
114:       IF( N.LT.0 ) THEN
115:          INFO = -1
116:       END IF
117:       IF( INFO.NE.0 ) THEN
118:          CALL XERBLA( 'DLAEDA', -INFO )
119:          RETURN
120:       END IF
121: *
122: *     Quick return if possible
123: *
124:       IF( N.EQ.0 )
125:      $   RETURN
126: *
127: *     Determine location of first number in second half.
128: *
129:       MID = N / 2 + 1
130: *
131: *     Gather last/first rows of appropriate eigenblocks into center of Z
132: *
133:       PTR = 1
134: *
135: *     Determine location of lowest level subproblem in the full storage
136: *     scheme
137: *
138:       CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
139: *
140: *     Determine size of these matrices.  We add HALF to the value of
141: *     the SQRT in case the machine underestimates one of these square
142: *     roots.
143: *
144:       BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
145:       BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
146:       DO 10 K = 1, MID - BSIZ1 - 1
147:          Z( K ) = ZERO
148:    10 CONTINUE
149:       CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
150:      $            Z( MID-BSIZ1 ), 1 )
151:       CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
152:       DO 20 K = MID + BSIZ2, N
153:          Z( K ) = ZERO
154:    20 CONTINUE
155: *
156: *     Loop thru remaining levels 1 -> CURLVL applying the Givens
157: *     rotations and permutation and then multiplying the center matrices
158: *     against the current Z.
159: *
160:       PTR = 2**TLVLS + 1
161:       DO 70 K = 1, CURLVL - 1
162:          CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
163:          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
164:          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
165:          ZPTR1 = MID - PSIZ1
166: *
167: *       Apply Givens at CURR and CURR+1
168: *
169:          DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
170:             CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
171:      $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
172:      $                 GIVNUM( 2, I ) )
173:    30    CONTINUE
174:          DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
175:             CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
176:      $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
177:      $                 GIVNUM( 2, I ) )
178:    40    CONTINUE
179:          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
180:          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
181:          DO 50 I = 0, PSIZ1 - 1
182:             ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
183:    50    CONTINUE
184:          DO 60 I = 0, PSIZ2 - 1
185:             ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
186:    60    CONTINUE
187: *
188: *        Multiply Blocks at CURR and CURR+1
189: *
190: *        Determine size of these matrices.  We add HALF to the value of
191: *        the SQRT in case the machine underestimates one of these
192: *        square roots.
193: *
194:          BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
195:          BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
196:      $           1 ) ) ) )
197:          IF( BSIZ1.GT.0 ) THEN
198:             CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
199:      $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
200:          END IF
201:          CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
202:      $               1 )
203:          IF( BSIZ2.GT.0 ) THEN
204:             CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
205:      $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
206:          END IF
207:          CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
208:      $               Z( MID+BSIZ2 ), 1 )
209: *
210:          PTR = PTR + 2**( TLVLS-K )
211:    70 CONTINUE
212: *
213:       RETURN
214: *
215: *     End of DLAEDA
216: *
217:       END
218: