LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaeda.f
Go to the documentation of this file.
1*> \brief \b DLAEDA used by DSTEDC. 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 DLAEDA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaeda.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaeda.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaeda.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAEDA( 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* DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLAEDA 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 laeda
156*
157*> \par Contributors:
158* ==================
159*>
160*> Jeff Rutter, Computer Science Division, University of California
161*> at Berkeley, USA
162*
163* =====================================================================
164 SUBROUTINE dlaeda( 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 DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 DOUBLE PRECISION ZERO, HALF, ONE
184 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
185* ..
186* .. Local Scalars ..
187 INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
188 $ ptr, zptr1
189* ..
190* .. External Subroutines ..
191 EXTERNAL dcopy, dgemv, drot, xerbla
192* ..
193* .. Intrinsic Functions ..
194 INTRINSIC dble, int, 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( 'DLAEDA', -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( dble( qptr( curr+1 )-qptr( curr ) ) ) )
233 bsiz2 = int( half+sqrt( dble( qptr( curr+2 )-qptr( curr+1 ) ) ) )
234 DO 10 k = 1, mid - bsiz1 - 1
235 z( k ) = zero
236 10 CONTINUE
237 CALL dcopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,
238 $ z( mid-bsiz1 ), 1 )
239 CALL dcopy( 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 drot( 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 drot( 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( dble( qptr( curr+1 )-qptr( curr ) ) ) )
283 bsiz2 = int( half+sqrt( dble( qptr( curr+2 )-qptr( curr+
284 $ 1 ) ) ) )
285 IF( bsiz1.GT.0 ) THEN
286 CALL dgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),
287 $ bsiz1, ztemp( 1 ), 1, zero, z( zptr1 ), 1 )
288 END IF
289 CALL dcopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1, z( zptr1+bsiz1 ),
290 $ 1 )
291 IF( bsiz2.GT.0 ) THEN
292 CALL dgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),
293 $ bsiz2, ztemp( psiz1+1 ), 1, zero, z( mid ), 1 )
294 END IF
295 CALL dcopy( 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 DLAEDA
304*
305 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dlaeda(n, tlvls, curlvl, curpbm, prmptr, perm, givptr, givcol, givnum, q, qptr, z, ztemp, info)
DLAEDA used by DSTEDC. Computes the Z vector determining the rank-one modification of the diagonal ma...
Definition dlaeda.f:166
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92