Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / zlaed0.f
1 *> \brief \b ZLAED0 used by sstedc. 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 *> \date September 2012
141 *
142 *> \ingroup complex16OTHERcomputational
143 *
144 *  =====================================================================
145       SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
146      $                   IWORK, INFO )
147 *
148 *  -- LAPACK computational routine (version 3.4.2) --
149 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
150 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 *     September 2012
152 *
153 *     .. Scalar Arguments ..
154       INTEGER            INFO, LDQ, LDQS, N, QSIZ
155 *     ..
156 *     .. Array Arguments ..
157       INTEGER            IWORK( * )
158       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
159       COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
160 *     ..
161 *
162 *  =====================================================================
163 *
164 *  Warning:      N could be as big as QSIZ!
165 *
166 *     .. Parameters ..
167       DOUBLE PRECISION   TWO
168       PARAMETER          ( TWO = 2.D+0 )
169 *     ..
170 *     .. Local Scalars ..
171       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
172      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
173      $                   J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
174      $                   SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
175       DOUBLE PRECISION   TEMP
176 *     ..
177 *     .. External Subroutines ..
178       EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
179 *     ..
180 *     .. External Functions ..
181       INTEGER            ILAENV
182       EXTERNAL           ILAENV
183 *     ..
184 *     .. Intrinsic Functions ..
185       INTRINSIC          ABS, DBLE, INT, LOG, MAX
186 *     ..
187 *     .. Executable Statements ..
188 *
189 *     Test the input parameters.
190 *
191       INFO = 0
192 *
193 *     IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
194 *        INFO = -1
195 *     ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
196 *    $        THEN
197       IF( QSIZ.LT.MAX( 0, N ) ) THEN
198          INFO = -1
199       ELSE IF( N.LT.0 ) THEN
200          INFO = -2
201       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
202          INFO = -6
203       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
204          INFO = -8
205       END IF
206       IF( INFO.NE.0 ) THEN
207          CALL XERBLA( 'ZLAED0', -INFO )
208          RETURN
209       END IF
210 *
211 *     Quick return if possible
212 *
213       IF( N.EQ.0 )
214      $   RETURN
215 *
216       SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
217 *
218 *     Determine the size and placement of the submatrices, and save in
219 *     the leading elements of IWORK.
220 *
221       IWORK( 1 ) = N
222       SUBPBS = 1
223       TLVLS = 0
224    10 CONTINUE
225       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
226          DO 20 J = SUBPBS, 1, -1
227             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
228             IWORK( 2*J-1 ) = IWORK( J ) / 2
229    20    CONTINUE
230          TLVLS = TLVLS + 1
231          SUBPBS = 2*SUBPBS
232          GO TO 10
233       END IF
234       DO 30 J = 2, SUBPBS
235          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
236    30 CONTINUE
237 *
238 *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
239 *     using rank-1 modifications (cuts).
240 *
241       SPM1 = SUBPBS - 1
242       DO 40 I = 1, SPM1
243          SUBMAT = IWORK( I ) + 1
244          SMM1 = SUBMAT - 1
245          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
246          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
247    40 CONTINUE
248 *
249       INDXQ = 4*N + 3
250 *
251 *     Set up workspaces for eigenvalues only/accumulate new vectors
252 *     routine
253 *
254       TEMP = LOG( DBLE( N ) ) / LOG( TWO )
255       LGN = INT( TEMP )
256       IF( 2**LGN.LT.N )
257      $   LGN = LGN + 1
258       IF( 2**LGN.LT.N )
259      $   LGN = LGN + 1
260       IPRMPT = INDXQ + N + 1
261       IPERM = IPRMPT + N*LGN
262       IQPTR = IPERM + N*LGN
263       IGIVPT = IQPTR + N + 2
264       IGIVCL = IGIVPT + N*LGN
265 *
266       IGIVNM = 1
267       IQ = IGIVNM + 2*N*LGN
268       IWREM = IQ + N**2 + 1
269 *     Initialize pointers
270       DO 50 I = 0, SUBPBS
271          IWORK( IPRMPT+I ) = 1
272          IWORK( IGIVPT+I ) = 1
273    50 CONTINUE
274       IWORK( IQPTR ) = 1
275 *
276 *     Solve each submatrix eigenproblem at the bottom of the divide and
277 *     conquer tree.
278 *
279       CURR = 0
280       DO 70 I = 0, SPM1
281          IF( I.EQ.0 ) THEN
282             SUBMAT = 1
283             MATSIZ = IWORK( 1 )
284          ELSE
285             SUBMAT = IWORK( I ) + 1
286             MATSIZ = IWORK( I+1 ) - IWORK( I )
287          END IF
288          LL = IQ - 1 + IWORK( IQPTR+CURR )
289          CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
290      $                RWORK( LL ), MATSIZ, RWORK, INFO )
291          CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
292      $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
293      $                RWORK( IWREM ) )
294          IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
295          CURR = CURR + 1
296          IF( INFO.GT.0 ) THEN
297             INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
298             RETURN
299          END IF
300          K = 1
301          DO 60 J = SUBMAT, IWORK( I+1 )
302             IWORK( INDXQ+J ) = K
303             K = K + 1
304    60    CONTINUE
305    70 CONTINUE
306 *
307 *     Successively merge eigensystems of adjacent submatrices
308 *     into eigensystem for the corresponding larger matrix.
309 *
310 *     while ( SUBPBS > 1 )
311 *
312       CURLVL = 1
313    80 CONTINUE
314       IF( SUBPBS.GT.1 ) THEN
315          SPM2 = SUBPBS - 2
316          DO 90 I = 0, SPM2, 2
317             IF( I.EQ.0 ) THEN
318                SUBMAT = 1
319                MATSIZ = IWORK( 2 )
320                MSD2 = IWORK( 1 )
321                CURPRB = 0
322             ELSE
323                SUBMAT = IWORK( I ) + 1
324                MATSIZ = IWORK( I+2 ) - IWORK( I )
325                MSD2 = MATSIZ / 2
326                CURPRB = CURPRB + 1
327             END IF
328 *
329 *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
330 *     into an eigensystem of size MATSIZ.  ZLAED7 handles the case
331 *     when the eigenvectors of a full or band Hermitian matrix (which
332 *     was reduced to tridiagonal form) are desired.
333 *
334 *     I am free to use Q as a valuable working space until Loop 150.
335 *
336             CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
337      $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
338      $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
339      $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
340      $                   IWORK( IPERM ), IWORK( IGIVPT ),
341      $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
342      $                   Q( 1, SUBMAT ), RWORK( IWREM ),
343      $                   IWORK( SUBPBS+1 ), INFO )
344             IF( INFO.GT.0 ) THEN
345                INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
346                RETURN
347             END IF
348             IWORK( I / 2+1 ) = IWORK( I+2 )
349    90    CONTINUE
350          SUBPBS = SUBPBS / 2
351          CURLVL = CURLVL + 1
352          GO TO 80
353       END IF
354 *
355 *     end while
356 *
357 *     Re-merge the eigenvalues/vectors which were deflated at the final
358 *     merge step.
359 *
360       DO 100 I = 1, N
361          J = IWORK( INDXQ+I )
362          RWORK( I ) = D( J )
363          CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
364   100 CONTINUE
365       CALL DCOPY( N, RWORK, 1, D, 1 )
366 *
367       RETURN
368 *
369 *     End of ZLAED0
370 *
371       END