1 SUBROUTINE DLAED1( N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.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 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER CUTPNT, INFO, LDQ, N
11 DOUBLE PRECISION RHO
12 * ..
13 * .. Array Arguments ..
14 INTEGER INDXQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), Q( LDQ, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAED1 computes the updated eigensystem of a diagonal
22 * matrix after modification by a rank-one symmetric matrix. This
23 * routine is used only for the eigenproblem which requires all
24 * eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles
25 * the case in which eigenvalues only or eigenvalues and eigenvectors
26 * of a full symmetric matrix (which was reduced to tridiagonal form)
27 * are desired.
28 *
29 * T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
30 *
31 * where Z = Q**T*u, u is a vector of length N with ones in the
32 * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
33 *
34 * The eigenvectors of the original matrix are stored in Q, and the
35 * eigenvalues are in D. The algorithm consists of three stages:
36 *
37 * The first stage consists of deflating the size of the problem
38 * when there are multiple eigenvalues or if there is a zero in
39 * the Z vector. For each such occurence the dimension of the
40 * secular equation problem is reduced by one. This stage is
41 * performed by the routine DLAED2.
42 *
43 * The second stage consists of calculating the updated
44 * eigenvalues. This is done by finding the roots of the secular
45 * equation via the routine DLAED4 (as called by DLAED3).
46 * This routine also calculates the eigenvectors of the current
47 * problem.
48 *
49 * The final stage consists of computing the updated eigenvectors
50 * directly using the updated eigenvalues. The eigenvectors for
51 * the current problem are multiplied with the eigenvectors from
52 * the overall problem.
53 *
54 * Arguments
55 * =========
56 *
57 * N (input) INTEGER
58 * The dimension of the symmetric tridiagonal matrix. N >= 0.
59 *
60 * D (input/output) DOUBLE PRECISION array, dimension (N)
61 * On entry, the eigenvalues of the rank-1-perturbed matrix.
62 * On exit, the eigenvalues of the repaired matrix.
63 *
64 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
65 * On entry, the eigenvectors of the rank-1-perturbed matrix.
66 * On exit, the eigenvectors of the repaired tridiagonal matrix.
67 *
68 * LDQ (input) INTEGER
69 * The leading dimension of the array Q. LDQ >= max(1,N).
70 *
71 * INDXQ (input/output) INTEGER array, dimension (N)
72 * On entry, the permutation which separately sorts the two
73 * subproblems in D into ascending order.
74 * On exit, the permutation which will reintegrate the
75 * subproblems back into sorted order,
76 * i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
77 *
78 * RHO (input) DOUBLE PRECISION
79 * The subdiagonal entry used to create the rank-1 modification.
80 *
81 * CUTPNT (input) INTEGER
82 * The location of the last eigenvalue in the leading sub-matrix.
83 * min(1,N) <= CUTPNT <= N/2.
84 *
85 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)
86 *
87 * IWORK (workspace) INTEGER array, dimension (4*N)
88 *
89 * INFO (output) INTEGER
90 * = 0: successful exit.
91 * < 0: if INFO = -i, the i-th argument had an illegal value.
92 * > 0: if INFO = 1, an eigenvalue did not converge
93 *
94 * Further Details
95 * ===============
96 *
97 * Based on contributions by
98 * Jeff Rutter, Computer Science Division, University of California
99 * at Berkeley, USA
100 * Modified by Francoise Tisseur, University of Tennessee.
101 *
102 * =====================================================================
103 *
104 * .. Local Scalars ..
105 INTEGER COLTYP, I, IDLMDA, INDX, INDXC, INDXP, IQ2, IS,
106 $ IW, IZ, K, N1, N2, ZPP1
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL DCOPY, DLAED2, DLAED3, DLAMRG, XERBLA
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC MAX, MIN
113 * ..
114 * .. Executable Statements ..
115 *
116 * Test the input parameters.
117 *
118 INFO = 0
119 *
120 IF( N.LT.0 ) THEN
121 INFO = -1
122 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
123 INFO = -4
124 ELSE IF( MIN( 1, N / 2 ).GT.CUTPNT .OR. ( N / 2 ).LT.CUTPNT ) THEN
125 INFO = -7
126 END IF
127 IF( INFO.NE.0 ) THEN
128 CALL XERBLA( 'DLAED1', -INFO )
129 RETURN
130 END IF
131 *
132 * Quick return if possible
133 *
134 IF( N.EQ.0 )
135 $ RETURN
136 *
137 * The following values are integer pointers which indicate
138 * the portion of the workspace
139 * used by a particular array in DLAED2 and DLAED3.
140 *
141 IZ = 1
142 IDLMDA = IZ + N
143 IW = IDLMDA + N
144 IQ2 = IW + N
145 *
146 INDX = 1
147 INDXC = INDX + N
148 COLTYP = INDXC + N
149 INDXP = COLTYP + N
150 *
151 *
152 * Form the z-vector which consists of the last row of Q_1 and the
153 * first row of Q_2.
154 *
155 CALL DCOPY( CUTPNT, Q( CUTPNT, 1 ), LDQ, WORK( IZ ), 1 )
156 ZPP1 = CUTPNT + 1
157 CALL DCOPY( N-CUTPNT, Q( ZPP1, ZPP1 ), LDQ, WORK( IZ+CUTPNT ), 1 )
158 *
159 * Deflate eigenvalues.
160 *
161 CALL DLAED2( K, N, CUTPNT, D, Q, LDQ, INDXQ, RHO, WORK( IZ ),
162 $ WORK( IDLMDA ), WORK( IW ), WORK( IQ2 ),
163 $ IWORK( INDX ), IWORK( INDXC ), IWORK( INDXP ),
164 $ IWORK( COLTYP ), INFO )
165 *
166 IF( INFO.NE.0 )
167 $ GO TO 20
168 *
169 * Solve Secular Equation.
170 *
171 IF( K.NE.0 ) THEN
172 IS = ( IWORK( COLTYP )+IWORK( COLTYP+1 ) )*CUTPNT +
173 $ ( IWORK( COLTYP+1 )+IWORK( COLTYP+2 ) )*( N-CUTPNT ) + IQ2
174 CALL DLAED3( K, N, CUTPNT, D, Q, LDQ, RHO, WORK( IDLMDA ),
175 $ WORK( IQ2 ), IWORK( INDXC ), IWORK( COLTYP ),
176 $ WORK( IW ), WORK( IS ), INFO )
177 IF( INFO.NE.0 )
178 $ GO TO 20
179 *
180 * Prepare the INDXQ sorting permutation.
181 *
182 N1 = K
183 N2 = N - K
184 CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
185 ELSE
186 DO 10 I = 1, N
187 INDXQ( I ) = I
188 10 CONTINUE
189 END IF
190 *
191 20 CONTINUE
192 RETURN
193 *
194 * End of DLAED1
195 *
196 END
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.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 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER CUTPNT, INFO, LDQ, N
11 DOUBLE PRECISION RHO
12 * ..
13 * .. Array Arguments ..
14 INTEGER INDXQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), Q( LDQ, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAED1 computes the updated eigensystem of a diagonal
22 * matrix after modification by a rank-one symmetric matrix. This
23 * routine is used only for the eigenproblem which requires all
24 * eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles
25 * the case in which eigenvalues only or eigenvalues and eigenvectors
26 * of a full symmetric matrix (which was reduced to tridiagonal form)
27 * are desired.
28 *
29 * T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
30 *
31 * where Z = Q**T*u, u is a vector of length N with ones in the
32 * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
33 *
34 * The eigenvectors of the original matrix are stored in Q, and the
35 * eigenvalues are in D. The algorithm consists of three stages:
36 *
37 * The first stage consists of deflating the size of the problem
38 * when there are multiple eigenvalues or if there is a zero in
39 * the Z vector. For each such occurence the dimension of the
40 * secular equation problem is reduced by one. This stage is
41 * performed by the routine DLAED2.
42 *
43 * The second stage consists of calculating the updated
44 * eigenvalues. This is done by finding the roots of the secular
45 * equation via the routine DLAED4 (as called by DLAED3).
46 * This routine also calculates the eigenvectors of the current
47 * problem.
48 *
49 * The final stage consists of computing the updated eigenvectors
50 * directly using the updated eigenvalues. The eigenvectors for
51 * the current problem are multiplied with the eigenvectors from
52 * the overall problem.
53 *
54 * Arguments
55 * =========
56 *
57 * N (input) INTEGER
58 * The dimension of the symmetric tridiagonal matrix. N >= 0.
59 *
60 * D (input/output) DOUBLE PRECISION array, dimension (N)
61 * On entry, the eigenvalues of the rank-1-perturbed matrix.
62 * On exit, the eigenvalues of the repaired matrix.
63 *
64 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
65 * On entry, the eigenvectors of the rank-1-perturbed matrix.
66 * On exit, the eigenvectors of the repaired tridiagonal matrix.
67 *
68 * LDQ (input) INTEGER
69 * The leading dimension of the array Q. LDQ >= max(1,N).
70 *
71 * INDXQ (input/output) INTEGER array, dimension (N)
72 * On entry, the permutation which separately sorts the two
73 * subproblems in D into ascending order.
74 * On exit, the permutation which will reintegrate the
75 * subproblems back into sorted order,
76 * i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
77 *
78 * RHO (input) DOUBLE PRECISION
79 * The subdiagonal entry used to create the rank-1 modification.
80 *
81 * CUTPNT (input) INTEGER
82 * The location of the last eigenvalue in the leading sub-matrix.
83 * min(1,N) <= CUTPNT <= N/2.
84 *
85 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)
86 *
87 * IWORK (workspace) INTEGER array, dimension (4*N)
88 *
89 * INFO (output) INTEGER
90 * = 0: successful exit.
91 * < 0: if INFO = -i, the i-th argument had an illegal value.
92 * > 0: if INFO = 1, an eigenvalue did not converge
93 *
94 * Further Details
95 * ===============
96 *
97 * Based on contributions by
98 * Jeff Rutter, Computer Science Division, University of California
99 * at Berkeley, USA
100 * Modified by Francoise Tisseur, University of Tennessee.
101 *
102 * =====================================================================
103 *
104 * .. Local Scalars ..
105 INTEGER COLTYP, I, IDLMDA, INDX, INDXC, INDXP, IQ2, IS,
106 $ IW, IZ, K, N1, N2, ZPP1
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL DCOPY, DLAED2, DLAED3, DLAMRG, XERBLA
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC MAX, MIN
113 * ..
114 * .. Executable Statements ..
115 *
116 * Test the input parameters.
117 *
118 INFO = 0
119 *
120 IF( N.LT.0 ) THEN
121 INFO = -1
122 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
123 INFO = -4
124 ELSE IF( MIN( 1, N / 2 ).GT.CUTPNT .OR. ( N / 2 ).LT.CUTPNT ) THEN
125 INFO = -7
126 END IF
127 IF( INFO.NE.0 ) THEN
128 CALL XERBLA( 'DLAED1', -INFO )
129 RETURN
130 END IF
131 *
132 * Quick return if possible
133 *
134 IF( N.EQ.0 )
135 $ RETURN
136 *
137 * The following values are integer pointers which indicate
138 * the portion of the workspace
139 * used by a particular array in DLAED2 and DLAED3.
140 *
141 IZ = 1
142 IDLMDA = IZ + N
143 IW = IDLMDA + N
144 IQ2 = IW + N
145 *
146 INDX = 1
147 INDXC = INDX + N
148 COLTYP = INDXC + N
149 INDXP = COLTYP + N
150 *
151 *
152 * Form the z-vector which consists of the last row of Q_1 and the
153 * first row of Q_2.
154 *
155 CALL DCOPY( CUTPNT, Q( CUTPNT, 1 ), LDQ, WORK( IZ ), 1 )
156 ZPP1 = CUTPNT + 1
157 CALL DCOPY( N-CUTPNT, Q( ZPP1, ZPP1 ), LDQ, WORK( IZ+CUTPNT ), 1 )
158 *
159 * Deflate eigenvalues.
160 *
161 CALL DLAED2( K, N, CUTPNT, D, Q, LDQ, INDXQ, RHO, WORK( IZ ),
162 $ WORK( IDLMDA ), WORK( IW ), WORK( IQ2 ),
163 $ IWORK( INDX ), IWORK( INDXC ), IWORK( INDXP ),
164 $ IWORK( COLTYP ), INFO )
165 *
166 IF( INFO.NE.0 )
167 $ GO TO 20
168 *
169 * Solve Secular Equation.
170 *
171 IF( K.NE.0 ) THEN
172 IS = ( IWORK( COLTYP )+IWORK( COLTYP+1 ) )*CUTPNT +
173 $ ( IWORK( COLTYP+1 )+IWORK( COLTYP+2 ) )*( N-CUTPNT ) + IQ2
174 CALL DLAED3( K, N, CUTPNT, D, Q, LDQ, RHO, WORK( IDLMDA ),
175 $ WORK( IQ2 ), IWORK( INDXC ), IWORK( COLTYP ),
176 $ WORK( IW ), WORK( IS ), INFO )
177 IF( INFO.NE.0 )
178 $ GO TO 20
179 *
180 * Prepare the INDXQ sorting permutation.
181 *
182 N1 = K
183 N2 = N - K
184 CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
185 ELSE
186 DO 10 I = 1, N
187 INDXQ( I ) = I
188 10 CONTINUE
189 END IF
190 *
191 20 CONTINUE
192 RETURN
193 *
194 * End of DLAED1
195 *
196 END