1: SUBROUTINE DSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
2: $ LWORK, IWORK, LIWORK, INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER JOBZ, UPLO
11: INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: INTEGER IWORK( * )
15: DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
22: * a real symmetric band matrix A. If eigenvectors are desired, it uses
23: * a divide and conquer algorithm.
24: *
25: * The divide and conquer algorithm makes very mild assumptions about
26: * floating point arithmetic. It will work on machines with a guard
27: * digit in add/subtract, or on those binary machines without guard
28: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30: * without guard digits, but we know of none.
31: *
32: * Arguments
33: * =========
34: *
35: * JOBZ (input) CHARACTER*1
36: * = 'N': Compute eigenvalues only;
37: * = 'V': Compute eigenvalues and eigenvectors.
38: *
39: * UPLO (input) CHARACTER*1
40: * = 'U': Upper triangle of A is stored;
41: * = 'L': Lower triangle of A is stored.
42: *
43: * N (input) INTEGER
44: * The order of the matrix A. N >= 0.
45: *
46: * KD (input) INTEGER
47: * The number of superdiagonals of the matrix A if UPLO = 'U',
48: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
49: *
50: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
51: * On entry, the upper or lower triangle of the symmetric band
52: * matrix A, stored in the first KD+1 rows of the array. The
53: * j-th column of A is stored in the j-th column of the array AB
54: * as follows:
55: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
56: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
57: *
58: * On exit, AB is overwritten by values generated during the
59: * reduction to tridiagonal form. If UPLO = 'U', the first
60: * superdiagonal and the diagonal of the tridiagonal matrix T
61: * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
62: * the diagonal and first subdiagonal of T are returned in the
63: * first two rows of AB.
64: *
65: * LDAB (input) INTEGER
66: * The leading dimension of the array AB. LDAB >= KD + 1.
67: *
68: * W (output) DOUBLE PRECISION array, dimension (N)
69: * If INFO = 0, the eigenvalues in ascending order.
70: *
71: * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
72: * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
73: * eigenvectors of the matrix A, with the i-th column of Z
74: * holding the eigenvector associated with W(i).
75: * If JOBZ = 'N', then Z is not referenced.
76: *
77: * LDZ (input) INTEGER
78: * The leading dimension of the array Z. LDZ >= 1, and if
79: * JOBZ = 'V', LDZ >= max(1,N).
80: *
81: * WORK (workspace/output) DOUBLE PRECISION array,
82: * dimension (LWORK)
83: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84: *
85: * LWORK (input) INTEGER
86: * The dimension of the array WORK.
87: * IF N <= 1, LWORK must be at least 1.
88: * If JOBZ = 'N' and N > 2, LWORK must be at least 2*N.
89: * If JOBZ = 'V' and N > 2, LWORK must be at least
90: * ( 1 + 5*N + 2*N**2 ).
91: *
92: * If LWORK = -1, then a workspace query is assumed; the routine
93: * only calculates the optimal sizes of the WORK and IWORK
94: * arrays, returns these values as the first entries of the WORK
95: * and IWORK arrays, and no error message related to LWORK or
96: * LIWORK is issued by XERBLA.
97: *
98: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
99: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
100: *
101: * LIWORK (input) INTEGER
102: * The dimension of the array LIWORK.
103: * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
104: * If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
105: *
106: * If LIWORK = -1, then a workspace query is assumed; the
107: * routine only calculates the optimal sizes of the WORK and
108: * IWORK arrays, returns these values as the first entries of
109: * the WORK and IWORK arrays, and no error message related to
110: * LWORK or LIWORK is issued by XERBLA.
111: *
112: * INFO (output) INTEGER
113: * = 0: successful exit
114: * < 0: if INFO = -i, the i-th argument had an illegal value
115: * > 0: if INFO = i, the algorithm failed to converge; i
116: * off-diagonal elements of an intermediate tridiagonal
117: * form did not converge to zero.
118: *
119: * =====================================================================
120: *
121: * .. Parameters ..
122: DOUBLE PRECISION ZERO, ONE
123: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124: * ..
125: * .. Local Scalars ..
126: LOGICAL LOWER, LQUERY, WANTZ
127: INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
128: $ LLWRK2, LWMIN
129: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
130: $ SMLNUM
131: * ..
132: * .. External Functions ..
133: LOGICAL LSAME
134: DOUBLE PRECISION DLAMCH, DLANSB
135: EXTERNAL LSAME, DLAMCH, DLANSB
136: * ..
137: * .. External Subroutines ..
138: EXTERNAL DGEMM, DLACPY, DLASCL, DSBTRD, DSCAL, DSTEDC,
139: $ DSTERF, XERBLA
140: * ..
141: * .. Intrinsic Functions ..
142: INTRINSIC SQRT
143: * ..
144: * .. Executable Statements ..
145: *
146: * Test the input parameters.
147: *
148: WANTZ = LSAME( JOBZ, 'V' )
149: LOWER = LSAME( UPLO, 'L' )
150: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
151: *
152: INFO = 0
153: IF( N.LE.1 ) THEN
154: LIWMIN = 1
155: LWMIN = 1
156: ELSE
157: IF( WANTZ ) THEN
158: LIWMIN = 3 + 5*N
159: LWMIN = 1 + 5*N + 2*N**2
160: ELSE
161: LIWMIN = 1
162: LWMIN = 2*N
163: END IF
164: END IF
165: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
166: INFO = -1
167: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
168: INFO = -2
169: ELSE IF( N.LT.0 ) THEN
170: INFO = -3
171: ELSE IF( KD.LT.0 ) THEN
172: INFO = -4
173: ELSE IF( LDAB.LT.KD+1 ) THEN
174: INFO = -6
175: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
176: INFO = -9
177: END IF
178: *
179: IF( INFO.EQ.0 ) THEN
180: WORK( 1 ) = LWMIN
181: IWORK( 1 ) = LIWMIN
182: *
183: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
184: INFO = -11
185: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
186: INFO = -13
187: END IF
188: END IF
189: *
190: IF( INFO.NE.0 ) THEN
191: CALL XERBLA( 'DSBEVD', -INFO )
192: RETURN
193: ELSE IF( LQUERY ) THEN
194: RETURN
195: END IF
196: *
197: * Quick return if possible
198: *
199: IF( N.EQ.0 )
200: $ RETURN
201: *
202: IF( N.EQ.1 ) THEN
203: W( 1 ) = AB( 1, 1 )
204: IF( WANTZ )
205: $ Z( 1, 1 ) = ONE
206: RETURN
207: END IF
208: *
209: * Get machine constants.
210: *
211: SAFMIN = DLAMCH( 'Safe minimum' )
212: EPS = DLAMCH( 'Precision' )
213: SMLNUM = SAFMIN / EPS
214: BIGNUM = ONE / SMLNUM
215: RMIN = SQRT( SMLNUM )
216: RMAX = SQRT( BIGNUM )
217: *
218: * Scale matrix to allowable range, if necessary.
219: *
220: ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
221: ISCALE = 0
222: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
223: ISCALE = 1
224: SIGMA = RMIN / ANRM
225: ELSE IF( ANRM.GT.RMAX ) THEN
226: ISCALE = 1
227: SIGMA = RMAX / ANRM
228: END IF
229: IF( ISCALE.EQ.1 ) THEN
230: IF( LOWER ) THEN
231: CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
232: ELSE
233: CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
234: END IF
235: END IF
236: *
237: * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
238: *
239: INDE = 1
240: INDWRK = INDE + N
241: INDWK2 = INDWRK + N*N
242: LLWRK2 = LWORK - INDWK2 + 1
243: CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
244: $ WORK( INDWRK ), IINFO )
245: *
246: * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
247: *
248: IF( .NOT.WANTZ ) THEN
249: CALL DSTERF( N, W, WORK( INDE ), INFO )
250: ELSE
251: CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
252: $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
253: CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
254: $ ZERO, WORK( INDWK2 ), N )
255: CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
256: END IF
257: *
258: * If matrix was scaled, then rescale eigenvalues appropriately.
259: *
260: IF( ISCALE.EQ.1 )
261: $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
262: *
263: WORK( 1 ) = LWMIN
264: IWORK( 1 ) = LIWMIN
265: RETURN
266: *
267: * End of DSBEVD
268: *
269: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsbevd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack