LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssbevd.f
Go to the documentation of this file.
1 *> \brief <b> SSBEVD 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 SSBEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
22 * LWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SSBEVD computes all the eigenvalues and, optionally, eigenvectors of
40 *> a real symmetric band matrix A. If eigenvectors are desired, it uses
41 *> 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] KD
75 *> \verbatim
76 *> KD is INTEGER
77 *> The number of superdiagonals of the matrix A if UPLO = 'U',
78 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
79 *> \endverbatim
80 *>
81 *> \param[in,out] AB
82 *> \verbatim
83 *> AB is REAL array, dimension (LDAB, N)
84 *> On entry, the upper or lower triangle of the symmetric band
85 *> matrix A, stored in the first KD+1 rows of the array. The
86 *> j-th column of A is stored in the j-th column of the array AB
87 *> as follows:
88 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
89 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
90 *>
91 *> On exit, AB is overwritten by values generated during the
92 *> reduction to tridiagonal form. If UPLO = 'U', the first
93 *> superdiagonal and the diagonal of the tridiagonal matrix T
94 *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
95 *> the diagonal and first subdiagonal of T are returned in the
96 *> first two rows of AB.
97 *> \endverbatim
98 *>
99 *> \param[in] LDAB
100 *> \verbatim
101 *> LDAB is INTEGER
102 *> The leading dimension of the array AB. LDAB >= KD + 1.
103 *> \endverbatim
104 *>
105 *> \param[out] W
106 *> \verbatim
107 *> W is REAL array, dimension (N)
108 *> If INFO = 0, the eigenvalues in ascending order.
109 *> \endverbatim
110 *>
111 *> \param[out] Z
112 *> \verbatim
113 *> Z is REAL array, dimension (LDZ, N)
114 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
115 *> eigenvectors of the matrix A, with the i-th column of Z
116 *> holding the eigenvector associated with W(i).
117 *> If JOBZ = 'N', then Z is not referenced.
118 *> \endverbatim
119 *>
120 *> \param[in] LDZ
121 *> \verbatim
122 *> LDZ is INTEGER
123 *> The leading dimension of the array Z. LDZ >= 1, and if
124 *> JOBZ = 'V', LDZ >= max(1,N).
125 *> \endverbatim
126 *>
127 *> \param[out] WORK
128 *> \verbatim
129 *> WORK is REAL array,
130 *> dimension (LWORK)
131 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
132 *> \endverbatim
133 *>
134 *> \param[in] LWORK
135 *> \verbatim
136 *> LWORK is INTEGER
137 *> The dimension of the array WORK.
138 *> IF N <= 1, LWORK must be at least 1.
139 *> If JOBZ = 'N' and N > 2, LWORK must be at least 2*N.
140 *> If JOBZ = 'V' and N > 2, LWORK must be at least
141 *> ( 1 + 5*N + 2*N**2 ).
142 *>
143 *> If LWORK = -1, then a workspace query is assumed; the routine
144 *> only calculates the optimal sizes of the WORK and IWORK
145 *> arrays, returns these values as the first entries of the WORK
146 *> and IWORK arrays, and no error message related to LWORK or
147 *> LIWORK is issued by XERBLA.
148 *> \endverbatim
149 *>
150 *> \param[out] IWORK
151 *> \verbatim
152 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
153 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
154 *> \endverbatim
155 *>
156 *> \param[in] LIWORK
157 *> \verbatim
158 *> LIWORK is INTEGER
159 *> The dimension of the array IWORK.
160 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
161 *> If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
162 *>
163 *> If LIWORK = -1, then a workspace query is assumed; the
164 *> routine only calculates the optimal sizes of the WORK and
165 *> IWORK arrays, returns these values as the first entries of
166 *> the WORK and IWORK arrays, and no error message related to
167 *> LWORK or LIWORK is issued by XERBLA.
168 *> \endverbatim
169 *>
170 *> \param[out] INFO
171 *> \verbatim
172 *> INFO is INTEGER
173 *> = 0: successful exit
174 *> < 0: if INFO = -i, the i-th argument had an illegal value
175 *> > 0: if INFO = i, the algorithm failed to converge; i
176 *> off-diagonal elements of an intermediate tridiagonal
177 *> form did not converge to zero.
178 *> \endverbatim
179 *
180 * Authors:
181 * ========
182 *
183 *> \author Univ. of Tennessee
184 *> \author Univ. of California Berkeley
185 *> \author Univ. of Colorado Denver
186 *> \author NAG Ltd.
187 *
188 *> \ingroup realOTHEReigen
189 *
190 * =====================================================================
191  SUBROUTINE ssbevd( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
192  $ LWORK, IWORK, LIWORK, INFO )
193 *
194 * -- LAPACK driver routine --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 *
198 * .. Scalar Arguments ..
199  CHARACTER JOBZ, UPLO
200  INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
201 * ..
202 * .. Array Arguments ..
203  INTEGER IWORK( * )
204  REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
205 * ..
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210  REAL ZERO, ONE
211  parameter( zero = 0.0e+0, one = 1.0e+0 )
212 * ..
213 * .. Local Scalars ..
214  LOGICAL LOWER, LQUERY, WANTZ
215  INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
216  $ llwrk2, lwmin
217  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218  $ smlnum
219 * ..
220 * .. External Functions ..
221  LOGICAL LSAME
222  REAL SLAMCH, SLANSB
223  EXTERNAL lsame, slamch, slansb
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL sgemm, slacpy, slascl, ssbtrd, sscal, sstedc,
227  $ ssterf, xerbla
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC sqrt
231 * ..
232 * .. Executable Statements ..
233 *
234 * Test the input parameters.
235 *
236  wantz = lsame( jobz, 'V' )
237  lower = lsame( uplo, 'L' )
238  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
239 *
240  info = 0
241  IF( n.LE.1 ) THEN
242  liwmin = 1
243  lwmin = 1
244  ELSE
245  IF( wantz ) THEN
246  liwmin = 3 + 5*n
247  lwmin = 1 + 5*n + 2*n**2
248  ELSE
249  liwmin = 1
250  lwmin = 2*n
251  END IF
252  END IF
253  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
254  info = -1
255  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
256  info = -2
257  ELSE IF( n.LT.0 ) THEN
258  info = -3
259  ELSE IF( kd.LT.0 ) THEN
260  info = -4
261  ELSE IF( ldab.LT.kd+1 ) THEN
262  info = -6
263  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
264  info = -9
265  END IF
266 *
267  IF( info.EQ.0 ) THEN
268  work( 1 ) = lwmin
269  iwork( 1 ) = liwmin
270 *
271  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272  info = -11
273  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
274  info = -13
275  END IF
276  END IF
277 *
278  IF( info.NE.0 ) THEN
279  CALL xerbla( 'SSBEVD', -info )
280  RETURN
281  ELSE IF( lquery ) THEN
282  RETURN
283  END IF
284 *
285 * Quick return if possible
286 *
287  IF( n.EQ.0 )
288  $ RETURN
289 *
290  IF( n.EQ.1 ) THEN
291  w( 1 ) = ab( 1, 1 )
292  IF( wantz )
293  $ z( 1, 1 ) = one
294  RETURN
295  END IF
296 *
297 * Get machine constants.
298 *
299  safmin = slamch( 'Safe minimum' )
300  eps = slamch( 'Precision' )
301  smlnum = safmin / eps
302  bignum = one / smlnum
303  rmin = sqrt( smlnum )
304  rmax = sqrt( bignum )
305 *
306 * Scale matrix to allowable range, if necessary.
307 *
308  anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
309  iscale = 0
310  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
311  iscale = 1
312  sigma = rmin / anrm
313  ELSE IF( anrm.GT.rmax ) THEN
314  iscale = 1
315  sigma = rmax / anrm
316  END IF
317  IF( iscale.EQ.1 ) THEN
318  IF( lower ) THEN
319  CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
320  ELSE
321  CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
322  END IF
323  END IF
324 *
325 * Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
326 *
327  inde = 1
328  indwrk = inde + n
329  indwk2 = indwrk + n*n
330  llwrk2 = lwork - indwk2 + 1
331  CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
332  $ work( indwrk ), iinfo )
333 *
334 * For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC.
335 *
336  IF( .NOT.wantz ) THEN
337  CALL ssterf( n, w, work( inde ), info )
338  ELSE
339  CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
340  $ work( indwk2 ), llwrk2, iwork, liwork, info )
341  CALL sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
342  $ zero, work( indwk2 ), n )
343  CALL slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
344  END IF
345 *
346 * If matrix was scaled, then rescale eigenvalues appropriately.
347 *
348  IF( iscale.EQ.1 )
349  $ CALL sscal( n, one / sigma, w, 1 )
350 *
351  work( 1 ) = lwmin
352  iwork( 1 ) = liwmin
353  RETURN
354 *
355 * End of SSBEVD
356 *
357  END
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 ssbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
SSBTRD
Definition: ssbtrd.f:163
subroutine ssbevd(JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: ssbevd.f:193
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187