LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sspevd.f
Go to the documentation of this file.
1 *> \brief <b> SSPEVD 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 SSPEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sspevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sspevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sspevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSPEVD( 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 * REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SSPEVD 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 REAL 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 REAL array, dimension (N)
94 *> If INFO = 0, the eigenvalues in ascending order.
95 *> \endverbatim
96 *>
97 *> \param[out] Z
98 *> \verbatim
99 *> Z is REAL 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 REAL array, dimension (MAX(1,LWORK))
116 *> On exit, if INFO = 0, WORK(1) returns the required LWORK.
117 *> \endverbatim
118 *>
119 *> \param[in] LWORK
120 *> \verbatim
121 *> LWORK is INTEGER
122 *> The dimension of the array WORK.
123 *> If N <= 1, LWORK must be at least 1.
124 *> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
125 *> If JOBZ = 'V' and N > 1, LWORK must be at least
126 *> 1 + 6*N + N**2.
127 *>
128 *> If LWORK = -1, then a workspace query is assumed; the routine
129 *> only calculates the required sizes of the WORK and IWORK
130 *> arrays, returns these values as the first entries of the WORK
131 *> and IWORK arrays, and no error message related to LWORK or
132 *> LIWORK is issued by XERBLA.
133 *> \endverbatim
134 *>
135 *> \param[out] IWORK
136 *> \verbatim
137 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
138 *> On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
139 *> \endverbatim
140 *>
141 *> \param[in] LIWORK
142 *> \verbatim
143 *> LIWORK is INTEGER
144 *> The dimension of the array IWORK.
145 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
146 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
147 *>
148 *> If LIWORK = -1, then a workspace query is assumed; the
149 *> routine only calculates the required sizes of the WORK and
150 *> IWORK arrays, returns these values as the first entries of
151 *> the WORK and IWORK arrays, and no error message related to
152 *> LWORK or LIWORK is issued by XERBLA.
153 *> \endverbatim
154 *>
155 *> \param[out] INFO
156 *> \verbatim
157 *> INFO is INTEGER
158 *> = 0: successful exit
159 *> < 0: if INFO = -i, the i-th argument had an illegal value.
160 *> > 0: if INFO = i, the algorithm failed to converge; i
161 *> off-diagonal elements of an intermediate tridiagonal
162 *> form did not converge to zero.
163 *> \endverbatim
164 *
165 * Authors:
166 * ========
167 *
168 *> \author Univ. of Tennessee
169 *> \author Univ. of California Berkeley
170 *> \author Univ. of Colorado Denver
171 *> \author NAG Ltd.
172 *
173 *> \ingroup realOTHEReigen
174 *
175 * =====================================================================
176  SUBROUTINE sspevd( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
177  $ IWORK, LIWORK, INFO )
178 *
179 * -- LAPACK driver routine --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 *
183 * .. Scalar Arguments ..
184  CHARACTER JOBZ, UPLO
185  INTEGER INFO, LDZ, LIWORK, LWORK, N
186 * ..
187 * .. Array Arguments ..
188  INTEGER IWORK( * )
189  REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  REAL ZERO, ONE
196  parameter( zero = 0.0e+0, one = 1.0e+0 )
197 * ..
198 * .. Local Scalars ..
199  LOGICAL LQUERY, WANTZ
200  INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
201  $ llwork, lwmin
202  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
203  $ smlnum
204 * ..
205 * .. External Functions ..
206  LOGICAL LSAME
207  REAL SLAMCH, SLANSP
208  EXTERNAL lsame, slamch, slansp
209 * ..
210 * .. External Subroutines ..
211  EXTERNAL sopmtr, sscal, ssptrd, sstedc, ssterf, xerbla
212 * ..
213 * .. Intrinsic Functions ..
214  INTRINSIC sqrt
215 * ..
216 * .. Executable Statements ..
217 *
218 * Test the input parameters.
219 *
220  wantz = lsame( jobz, 'V' )
221  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
222 *
223  info = 0
224  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
225  info = -1
226  ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
227  $ THEN
228  info = -2
229  ELSE IF( n.LT.0 ) THEN
230  info = -3
231  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
232  info = -7
233  END IF
234 *
235  IF( info.EQ.0 ) THEN
236  IF( n.LE.1 ) THEN
237  liwmin = 1
238  lwmin = 1
239  ELSE
240  IF( wantz ) THEN
241  liwmin = 3 + 5*n
242  lwmin = 1 + 6*n + n**2
243  ELSE
244  liwmin = 1
245  lwmin = 2*n
246  END IF
247  END IF
248  iwork( 1 ) = liwmin
249  work( 1 ) = lwmin
250 *
251  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
252  info = -9
253  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
254  info = -11
255  END IF
256  END IF
257 *
258  IF( info.NE.0 ) THEN
259  CALL xerbla( 'SSPEVD', -info )
260  RETURN
261  ELSE IF( lquery ) THEN
262  RETURN
263  END IF
264 *
265 * Quick return if possible
266 *
267  IF( n.EQ.0 )
268  $ RETURN
269 *
270  IF( n.EQ.1 ) THEN
271  w( 1 ) = ap( 1 )
272  IF( wantz )
273  $ z( 1, 1 ) = one
274  RETURN
275  END IF
276 *
277 * Get machine constants.
278 *
279  safmin = slamch( 'Safe minimum' )
280  eps = slamch( 'Precision' )
281  smlnum = safmin / eps
282  bignum = one / smlnum
283  rmin = sqrt( smlnum )
284  rmax = sqrt( bignum )
285 *
286 * Scale matrix to allowable range, if necessary.
287 *
288  anrm = slansp( 'M', uplo, n, ap, work )
289  iscale = 0
290  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
291  iscale = 1
292  sigma = rmin / anrm
293  ELSE IF( anrm.GT.rmax ) THEN
294  iscale = 1
295  sigma = rmax / anrm
296  END IF
297  IF( iscale.EQ.1 ) THEN
298  CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
299  END IF
300 *
301 * Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
302 *
303  inde = 1
304  indtau = inde + n
305  CALL ssptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
306 *
307 * For eigenvalues only, call SSTERF. For eigenvectors, first call
308 * SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
309 * tridiagonal matrix, then call SOPMTR to multiply it by the
310 * Householder transformations represented in AP.
311 *
312  IF( .NOT.wantz ) THEN
313  CALL ssterf( n, w, work( inde ), info )
314  ELSE
315  indwrk = indtau + n
316  llwork = lwork - indwrk + 1
317  CALL sstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
318  $ llwork, iwork, liwork, info )
319  CALL sopmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
320  $ work( indwrk ), iinfo )
321  END IF
322 *
323 * If matrix was scaled, then rescale eigenvalues appropriately.
324 *
325  IF( iscale.EQ.1 )
326  $ CALL sscal( n, one / sigma, w, 1 )
327 *
328  work( 1 ) = lwmin
329  iwork( 1 ) = liwmin
330  RETURN
331 *
332 * End of SSPEVD
333 *
334  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine ssptrd(UPLO, N, AP, D, E, TAU, INFO)
SSPTRD
Definition: ssptrd.f:150
subroutine sopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
SOPMTR
Definition: sopmtr.f:150
subroutine sspevd(JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: sspevd.f:178
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79