Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / dlaed1.f
1 *> \brief \b DLAED1 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAED1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE DLAED1( N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK,
22 *                          INFO )
23 *
24 *       .. Scalar Arguments ..
25 *       INTEGER            CUTPNT, INFO, LDQ, N
26 *       DOUBLE PRECISION   RHO
27 *       ..
28 *       .. Array Arguments ..
29 *       INTEGER            INDXQ( * ), IWORK( * )
30 *       DOUBLE PRECISION   D( * ), Q( LDQ, * ), WORK( * )
31 *       ..
32 *
33 *
34 *> \par Purpose:
35 *  =============
36 *>
37 *> \verbatim
38 *>
39 *> DLAED1 computes the updated eigensystem of a diagonal
40 *> matrix after modification by a rank-one symmetric matrix.  This
41 *> routine is used only for the eigenproblem which requires all
42 *> eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles
43 *> the case in which eigenvalues only or eigenvalues and eigenvectors
44 *> of a full symmetric matrix (which was reduced to tridiagonal form)
45 *> are desired.
46 *>
47 *>   T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
48 *>
49 *>    where Z = Q**T*u, u is a vector of length N with ones in the
50 *>    CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
51 *>
52 *>    The eigenvectors of the original matrix are stored in Q, and the
53 *>    eigenvalues are in D.  The algorithm consists of three stages:
54 *>
55 *>       The first stage consists of deflating the size of the problem
56 *>       when there are multiple eigenvalues or if there is a zero in
57 *>       the Z vector.  For each such occurrence the dimension of the
58 *>       secular equation problem is reduced by one.  This stage is
59 *>       performed by the routine DLAED2.
60 *>
61 *>       The second stage consists of calculating the updated
62 *>       eigenvalues. This is done by finding the roots of the secular
63 *>       equation via the routine DLAED4 (as called by DLAED3).
64 *>       This routine also calculates the eigenvectors of the current
65 *>       problem.
66 *>
67 *>       The final stage consists of computing the updated eigenvectors
68 *>       directly using the updated eigenvalues.  The eigenvectors for
69 *>       the current problem are multiplied with the eigenvectors from
70 *>       the overall problem.
71 *> \endverbatim
72 *
73 *  Arguments:
74 *  ==========
75 *
76 *> \param[in] N
77 *> \verbatim
78 *>          N is INTEGER
79 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in,out] D
83 *> \verbatim
84 *>          D is DOUBLE PRECISION array, dimension (N)
85 *>         On entry, the eigenvalues of the rank-1-perturbed matrix.
86 *>         On exit, the eigenvalues of the repaired matrix.
87 *> \endverbatim
88 *>
89 *> \param[in,out] Q
90 *> \verbatim
91 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
92 *>         On entry, the eigenvectors of the rank-1-perturbed matrix.
93 *>         On exit, the eigenvectors of the repaired tridiagonal matrix.
94 *> \endverbatim
95 *>
96 *> \param[in] LDQ
97 *> \verbatim
98 *>          LDQ is INTEGER
99 *>         The leading dimension of the array Q.  LDQ >= max(1,N).
100 *> \endverbatim
101 *>
102 *> \param[in,out] INDXQ
103 *> \verbatim
104 *>          INDXQ is INTEGER array, dimension (N)
105 *>         On entry, the permutation which separately sorts the two
106 *>         subproblems in D into ascending order.
107 *>         On exit, the permutation which will reintegrate the
108 *>         subproblems back into sorted order,
109 *>         i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
110 *> \endverbatim
111 *>
112 *> \param[in] RHO
113 *> \verbatim
114 *>          RHO is DOUBLE PRECISION
115 *>         The subdiagonal entry used to create the rank-1 modification.
116 *> \endverbatim
117 *>
118 *> \param[in] CUTPNT
119 *> \verbatim
120 *>          CUTPNT is INTEGER
121 *>         The location of the last eigenvalue in the leading sub-matrix.
122 *>         min(1,N) <= CUTPNT <= N/2.
123 *> \endverbatim
124 *>
125 *> \param[out] WORK
126 *> \verbatim
127 *>          WORK is DOUBLE PRECISION array, dimension (4*N + N**2)
128 *> \endverbatim
129 *>
130 *> \param[out] IWORK
131 *> \verbatim
132 *>          IWORK is INTEGER array, dimension (4*N)
133 *> \endverbatim
134 *>
135 *> \param[out] INFO
136 *> \verbatim
137 *>          INFO is INTEGER
138 *>          = 0:  successful exit.
139 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
140 *>          > 0:  if INFO = 1, an eigenvalue did not converge
141 *> \endverbatim
142 *
143 *  Authors:
144 *  ========
145 *
146 *> \author Univ. of Tennessee
147 *> \author Univ. of California Berkeley
148 *> \author Univ. of Colorado Denver
149 *> \author NAG Ltd.
150 *
151 *> \date June 2016
152 *
153 *> \ingroup auxOTHERcomputational
154 *
155 *> \par Contributors:
156 *  ==================
157 *>
158 *> Jeff Rutter, Computer Science Division, University of California
159 *> at Berkeley, USA \n
160 *>  Modified by Francoise Tisseur, University of Tennessee
161 *>
162 *  =====================================================================
163       SUBROUTINE DLAED1( N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK,
164      $                   INFO )
165 *
166 *  -- LAPACK computational routine (version 3.6.1) --
167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 *     June 2016
170 *
171 *     .. Scalar Arguments ..
172       INTEGER            CUTPNT, INFO, LDQ, N
173       DOUBLE PRECISION   RHO
174 *     ..
175 *     .. Array Arguments ..
176       INTEGER            INDXQ( * ), IWORK( * )
177       DOUBLE PRECISION   D( * ), Q( LDQ, * ), WORK( * )
178 *     ..
179 *
180 *  =====================================================================
181 *
182 *     .. Local Scalars ..
183       INTEGER            COLTYP, I, IDLMDA, INDX, INDXC, INDXP, IQ2, IS,
184      $                   IW, IZ, K, N1, N2, ZPP1
185 *     ..
186 *     .. External Subroutines ..
187       EXTERNAL           DCOPY, DLAED2, DLAED3, DLAMRG, XERBLA
188 *     ..
189 *     .. Intrinsic Functions ..
190       INTRINSIC          MAX, MIN
191 *     ..
192 *     .. Executable Statements ..
193 *
194 *     Test the input parameters.
195 *
196       INFO = 0
197 *
198       IF( N.LT.0 ) THEN
199          INFO = -1
200       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
201          INFO = -4
202       ELSE IF( MIN( 1, N / 2 ).GT.CUTPNT .OR. ( N / 2 ).LT.CUTPNT ) THEN
203          INFO = -7
204       END IF
205       IF( INFO.NE.0 ) THEN
206          CALL XERBLA( 'DLAED1', -INFO )
207          RETURN
208       END IF
209 *
210 *     Quick return if possible
211 *
212       IF( N.EQ.0 )
213      $   RETURN
214 *
215 *     The following values are integer pointers which indicate
216 *     the portion of the workspace
217 *     used by a particular array in DLAED2 and DLAED3.
218 *
219       IZ = 1
220       IDLMDA = IZ + N
221       IW = IDLMDA + N
222       IQ2 = IW + N
223 *
224       INDX = 1
225       INDXC = INDX + N
226       COLTYP = INDXC + N
227       INDXP = COLTYP + N
228 *
229 *
230 *     Form the z-vector which consists of the last row of Q_1 and the
231 *     first row of Q_2.
232 *
233       CALL DCOPY( CUTPNT, Q( CUTPNT, 1 ), LDQ, WORK( IZ ), 1 )
234       ZPP1 = CUTPNT + 1
235       CALL DCOPY( N-CUTPNT, Q( ZPP1, ZPP1 ), LDQ, WORK( IZ+CUTPNT ), 1 )
236 *
237 *     Deflate eigenvalues.
238 *
239       CALL DLAED2( K, N, CUTPNT, D, Q, LDQ, INDXQ, RHO, WORK( IZ ),
240      $             WORK( IDLMDA ), WORK( IW ), WORK( IQ2 ),
241      $             IWORK( INDX ), IWORK( INDXC ), IWORK( INDXP ),
242      $             IWORK( COLTYP ), INFO )
243 *
244       IF( INFO.NE.0 )
245      $   GO TO 20
246 *
247 *     Solve Secular Equation.
248 *
249       IF( K.NE.0 ) THEN
250          IS = ( IWORK( COLTYP )+IWORK( COLTYP+1 ) )*CUTPNT +
251      $        ( IWORK( COLTYP+1 )+IWORK( COLTYP+2 ) )*( N-CUTPNT ) + IQ2
252          CALL DLAED3( K, N, CUTPNT, D, Q, LDQ, RHO, WORK( IDLMDA ),
253      $                WORK( IQ2 ), IWORK( INDXC ), IWORK( COLTYP ),
254      $                WORK( IW ), WORK( IS ), INFO )
255          IF( INFO.NE.0 )
256      $      GO TO 20
257 *
258 *     Prepare the INDXQ sorting permutation.
259 *
260          N1 = K
261          N2 = N - K
262          CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
263       ELSE
264          DO 10 I = 1, N
265             INDXQ( I ) = I
266    10    CONTINUE
267       END IF
268 *
269    20 CONTINUE
270       RETURN
271 *
272 *     End of DLAED1
273 *
274       END