LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaeda.f
Go to the documentation of this file.
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 *> \date September 2012
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 slaeda( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
167  $ givcol, givnum, q, qptr, z, ztemp, info )
168 *
169 * -- LAPACK computational routine (version 3.4.2) --
170 * -- LAPACK is a software package provided by Univ. of Tennessee, --
171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172 * September 2012
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  REAL givnum( 2, * ), q( * ), z( * ), ztemp( * )
181 * ..
182 *
183 * =====================================================================
184 *
185 * .. Parameters ..
186  REAL zero, half, one
187  parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0 )
188 * ..
189 * .. Local Scalars ..
190  INTEGER bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2,
191  $ ptr, zptr1
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL scopy, sgemv, srot, xerbla
195 * ..
196 * .. Intrinsic Functions ..
197  INTRINSIC int, REAL, 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( 'SLAEDA', -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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
236  bsiz2 = int( half+sqrt( REAL( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
237  DO 10 k = 1, mid - bsiz1 - 1
238  z( k ) = zero
239  10 continue
240  CALL scopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,
241  $ z( mid-bsiz1 ), 1 )
242  CALL scopy( 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 srot( 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 srot( 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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
286  bsiz2 = int( half+sqrt( REAL( QPTR( CURR+2 )-QPTR( CURR+ $ 1 ) ) ) )
287  IF( bsiz1.GT.0 ) THEN
288  CALL sgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),
289  $ bsiz1, ztemp( 1 ), 1, zero, z( zptr1 ), 1 )
290  END IF
291  CALL scopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1, z( zptr1+bsiz1 ),
292  $ 1 )
293  IF( bsiz2.GT.0 ) THEN
294  CALL sgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),
295  $ bsiz2, ztemp( psiz1+1 ), 1, zero, z( mid ), 1 )
296  END IF
297  CALL scopy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1,
298  $ z( mid+bsiz2 ), 1 )
299 *
300  ptr = ptr + 2**( tlvls-k )
301  70 continue
302 *
303  return
304 *
305 * End of SLAEDA
306 *
307  END
308