Annotation of rpl/lapack/lapack/dlaeda.f, revision 1.19
1.19 ! bertrand 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.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
22: * GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
1.16 bertrand 23: *
1.9 bertrand 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: * ..
1.16 bertrand 32: *
1.9 bertrand 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: *
1.16 bertrand 150: *> \author Univ. of Tennessee
151: *> \author Univ. of California Berkeley
152: *> \author Univ. of Colorado Denver
153: *> \author NAG Ltd.
1.9 bertrand 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: * =====================================================================
1.1 bertrand 164: SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
165: $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
166: *
1.19 ! bertrand 167: * -- LAPACK computational routine --
1.1 bertrand 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: *
1.5 bertrand 244: * Loop through remaining levels 1 -> CURLVL applying the Givens
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaeda.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack