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*> \ingroup auxOTHERcomputational
156*
157*> \par Contributors:
158*  ==================
159*>
160*> Jeff Rutter, Computer Science Division, University of California
161*> at Berkeley, USA
162*
163*  =====================================================================
164      SUBROUTINE SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
165     $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
166*
167*  -- LAPACK computational routine --
168*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
169*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*
171*     .. Scalar Arguments ..
172      INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
173*     ..
174*     .. Array Arguments ..
175      INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
176     $                   PRMPTR( * ), QPTR( * )
177      REAL               GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
178*     ..
179*
180*  =====================================================================
181*
182*     .. Parameters ..
183      REAL               ZERO, HALF, ONE
184      PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0 )
185*     ..
186*     .. Local Scalars ..
187      INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
188     $                   PTR, ZPTR1
189*     ..
190*     .. External Subroutines ..
191      EXTERNAL           SCOPY, SGEMV, SROT, XERBLA
192*     ..
193*     .. Intrinsic Functions ..
194      INTRINSIC          INT, REAL, SQRT
195*     ..
196*     .. Executable Statements ..
197*
198*     Test the input parameters.
199*
200      INFO = 0
201*
202      IF( N.LT.0 ) THEN
203         INFO = -1
204      END IF
205      IF( INFO.NE.0 ) THEN
206         CALL XERBLA( 'SLAEDA', -INFO )
207         RETURN
208      END IF
209*
210*     Quick return if possible
211*
212      IF( N.EQ.0 )
213     $   RETURN
214*
215*     Determine location of first number in second half.
216*
217      MID = N / 2 + 1
218*
219*     Gather last/first rows of appropriate eigenblocks into center of Z
220*
221      PTR = 1
222*
223*     Determine location of lowest level subproblem in the full storage
224*     scheme
225*
226      CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
227*
228*     Determine size of these matrices.  We add HALF to the value of
229*     the SQRT in case the machine underestimates one of these square
230*     roots.
231*
232      BSIZ1 = INT( HALF+SQRT( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
233      BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
234      DO 10 K = 1, MID - BSIZ1 - 1
235         Z( K ) = ZERO
236   10 CONTINUE
237      CALL SCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
238     $            Z( MID-BSIZ1 ), 1 )
239      CALL SCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
240      DO 20 K = MID + BSIZ2, N
241         Z( K ) = ZERO
242   20 CONTINUE
243*
244*     Loop through remaining levels 1 -> CURLVL applying the Givens
245*     rotations and permutation and then multiplying the center matrices
246*     against the current Z.
247*
248      PTR = 2**TLVLS + 1
249      DO 70 K = 1, CURLVL - 1
250         CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
251         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
252         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
253         ZPTR1 = MID - PSIZ1
254*
255*       Apply Givens at CURR and CURR+1
256*
257         DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
258            CALL SROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
259     $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
260     $                 GIVNUM( 2, I ) )
261   30    CONTINUE
262         DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
263            CALL SROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
264     $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
265     $                 GIVNUM( 2, I ) )
266   40    CONTINUE
267         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
268         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
269         DO 50 I = 0, PSIZ1 - 1
270            ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
271   50    CONTINUE
272         DO 60 I = 0, PSIZ2 - 1
273            ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
274   60    CONTINUE
275*
276*        Multiply Blocks at CURR and CURR+1
277*
278*        Determine size of these matrices.  We add HALF to the value of
279*        the SQRT in case the machine underestimates one of these
280*        square roots.
281*
282         BSIZ1 = INT( HALF+SQRT( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
283         BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+
284     $           1 ) ) ) )
285         IF( BSIZ1.GT.0 ) THEN
286            CALL SGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
287     $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
288         END IF
289         CALL SCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
290     $               1 )
291         IF( BSIZ2.GT.0 ) THEN
292            CALL SGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
293     $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
294         END IF
295         CALL SCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
296     $               Z( MID+BSIZ2 ), 1 )
297*
298         PTR = PTR + 2**( TLVLS-K )
299   70 CONTINUE
300*
301      RETURN
302*
303*     End of SLAEDA
304*
305      END
306