File:
[local] /
rpl /
lapack /
lapack /
dsbev.f
Revision
1.7:
download - view:
text,
annotated -
select for diffs -
revision graph
Tue Dec 21 13:53:37 2010 UTC (13 years, 6 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_3,
rpl-4_1_2,
rpl-4_1_1,
rpl-4_1_0,
rpl-4_0_24,
rpl-4_0_22,
rpl-4_0_21,
rpl-4_0_20,
rpl-4_0,
HEAD
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DSBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
2: $ 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, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DSBEV computes all the eigenvalues and, optionally, eigenvectors of
21: * a real symmetric band matrix A.
22: *
23: * Arguments
24: * =========
25: *
26: * JOBZ (input) CHARACTER*1
27: * = 'N': Compute eigenvalues only;
28: * = 'V': Compute eigenvalues and eigenvectors.
29: *
30: * UPLO (input) CHARACTER*1
31: * = 'U': Upper triangle of A is stored;
32: * = 'L': Lower triangle of A is stored.
33: *
34: * N (input) INTEGER
35: * The order of the matrix A. N >= 0.
36: *
37: * KD (input) INTEGER
38: * The number of superdiagonals of the matrix A if UPLO = 'U',
39: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
40: *
41: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
42: * On entry, the upper or lower triangle of the symmetric band
43: * matrix A, stored in the first KD+1 rows of the array. The
44: * j-th column of A is stored in the j-th column of the array AB
45: * as follows:
46: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
47: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
48: *
49: * On exit, AB is overwritten by values generated during the
50: * reduction to tridiagonal form. If UPLO = 'U', the first
51: * superdiagonal and the diagonal of the tridiagonal matrix T
52: * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
53: * the diagonal and first subdiagonal of T are returned in the
54: * first two rows of AB.
55: *
56: * LDAB (input) INTEGER
57: * The leading dimension of the array AB. LDAB >= KD + 1.
58: *
59: * W (output) DOUBLE PRECISION array, dimension (N)
60: * If INFO = 0, the eigenvalues in ascending order.
61: *
62: * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
63: * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
64: * eigenvectors of the matrix A, with the i-th column of Z
65: * holding the eigenvector associated with W(i).
66: * If JOBZ = 'N', then Z is not referenced.
67: *
68: * LDZ (input) INTEGER
69: * The leading dimension of the array Z. LDZ >= 1, and if
70: * JOBZ = 'V', LDZ >= max(1,N).
71: *
72: * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
73: *
74: * INFO (output) INTEGER
75: * = 0: successful exit
76: * < 0: if INFO = -i, the i-th argument had an illegal value
77: * > 0: if INFO = i, the algorithm failed to converge; i
78: * off-diagonal elements of an intermediate tridiagonal
79: * form did not converge to zero.
80: *
81: * =====================================================================
82: *
83: * .. Parameters ..
84: DOUBLE PRECISION ZERO, ONE
85: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
86: * ..
87: * .. Local Scalars ..
88: LOGICAL LOWER, WANTZ
89: INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE
90: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
91: $ SMLNUM
92: * ..
93: * .. External Functions ..
94: LOGICAL LSAME
95: DOUBLE PRECISION DLAMCH, DLANSB
96: EXTERNAL LSAME, DLAMCH, DLANSB
97: * ..
98: * .. External Subroutines ..
99: EXTERNAL DLASCL, DSBTRD, DSCAL, DSTEQR, DSTERF, XERBLA
100: * ..
101: * .. Intrinsic Functions ..
102: INTRINSIC SQRT
103: * ..
104: * .. Executable Statements ..
105: *
106: * Test the input parameters.
107: *
108: WANTZ = LSAME( JOBZ, 'V' )
109: LOWER = LSAME( UPLO, 'L' )
110: *
111: INFO = 0
112: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
113: INFO = -1
114: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
115: INFO = -2
116: ELSE IF( N.LT.0 ) THEN
117: INFO = -3
118: ELSE IF( KD.LT.0 ) THEN
119: INFO = -4
120: ELSE IF( LDAB.LT.KD+1 ) THEN
121: INFO = -6
122: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
123: INFO = -9
124: END IF
125: *
126: IF( INFO.NE.0 ) THEN
127: CALL XERBLA( 'DSBEV ', -INFO )
128: RETURN
129: END IF
130: *
131: * Quick return if possible
132: *
133: IF( N.EQ.0 )
134: $ RETURN
135: *
136: IF( N.EQ.1 ) THEN
137: IF( LOWER ) THEN
138: W( 1 ) = AB( 1, 1 )
139: ELSE
140: W( 1 ) = AB( KD+1, 1 )
141: END IF
142: IF( WANTZ )
143: $ Z( 1, 1 ) = ONE
144: RETURN
145: END IF
146: *
147: * Get machine constants.
148: *
149: SAFMIN = DLAMCH( 'Safe minimum' )
150: EPS = DLAMCH( 'Precision' )
151: SMLNUM = SAFMIN / EPS
152: BIGNUM = ONE / SMLNUM
153: RMIN = SQRT( SMLNUM )
154: RMAX = SQRT( BIGNUM )
155: *
156: * Scale matrix to allowable range, if necessary.
157: *
158: ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
159: ISCALE = 0
160: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
161: ISCALE = 1
162: SIGMA = RMIN / ANRM
163: ELSE IF( ANRM.GT.RMAX ) THEN
164: ISCALE = 1
165: SIGMA = RMAX / ANRM
166: END IF
167: IF( ISCALE.EQ.1 ) THEN
168: IF( LOWER ) THEN
169: CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
170: ELSE
171: CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
172: END IF
173: END IF
174: *
175: * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
176: *
177: INDE = 1
178: INDWRK = INDE + N
179: CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
180: $ WORK( INDWRK ), IINFO )
181: *
182: * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEQR.
183: *
184: IF( .NOT.WANTZ ) THEN
185: CALL DSTERF( N, W, WORK( INDE ), INFO )
186: ELSE
187: CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
188: $ INFO )
189: END IF
190: *
191: * If matrix was scaled, then rescale eigenvalues appropriately.
192: *
193: IF( ISCALE.EQ.1 ) THEN
194: IF( INFO.EQ.0 ) THEN
195: IMAX = N
196: ELSE
197: IMAX = INFO - 1
198: END IF
199: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
200: END IF
201: *
202: RETURN
203: *
204: * End of DSBEV
205: *
206: END
CVSweb interface <joel.bertrand@systella.fr>