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