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