LAPACK  3.6.1 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
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 *> \date November 2011
174 *
175 *> \ingroup realOTHEReigen
176 *
177 * =====================================================================
178  SUBROUTINE sspevd( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
179  \$ iwork, liwork, info )
180 *
181 * -- LAPACK driver routine (version 3.4.0) --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 * November 2011
185 *
186 * .. Scalar Arguments ..
187  CHARACTER JOBZ, UPLO
188  INTEGER INFO, LDZ, LIWORK, LWORK, N
189 * ..
190 * .. Array Arguments ..
191  INTEGER IWORK( * )
192  REAL AP( * ), W( * ), WORK( * ), Z( ldz, * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  REAL ZERO, ONE
199  parameter ( zero = 0.0e+0, one = 1.0e+0 )
200 * ..
201 * .. Local Scalars ..
202  LOGICAL LQUERY, WANTZ
203  INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
204  \$ llwork, lwmin
205  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
206  \$ smlnum
207 * ..
208 * .. External Functions ..
209  LOGICAL LSAME
210  REAL SLAMCH, SLANSP
211  EXTERNAL lsame, slamch, slansp
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL sopmtr, sscal, ssptrd, sstedc, ssterf, xerbla
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC sqrt
218 * ..
219 * .. Executable Statements ..
220 *
221 * Test the input parameters.
222 *
223  wantz = lsame( jobz, 'V' )
224  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
225 *
226  info = 0
227  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
228  info = -1
229  ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
230  \$ THEN
231  info = -2
232  ELSE IF( n.LT.0 ) THEN
233  info = -3
234  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
235  info = -7
236  END IF
237 *
238  IF( info.EQ.0 ) THEN
239  IF( n.LE.1 ) THEN
240  liwmin = 1
241  lwmin = 1
242  ELSE
243  IF( wantz ) THEN
244  liwmin = 3 + 5*n
245  lwmin = 1 + 6*n + n**2
246  ELSE
247  liwmin = 1
248  lwmin = 2*n
249  END IF
250  END IF
251  iwork( 1 ) = liwmin
252  work( 1 ) = lwmin
253 *
254  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
255  info = -9
256  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
257  info = -11
258  END IF
259  END IF
260 *
261  IF( info.NE.0 ) THEN
262  CALL xerbla( 'SSPEVD', -info )
263  RETURN
264  ELSE IF( lquery ) THEN
265  RETURN
266  END IF
267 *
268 * Quick return if possible
269 *
270  IF( n.EQ.0 )
271  \$ RETURN
272 *
273  IF( n.EQ.1 ) THEN
274  w( 1 ) = ap( 1 )
275  IF( wantz )
276  \$ z( 1, 1 ) = one
277  RETURN
278  END IF
279 *
280 * Get machine constants.
281 *
282  safmin = slamch( 'Safe minimum' )
283  eps = slamch( 'Precision' )
284  smlnum = safmin / eps
285  bignum = one / smlnum
286  rmin = sqrt( smlnum )
287  rmax = sqrt( bignum )
288 *
289 * Scale matrix to allowable range, if necessary.
290 *
291  anrm = slansp( 'M', uplo, n, ap, work )
292  iscale = 0
293  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
294  iscale = 1
295  sigma = rmin / anrm
296  ELSE IF( anrm.GT.rmax ) THEN
297  iscale = 1
298  sigma = rmax / anrm
299  END IF
300  IF( iscale.EQ.1 ) THEN
301  CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
302  END IF
303 *
304 * Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
305 *
306  inde = 1
307  indtau = inde + n
308  CALL ssptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
309 *
310 * For eigenvalues only, call SSTERF. For eigenvectors, first call
311 * SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
312 * tridiagonal matrix, then call SOPMTR to multiply it by the
313 * Householder transformations represented in AP.
314 *
315  IF( .NOT.wantz ) THEN
316  CALL ssterf( n, w, work( inde ), info )
317  ELSE
318  indwrk = indtau + n
319  llwork = lwork - indwrk + 1
320  CALL sstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
321  \$ llwork, iwork, liwork, info )
322  CALL sopmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
323  \$ work( indwrk ), iinfo )
324  END IF
325 *
326 * If matrix was scaled, then rescale eigenvalues appropriately.
327 *
328  IF( iscale.EQ.1 )
329  \$ CALL sscal( n, one / sigma, w, 1 )
330 *
331  work( 1 ) = lwmin
332  iwork( 1 ) = liwmin
333  RETURN
334 *
335 * End of SSPEVD
336 *
337  END
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 matric...
Definition: sspevd.f:180
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
SOPMTR
Definition: sopmtr.f:152
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:190
subroutine ssptrd(UPLO, N, AP, D, E, TAU, INFO)
SSPTRD
Definition: ssptrd.f:152