1 SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
2 $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
3 *
4 * -- LAPACK routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
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 *
157 * Loop through remaining levels 1 -> CURLVL applying the Givens
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
2 $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
3 *
4 * -- LAPACK routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
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 *
157 * Loop through remaining levels 1 -> CURLVL applying the Givens
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