Annotation of rpl/lapack/lapack/dlaeda.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
2: $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
3: *
1.5 bertrand 4: * -- LAPACK routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 7: * June 2010
1.1 bertrand 8: *
9: * .. Scalar Arguments ..
10: INTEGER CURLVL, CURPBM, INFO, N, TLVLS
11: * ..
12: * .. Array Arguments ..
13: INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
14: $ PRMPTR( * ), QPTR( * )
15: DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLAEDA computes the Z vector corresponding to the merge step in the
22: * CURLVLth step of the merge process with TLVLS steps for the CURPBMth
23: * problem.
24: *
25: * Arguments
26: * =========
27: *
28: * N (input) INTEGER
29: * The dimension of the symmetric tridiagonal matrix. N >= 0.
30: *
31: * TLVLS (input) INTEGER
32: * The total number of merging levels in the overall divide and
33: * conquer tree.
34: *
35: * CURLVL (input) INTEGER
36: * The current level in the overall merge routine,
37: * 0 <= curlvl <= tlvls.
38: *
39: * CURPBM (input) INTEGER
40: * The current problem in the current level in the overall
41: * merge routine (counting from upper left to lower right).
42: *
43: * PRMPTR (input) INTEGER array, dimension (N lg N)
44: * Contains a list of pointers which indicate where in PERM a
45: * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
46: * indicates the size of the permutation and incidentally the
47: * size of the full, non-deflated problem.
48: *
49: * PERM (input) INTEGER array, dimension (N lg N)
50: * Contains the permutations (from deflation and sorting) to be
51: * applied to each eigenblock.
52: *
53: * GIVPTR (input) INTEGER array, dimension (N lg N)
54: * Contains a list of pointers which indicate where in GIVCOL a
55: * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
56: * indicates the number of Givens rotations.
57: *
58: * GIVCOL (input) INTEGER array, dimension (2, N lg N)
59: * Each pair of numbers indicates a pair of columns to take place
60: * in a Givens rotation.
61: *
62: * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
63: * Each number indicates the S value to be used in the
64: * corresponding Givens rotation.
65: *
66: * Q (input) DOUBLE PRECISION array, dimension (N**2)
67: * Contains the square eigenblocks from previous levels, the
68: * starting positions for blocks are given by QPTR.
69: *
70: * QPTR (input) INTEGER array, dimension (N+2)
71: * Contains a list of pointers which indicate where in Q an
72: * eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates
73: * the size of the block.
74: *
75: * Z (output) DOUBLE PRECISION array, dimension (N)
76: * On output this vector contains the updating vector (the last
77: * row of the first sub-eigenvector matrix and the first row of
78: * the second sub-eigenvector matrix).
79: *
80: * ZTEMP (workspace) DOUBLE PRECISION array, dimension (N)
81: *
82: * INFO (output) INTEGER
83: * = 0: successful exit.
84: * < 0: if INFO = -i, the i-th argument had an illegal value.
85: *
86: * Further Details
87: * ===============
88: *
89: * Based on contributions by
90: * Jeff Rutter, Computer Science Division, University of California
91: * at Berkeley, USA
92: *
93: * =====================================================================
94: *
95: * .. Parameters ..
96: DOUBLE PRECISION ZERO, HALF, ONE
97: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
98: * ..
99: * .. Local Scalars ..
100: INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
101: $ PTR, ZPTR1
102: * ..
103: * .. External Subroutines ..
104: EXTERNAL DCOPY, DGEMV, DROT, XERBLA
105: * ..
106: * .. Intrinsic Functions ..
107: INTRINSIC DBLE, INT, 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( 'DLAEDA', -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( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
146: BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
147: DO 10 K = 1, MID - BSIZ1 - 1
148: Z( K ) = ZERO
149: 10 CONTINUE
150: CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
151: $ Z( MID-BSIZ1 ), 1 )
152: CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
153: DO 20 K = MID + BSIZ2, N
154: Z( K ) = ZERO
155: 20 CONTINUE
156: *
1.5 bertrand 157: * Loop through remaining levels 1 -> CURLVL applying the Givens
1.1 bertrand 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 DROT( 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 DROT( 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( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
196: BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
197: $ 1 ) ) ) )
198: IF( BSIZ1.GT.0 ) THEN
199: CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
200: $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
201: END IF
202: CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
203: $ 1 )
204: IF( BSIZ2.GT.0 ) THEN
205: CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
206: $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
207: END IF
208: CALL DCOPY( 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 DLAEDA
217: *
218: 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