1*> \brief <b> DSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSPEVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspevd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspevd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspevd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
22*                          IWORK, LIWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBZ, UPLO
26*       INTEGER            INFO, LDZ, LIWORK, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       DOUBLE PRECISION   AP( * ), W( * ), WORK( * ), Z( LDZ, * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> DSPEVD computes all the eigenvalues and, optionally, eigenvectors
40*> of a real symmetric matrix A in packed storage. If eigenvectors are
41*> desired, it uses a divide and conquer algorithm.
42*>
43*> The divide and conquer algorithm makes very mild assumptions about
44*> floating point arithmetic. It will work on machines with a guard
45*> digit in add/subtract, or on those binary machines without guard
46*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
47*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
48*> without guard digits, but we know of none.
49*> \endverbatim
50*
51*  Arguments:
52*  ==========
53*
54*> \param[in] JOBZ
55*> \verbatim
56*>          JOBZ is CHARACTER*1
57*>          = 'N':  Compute eigenvalues only;
58*>          = 'V':  Compute eigenvalues and eigenvectors.
59*> \endverbatim
60*>
61*> \param[in] UPLO
62*> \verbatim
63*>          UPLO is CHARACTER*1
64*>          = 'U':  Upper triangle of A is stored;
65*>          = 'L':  Lower triangle of A is stored.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*>          N is INTEGER
71*>          The order of the matrix A.  N >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] AP
75*> \verbatim
76*>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
77*>          On entry, the upper or lower triangle of the symmetric matrix
78*>          A, packed columnwise in a linear array.  The j-th column of A
79*>          is stored in the array AP as follows:
80*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
81*>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
82*>
83*>          On exit, AP is overwritten by values generated during the
84*>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
85*>          and first superdiagonal of the tridiagonal matrix T overwrite
86*>          the corresponding elements of A, and if UPLO = 'L', the
87*>          diagonal and first subdiagonal of T overwrite the
88*>          corresponding elements of A.
89*> \endverbatim
90*>
91*> \param[out] W
92*> \verbatim
93*>          W is DOUBLE PRECISION array, dimension (N)
94*>          If INFO = 0, the eigenvalues in ascending order.
95*> \endverbatim
96*>
97*> \param[out] Z
98*> \verbatim
99*>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
100*>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
101*>          eigenvectors of the matrix A, with the i-th column of Z
102*>          holding the eigenvector associated with W(i).
103*>          If JOBZ = 'N', then Z is not referenced.
104*> \endverbatim
105*>
106*> \param[in] LDZ
107*> \verbatim
108*>          LDZ is INTEGER
109*>          The leading dimension of the array Z.  LDZ >= 1, and if
110*>          JOBZ = 'V', LDZ >= max(1,N).
111*> \endverbatim
112*>
113*> \param[out] WORK
114*> \verbatim
115*>          WORK is DOUBLE PRECISION array,
116*>                                         dimension (LWORK)
117*>          On exit, if INFO = 0, WORK(1) returns the required LWORK.
118*> \endverbatim
119*>
120*> \param[in] LWORK
121*> \verbatim
122*>          LWORK is INTEGER
123*>          The dimension of the array WORK.
124*>          If N <= 1,               LWORK must be at least 1.
125*>          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
126*>          If JOBZ = 'V' and N > 1, LWORK must be at least
127*>                                                 1 + 6*N + N**2.
128*>
129*>          If LWORK = -1, then a workspace query is assumed; the routine
130*>          only calculates the required sizes of the WORK and IWORK
131*>          arrays, returns these values as the first entries of the WORK
132*>          and IWORK arrays, and no error message related to LWORK or
133*>          LIWORK is issued by XERBLA.
134*> \endverbatim
135*>
136*> \param[out] IWORK
137*> \verbatim
138*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
139*>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
140*> \endverbatim
141*>
142*> \param[in] LIWORK
143*> \verbatim
144*>          LIWORK is INTEGER
145*>          The dimension of the array IWORK.
146*>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
147*>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
148*>
149*>          If LIWORK = -1, then a workspace query is assumed; the
150*>          routine only calculates the required sizes of the WORK and
151*>          IWORK arrays, returns these values as the first entries of
152*>          the WORK and IWORK arrays, and no error message related to
153*>          LWORK or LIWORK is issued by XERBLA.
154*> \endverbatim
155*>
156*> \param[out] INFO
157*> \verbatim
158*>          INFO is INTEGER
159*>          = 0:  successful exit
160*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
161*>          > 0:  if INFO = i, the algorithm failed to converge; i
162*>                off-diagonal elements of an intermediate tridiagonal
163*>                form did not converge to zero.
164*> \endverbatim
165*
166*  Authors:
167*  ========
168*
169*> \author Univ. of Tennessee
170*> \author Univ. of California Berkeley
171*> \author Univ. of Colorado Denver
172*> \author NAG Ltd.
173*
174*> \date November 2011
175*
176*> \ingroup doubleOTHEReigen
177*
178*  =====================================================================
179      SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
180     $                   IWORK, LIWORK, INFO )
181*
182*  -- LAPACK driver routine (version 3.4.0) --
183*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
184*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185*     November 2011
186*
187*     .. Scalar Arguments ..
188      CHARACTER          JOBZ, UPLO
189      INTEGER            INFO, LDZ, LIWORK, LWORK, N
190*     ..
191*     .. Array Arguments ..
192      INTEGER            IWORK( * )
193      DOUBLE PRECISION   AP( * ), W( * ), WORK( * ), Z( LDZ, * )
194*     ..
195*
196*  =====================================================================
197*
198*     .. Parameters ..
199      DOUBLE PRECISION   ZERO, ONE
200      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
201*     ..
202*     .. Local Scalars ..
203      LOGICAL            LQUERY, WANTZ
204      INTEGER            IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
205     $                   LLWORK, LWMIN
206      DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
207     $                   SMLNUM
208*     ..
209*     .. External Functions ..
210      LOGICAL            LSAME
211      DOUBLE PRECISION   DLAMCH, DLANSP
212      EXTERNAL           LSAME, DLAMCH, DLANSP
213*     ..
214*     .. External Subroutines ..
215      EXTERNAL           DOPMTR, DSCAL, DSPTRD, DSTEDC, DSTERF, XERBLA
216*     ..
217*     .. Intrinsic Functions ..
218      INTRINSIC          SQRT
219*     ..
220*     .. Executable Statements ..
221*
222*     Test the input parameters.
223*
224      WANTZ = LSAME( JOBZ, 'V' )
225      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
226*
227      INFO = 0
228      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
229         INFO = -1
230      ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
231     $          THEN
232         INFO = -2
233      ELSE IF( N.LT.0 ) THEN
234         INFO = -3
235      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
236         INFO = -7
237      END IF
238*
239      IF( INFO.EQ.0 ) THEN
240         IF( N.LE.1 ) THEN
241            LIWMIN = 1
242            LWMIN = 1
243         ELSE
244            IF( WANTZ ) THEN
245               LIWMIN = 3 + 5*N
246               LWMIN = 1 + 6*N + N**2
247            ELSE
248               LIWMIN = 1
249               LWMIN = 2*N
250            END IF
251         END IF
252         IWORK( 1 ) = LIWMIN
253         WORK( 1 ) = LWMIN
254*
255         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
256            INFO = -9
257         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
258            INFO = -11
259         END IF
260      END IF
261*
262      IF( INFO.NE.0 ) THEN
263         CALL XERBLA( 'DSPEVD', -INFO )
264         RETURN
265      ELSE IF( LQUERY ) THEN
266         RETURN
267      END IF
268*
269*     Quick return if possible
270*
271      IF( N.EQ.0 )
272     $   RETURN
273*
274      IF( N.EQ.1 ) THEN
275         W( 1 ) = AP( 1 )
276         IF( WANTZ )
277     $      Z( 1, 1 ) = ONE
278         RETURN
279      END IF
280*
281*     Get machine constants.
282*
283      SAFMIN = DLAMCH( 'Safe minimum' )
284      EPS = DLAMCH( 'Precision' )
285      SMLNUM = SAFMIN / EPS
286      BIGNUM = ONE / SMLNUM
287      RMIN = SQRT( SMLNUM )
288      RMAX = SQRT( BIGNUM )
289*
290*     Scale matrix to allowable range, if necessary.
291*
292      ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
293      ISCALE = 0
294      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
295         ISCALE = 1
296         SIGMA = RMIN / ANRM
297      ELSE IF( ANRM.GT.RMAX ) THEN
298         ISCALE = 1
299         SIGMA = RMAX / ANRM
300      END IF
301      IF( ISCALE.EQ.1 ) THEN
302         CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
303      END IF
304*
305*     Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
306*
307      INDE = 1
308      INDTAU = INDE + N
309      CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
310*
311*     For eigenvalues only, call DSTERF.  For eigenvectors, first call
312*     DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
313*     tridiagonal matrix, then call DOPMTR to multiply it by the
314*     Householder transformations represented in AP.
315*
316      IF( .NOT.WANTZ ) THEN
317         CALL DSTERF( N, W, WORK( INDE ), INFO )
318      ELSE
319         INDWRK = INDTAU + N
320         LLWORK = LWORK - INDWRK + 1
321         CALL DSTEDC( 'I', N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
322     $                LLWORK, IWORK, LIWORK, INFO )
323         CALL DOPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ,
324     $                WORK( INDWRK ), IINFO )
325      END IF
326*
327*     If matrix was scaled, then rescale eigenvalues appropriately.
328*
329      IF( ISCALE.EQ.1 )
330     $   CALL DSCAL( N, ONE / SIGMA, W, 1 )
331*
332      WORK( 1 ) = LWMIN
333      IWORK( 1 ) = LIWMIN
334      RETURN
335*
336*     End of DSPEVD
337*
338      END
339