1*> \brief \b SLAED0 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 SLAED0 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed0.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed0.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed0.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
22*                          WORK, IWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IWORK( * )
29*       REAL               D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
30*      $                   WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SLAED0 computes all eigenvalues and corresponding eigenvectors of a
40*> symmetric tridiagonal matrix using the divide and conquer method.
41*> \endverbatim
42*
43*  Arguments:
44*  ==========
45*
46*> \param[in] ICOMPQ
47*> \verbatim
48*>          ICOMPQ is INTEGER
49*>          = 0:  Compute eigenvalues only.
50*>          = 1:  Compute eigenvectors of original dense symmetric matrix
51*>                also.  On entry, Q contains the orthogonal matrix used
52*>                to reduce the original matrix to tridiagonal form.
53*>          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
54*>                matrix.
55*> \endverbatim
56*>
57*> \param[in] QSIZ
58*> \verbatim
59*>          QSIZ is INTEGER
60*>         The dimension of the orthogonal matrix used to reduce
61*>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*>          N is INTEGER
67*>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] D
71*> \verbatim
72*>          D is REAL array, dimension (N)
73*>         On entry, the main diagonal of the tridiagonal matrix.
74*>         On exit, its eigenvalues.
75*> \endverbatim
76*>
77*> \param[in] E
78*> \verbatim
79*>          E is REAL array, dimension (N-1)
80*>         The off-diagonal elements of the tridiagonal matrix.
81*>         On exit, E has been destroyed.
82*> \endverbatim
83*>
84*> \param[in,out] Q
85*> \verbatim
86*>          Q is REAL array, dimension (LDQ, N)
87*>         On entry, Q must contain an N-by-N orthogonal matrix.
88*>         If ICOMPQ = 0    Q is not referenced.
89*>         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
90*>                          orthogonal matrix used to reduce the full
91*>                          matrix to tridiagonal form corresponding to
92*>                          the subset of the full matrix which is being
93*>                          decomposed at this time.
94*>         If ICOMPQ = 2    On entry, Q will be the identity matrix.
95*>                          On exit, Q contains the eigenvectors of the
96*>                          tridiagonal matrix.
97*> \endverbatim
98*>
99*> \param[in] LDQ
100*> \verbatim
101*>          LDQ is INTEGER
102*>         The leading dimension of the array Q.  If eigenvectors are
103*>         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
104*> \endverbatim
105*>
106*> \param[out] QSTORE
107*> \verbatim
108*>          QSTORE is REAL array, dimension (LDQS, N)
109*>         Referenced only when ICOMPQ = 1.  Used to store parts of
110*>         the eigenvector matrix when the updating matrix multiplies
111*>         take place.
112*> \endverbatim
113*>
114*> \param[in] LDQS
115*> \verbatim
116*>          LDQS is INTEGER
117*>         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
118*>         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
119*> \endverbatim
120*>
121*> \param[out] WORK
122*> \verbatim
123*>          WORK is REAL array,
124*>         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
125*>                     1 + 3*N + 2*N*lg N + 3*N**2
126*>                     ( lg( N ) = smallest integer k
127*>                                 such that 2^k >= N )
128*>         If ICOMPQ = 2, the dimension of WORK must be at least
129*>                     4*N + N**2.
130*> \endverbatim
131*>
132*> \param[out] IWORK
133*> \verbatim
134*>          IWORK is INTEGER array,
135*>         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
136*>                        6 + 6*N + 5*N*lg N.
137*>                        ( lg( N ) = smallest integer k
138*>                                    such that 2^k >= N )
139*>         If ICOMPQ = 2, the dimension of IWORK must be at least
140*>                        3 + 5*N.
141*> \endverbatim
142*>
143*> \param[out] INFO
144*> \verbatim
145*>          INFO is INTEGER
146*>          = 0:  successful exit.
147*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
148*>          > 0:  The algorithm failed to compute an eigenvalue while
149*>                working on the submatrix lying in rows and columns
150*>                INFO/(N+1) through mod(INFO,N+1).
151*> \endverbatim
152*
153*  Authors:
154*  ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \date September 2012
162*
163*> \ingroup auxOTHERcomputational
164*
165*> \par Contributors:
166*  ==================
167*>
168*> Jeff Rutter, Computer Science Division, University of California
169*> at Berkeley, USA
170*
171*  =====================================================================
172      SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
173     $                   WORK, IWORK, INFO )
174*
175*  -- LAPACK computational routine (version 3.4.2) --
176*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
177*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*     September 2012
179*
180*     .. Scalar Arguments ..
181      INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
182*     ..
183*     .. Array Arguments ..
184      INTEGER            IWORK( * )
185      REAL               D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
186     $                   WORK( * )
187*     ..
188*
189*  =====================================================================
190*
191*     .. Parameters ..
192      REAL               ZERO, ONE, TWO
193      PARAMETER          ( ZERO = 0.E0, ONE = 1.E0, TWO = 2.E0 )
194*     ..
195*     .. Local Scalars ..
196      INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
197     $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
198     $                   J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
199     $                   SPM2, SUBMAT, SUBPBS, TLVLS
200      REAL               TEMP
201*     ..
202*     .. External Subroutines ..
203      EXTERNAL           SCOPY, SGEMM, SLACPY, SLAED1, SLAED7, SSTEQR,
204     $                   XERBLA
205*     ..
206*     .. External Functions ..
207      INTEGER            ILAENV
208      EXTERNAL           ILAENV
209*     ..
210*     .. Intrinsic Functions ..
211      INTRINSIC          ABS, INT, LOG, MAX, REAL
212*     ..
213*     .. Executable Statements ..
214*
215*     Test the input parameters.
216*
217      INFO = 0
218*
219      IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
220         INFO = -1
221      ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
222         INFO = -2
223      ELSE IF( N.LT.0 ) THEN
224         INFO = -3
225      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
226         INFO = -7
227      ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
228         INFO = -9
229      END IF
230      IF( INFO.NE.0 ) THEN
231         CALL XERBLA( 'SLAED0', -INFO )
232         RETURN
233      END IF
234*
235*     Quick return if possible
236*
237      IF( N.EQ.0 )
238     $   RETURN
239*
240      SMLSIZ = ILAENV( 9, 'SLAED0', ' ', 0, 0, 0, 0 )
241*
242*     Determine the size and placement of the submatrices, and save in
243*     the leading elements of IWORK.
244*
245      IWORK( 1 ) = N
246      SUBPBS = 1
247      TLVLS = 0
248   10 CONTINUE
249      IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
250         DO 20 J = SUBPBS, 1, -1
251            IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
252            IWORK( 2*J-1 ) = IWORK( J ) / 2
253   20    CONTINUE
254         TLVLS = TLVLS + 1
255         SUBPBS = 2*SUBPBS
256         GO TO 10
257      END IF
258      DO 30 J = 2, SUBPBS
259         IWORK( J ) = IWORK( J ) + IWORK( J-1 )
260   30 CONTINUE
261*
262*     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
263*     using rank-1 modifications (cuts).
264*
265      SPM1 = SUBPBS - 1
266      DO 40 I = 1, SPM1
267         SUBMAT = IWORK( I ) + 1
268         SMM1 = SUBMAT - 1
269         D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
270         D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
271   40 CONTINUE
272*
273      INDXQ = 4*N + 3
274      IF( ICOMPQ.NE.2 ) THEN
275*
276*        Set up workspaces for eigenvalues only/accumulate new vectors
277*        routine
278*
279         TEMP = LOG( REAL( N ) ) / LOG( TWO )
280         LGN = INT( TEMP )
281         IF( 2**LGN.LT.N )
282     $      LGN = LGN + 1
283         IF( 2**LGN.LT.N )
284     $      LGN = LGN + 1
285         IPRMPT = INDXQ + N + 1
286         IPERM = IPRMPT + N*LGN
287         IQPTR = IPERM + N*LGN
288         IGIVPT = IQPTR + N + 2
289         IGIVCL = IGIVPT + N*LGN
290*
291         IGIVNM = 1
292         IQ = IGIVNM + 2*N*LGN
293         IWREM = IQ + N**2 + 1
294*
295*        Initialize pointers
296*
297         DO 50 I = 0, SUBPBS
298            IWORK( IPRMPT+I ) = 1
299            IWORK( IGIVPT+I ) = 1
300   50    CONTINUE
301         IWORK( IQPTR ) = 1
302      END IF
303*
304*     Solve each submatrix eigenproblem at the bottom of the divide and
305*     conquer tree.
306*
307      CURR = 0
308      DO 70 I = 0, SPM1
309         IF( I.EQ.0 ) THEN
310            SUBMAT = 1
311            MATSIZ = IWORK( 1 )
312         ELSE
313            SUBMAT = IWORK( I ) + 1
314            MATSIZ = IWORK( I+1 ) - IWORK( I )
315         END IF
316         IF( ICOMPQ.EQ.2 ) THEN
317            CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
318     $                   Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
319            IF( INFO.NE.0 )
320     $         GO TO 130
321         ELSE
322            CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
323     $                   WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
324     $                   INFO )
325            IF( INFO.NE.0 )
326     $         GO TO 130
327            IF( ICOMPQ.EQ.1 ) THEN
328               CALL SGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
329     $                     Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
330     $                     CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
331     $                     LDQS )
332            END IF
333            IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
334            CURR = CURR + 1
335         END IF
336         K = 1
337         DO 60 J = SUBMAT, IWORK( I+1 )
338            IWORK( INDXQ+J ) = K
339            K = K + 1
340   60    CONTINUE
341   70 CONTINUE
342*
343*     Successively merge eigensystems of adjacent submatrices
344*     into eigensystem for the corresponding larger matrix.
345*
346*     while ( SUBPBS > 1 )
347*
348      CURLVL = 1
349   80 CONTINUE
350      IF( SUBPBS.GT.1 ) THEN
351         SPM2 = SUBPBS - 2
352         DO 90 I = 0, SPM2, 2
353            IF( I.EQ.0 ) THEN
354               SUBMAT = 1
355               MATSIZ = IWORK( 2 )
356               MSD2 = IWORK( 1 )
357               CURPRB = 0
358            ELSE
359               SUBMAT = IWORK( I ) + 1
360               MATSIZ = IWORK( I+2 ) - IWORK( I )
361               MSD2 = MATSIZ / 2
362               CURPRB = CURPRB + 1
363            END IF
364*
365*     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
366*     into an eigensystem of size MATSIZ.
367*     SLAED1 is used only for the full eigensystem of a tridiagonal
368*     matrix.
369*     SLAED7 handles the cases in which eigenvalues only or eigenvalues
370*     and eigenvectors of a full symmetric matrix (which was reduced to
371*     tridiagonal form) are desired.
372*
373            IF( ICOMPQ.EQ.2 ) THEN
374               CALL SLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
375     $                      LDQ, IWORK( INDXQ+SUBMAT ),
376     $                      E( SUBMAT+MSD2-1 ), MSD2, WORK,
377     $                      IWORK( SUBPBS+1 ), INFO )
378            ELSE
379               CALL SLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
380     $                      D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
381     $                      IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
382     $                      MSD2, WORK( IQ ), IWORK( IQPTR ),
383     $                      IWORK( IPRMPT ), IWORK( IPERM ),
384     $                      IWORK( IGIVPT ), IWORK( IGIVCL ),
385     $                      WORK( IGIVNM ), WORK( IWREM ),
386     $                      IWORK( SUBPBS+1 ), INFO )
387            END IF
388            IF( INFO.NE.0 )
389     $         GO TO 130
390            IWORK( I / 2+1 ) = IWORK( I+2 )
391   90    CONTINUE
392         SUBPBS = SUBPBS / 2
393         CURLVL = CURLVL + 1
394         GO TO 80
395      END IF
396*
397*     end while
398*
399*     Re-merge the eigenvalues/vectors which were deflated at the final
400*     merge step.
401*
402      IF( ICOMPQ.EQ.1 ) THEN
403         DO 100 I = 1, N
404            J = IWORK( INDXQ+I )
405            WORK( I ) = D( J )
406            CALL SCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
407  100    CONTINUE
408         CALL SCOPY( N, WORK, 1, D, 1 )
409      ELSE IF( ICOMPQ.EQ.2 ) THEN
410         DO 110 I = 1, N
411            J = IWORK( INDXQ+I )
412            WORK( I ) = D( J )
413            CALL SCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
414  110    CONTINUE
415         CALL SCOPY( N, WORK, 1, D, 1 )
416         CALL SLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
417      ELSE
418         DO 120 I = 1, N
419            J = IWORK( INDXQ+I )
420            WORK( I ) = D( J )
421  120    CONTINUE
422         CALL SCOPY( N, WORK, 1, D, 1 )
423      END IF
424      GO TO 140
425*
426  130 CONTINUE
427      INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
428*
429  140 CONTINUE
430      RETURN
431*
432*     End of SLAED0
433*
434      END
435