1*> \brief \b ZLAED0 used by ZSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLAED0 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed0.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed0.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed0.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
22*                          IWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, LDQ, LDQS, N, QSIZ
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IWORK( * )
29*       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
30*       COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> Using the divide and conquer method, ZLAED0 computes all eigenvalues
40*> of a symmetric tridiagonal matrix which is one diagonal block of
41*> those from reducing a dense or band Hermitian matrix and
42*> corresponding eigenvectors of the dense or band matrix.
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in] QSIZ
49*> \verbatim
50*>          QSIZ is INTEGER
51*>         The dimension of the unitary matrix used to reduce
52*>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*>          N is INTEGER
58*>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
59*> \endverbatim
60*>
61*> \param[in,out] D
62*> \verbatim
63*>          D is DOUBLE PRECISION array, dimension (N)
64*>         On entry, the diagonal elements of the tridiagonal matrix.
65*>         On exit, the eigenvalues in ascending order.
66*> \endverbatim
67*>
68*> \param[in,out] E
69*> \verbatim
70*>          E is DOUBLE PRECISION array, dimension (N-1)
71*>         On entry, the off-diagonal elements of the tridiagonal matrix.
72*>         On exit, E has been destroyed.
73*> \endverbatim
74*>
75*> \param[in,out] Q
76*> \verbatim
77*>          Q is COMPLEX*16 array, dimension (LDQ,N)
78*>         On entry, Q must contain an QSIZ x N matrix whose columns
79*>         unitarily orthonormal. It is a part of the unitary matrix
80*>         that reduces the full dense Hermitian matrix to a
81*>         (reducible) symmetric tridiagonal matrix.
82*> \endverbatim
83*>
84*> \param[in] LDQ
85*> \verbatim
86*>          LDQ is INTEGER
87*>         The leading dimension of the array Q.  LDQ >= max(1,N).
88*> \endverbatim
89*>
90*> \param[out] IWORK
91*> \verbatim
92*>          IWORK is INTEGER array,
93*>         the dimension of IWORK must be at least
94*>                      6 + 6*N + 5*N*lg N
95*>                      ( lg( N ) = smallest integer k
96*>                                  such that 2^k >= N )
97*> \endverbatim
98*>
99*> \param[out] RWORK
100*> \verbatim
101*>          RWORK is DOUBLE PRECISION array,
102*>                               dimension (1 + 3*N + 2*N*lg N + 3*N**2)
103*>                        ( lg( N ) = smallest integer k
104*>                                    such that 2^k >= N )
105*> \endverbatim
106*>
107*> \param[out] QSTORE
108*> \verbatim
109*>          QSTORE is COMPLEX*16 array, dimension (LDQS, N)
110*>         Used to store parts of
111*>         the eigenvector matrix when the updating matrix multiplies
112*>         take place.
113*> \endverbatim
114*>
115*> \param[in] LDQS
116*> \verbatim
117*>          LDQS is INTEGER
118*>         The leading dimension of the array QSTORE.
119*>         LDQS >= max(1,N).
120*> \endverbatim
121*>
122*> \param[out] INFO
123*> \verbatim
124*>          INFO is INTEGER
125*>          = 0:  successful exit.
126*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
127*>          > 0:  The algorithm failed to compute an eigenvalue while
128*>                working on the submatrix lying in rows and columns
129*>                INFO/(N+1) through mod(INFO,N+1).
130*> \endverbatim
131*
132*  Authors:
133*  ========
134*
135*> \author Univ. of Tennessee
136*> \author Univ. of California Berkeley
137*> \author Univ. of Colorado Denver
138*> \author NAG Ltd.
139*
140*> \ingroup complex16OTHERcomputational
141*
142*  =====================================================================
143      SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
144     $                   IWORK, INFO )
145*
146*  -- LAPACK computational routine --
147*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
148*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149*
150*     .. Scalar Arguments ..
151      INTEGER            INFO, LDQ, LDQS, N, QSIZ
152*     ..
153*     .. Array Arguments ..
154      INTEGER            IWORK( * )
155      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
156      COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
157*     ..
158*
159*  =====================================================================
160*
161*  Warning:      N could be as big as QSIZ!
162*
163*     .. Parameters ..
164      DOUBLE PRECISION   TWO
165      PARAMETER          ( TWO = 2.D+0 )
166*     ..
167*     .. Local Scalars ..
168      INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
169     $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
170     $                   J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
171     $                   SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
172      DOUBLE PRECISION   TEMP
173*     ..
174*     .. External Subroutines ..
175      EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
176*     ..
177*     .. External Functions ..
178      INTEGER            ILAENV
179      EXTERNAL           ILAENV
180*     ..
181*     .. Intrinsic Functions ..
182      INTRINSIC          ABS, DBLE, INT, LOG, MAX
183*     ..
184*     .. Executable Statements ..
185*
186*     Test the input parameters.
187*
188      INFO = 0
189*
190*     IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
191*        INFO = -1
192*     ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
193*    $        THEN
194      IF( QSIZ.LT.MAX( 0, N ) ) THEN
195         INFO = -1
196      ELSE IF( N.LT.0 ) THEN
197         INFO = -2
198      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
199         INFO = -6
200      ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
201         INFO = -8
202      END IF
203      IF( INFO.NE.0 ) THEN
204         CALL XERBLA( 'ZLAED0', -INFO )
205         RETURN
206      END IF
207*
208*     Quick return if possible
209*
210      IF( N.EQ.0 )
211     $   RETURN
212*
213      SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
214*
215*     Determine the size and placement of the submatrices, and save in
216*     the leading elements of IWORK.
217*
218      IWORK( 1 ) = N
219      SUBPBS = 1
220      TLVLS = 0
221   10 CONTINUE
222      IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
223         DO 20 J = SUBPBS, 1, -1
224            IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
225            IWORK( 2*J-1 ) = IWORK( J ) / 2
226   20    CONTINUE
227         TLVLS = TLVLS + 1
228         SUBPBS = 2*SUBPBS
229         GO TO 10
230      END IF
231      DO 30 J = 2, SUBPBS
232         IWORK( J ) = IWORK( J ) + IWORK( J-1 )
233   30 CONTINUE
234*
235*     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
236*     using rank-1 modifications (cuts).
237*
238      SPM1 = SUBPBS - 1
239      DO 40 I = 1, SPM1
240         SUBMAT = IWORK( I ) + 1
241         SMM1 = SUBMAT - 1
242         D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
243         D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
244   40 CONTINUE
245*
246      INDXQ = 4*N + 3
247*
248*     Set up workspaces for eigenvalues only/accumulate new vectors
249*     routine
250*
251      TEMP = LOG( DBLE( N ) ) / LOG( TWO )
252      LGN = INT( TEMP )
253      IF( 2**LGN.LT.N )
254     $   LGN = LGN + 1
255      IF( 2**LGN.LT.N )
256     $   LGN = LGN + 1
257      IPRMPT = INDXQ + N + 1
258      IPERM = IPRMPT + N*LGN
259      IQPTR = IPERM + N*LGN
260      IGIVPT = IQPTR + N + 2
261      IGIVCL = IGIVPT + N*LGN
262*
263      IGIVNM = 1
264      IQ = IGIVNM + 2*N*LGN
265      IWREM = IQ + N**2 + 1
266*     Initialize pointers
267      DO 50 I = 0, SUBPBS
268         IWORK( IPRMPT+I ) = 1
269         IWORK( IGIVPT+I ) = 1
270   50 CONTINUE
271      IWORK( IQPTR ) = 1
272*
273*     Solve each submatrix eigenproblem at the bottom of the divide and
274*     conquer tree.
275*
276      CURR = 0
277      DO 70 I = 0, SPM1
278         IF( I.EQ.0 ) THEN
279            SUBMAT = 1
280            MATSIZ = IWORK( 1 )
281         ELSE
282            SUBMAT = IWORK( I ) + 1
283            MATSIZ = IWORK( I+1 ) - IWORK( I )
284         END IF
285         LL = IQ - 1 + IWORK( IQPTR+CURR )
286         CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
287     $                RWORK( LL ), MATSIZ, RWORK, INFO )
288         CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
289     $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
290     $                RWORK( IWREM ) )
291         IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
292         CURR = CURR + 1
293         IF( INFO.GT.0 ) THEN
294            INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
295            RETURN
296         END IF
297         K = 1
298         DO 60 J = SUBMAT, IWORK( I+1 )
299            IWORK( INDXQ+J ) = K
300            K = K + 1
301   60    CONTINUE
302   70 CONTINUE
303*
304*     Successively merge eigensystems of adjacent submatrices
305*     into eigensystem for the corresponding larger matrix.
306*
307*     while ( SUBPBS > 1 )
308*
309      CURLVL = 1
310   80 CONTINUE
311      IF( SUBPBS.GT.1 ) THEN
312         SPM2 = SUBPBS - 2
313         DO 90 I = 0, SPM2, 2
314            IF( I.EQ.0 ) THEN
315               SUBMAT = 1
316               MATSIZ = IWORK( 2 )
317               MSD2 = IWORK( 1 )
318               CURPRB = 0
319            ELSE
320               SUBMAT = IWORK( I ) + 1
321               MATSIZ = IWORK( I+2 ) - IWORK( I )
322               MSD2 = MATSIZ / 2
323               CURPRB = CURPRB + 1
324            END IF
325*
326*     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
327*     into an eigensystem of size MATSIZ.  ZLAED7 handles the case
328*     when the eigenvectors of a full or band Hermitian matrix (which
329*     was reduced to tridiagonal form) are desired.
330*
331*     I am free to use Q as a valuable working space until Loop 150.
332*
333            CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
334     $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
335     $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
336     $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
337     $                   IWORK( IPERM ), IWORK( IGIVPT ),
338     $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
339     $                   Q( 1, SUBMAT ), RWORK( IWREM ),
340     $                   IWORK( SUBPBS+1 ), INFO )
341            IF( INFO.GT.0 ) THEN
342               INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
343               RETURN
344            END IF
345            IWORK( I / 2+1 ) = IWORK( I+2 )
346   90    CONTINUE
347         SUBPBS = SUBPBS / 2
348         CURLVL = CURLVL + 1
349         GO TO 80
350      END IF
351*
352*     end while
353*
354*     Re-merge the eigenvalues/vectors which were deflated at the final
355*     merge step.
356*
357      DO 100 I = 1, N
358         J = IWORK( INDXQ+I )
359         RWORK( I ) = D( J )
360         CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
361  100 CONTINUE
362      CALL DCOPY( N, RWORK, 1, D, 1 )
363*
364      RETURN
365*
366*     End of ZLAED0
367*
368      END
369