001:       SUBROUTINE DSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
002:      $                   LWORK, IWORK, LIWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          JOBZ, UPLO
011:       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
022: *  a real symmetric band matrix A. If eigenvectors are desired, it uses
023: *  a divide and conquer algorithm.
024: *
025: *  The divide and conquer algorithm makes very mild assumptions about
026: *  floating point arithmetic. It will work on machines with a guard
027: *  digit in add/subtract, or on those binary machines without guard
028: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
029: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
030: *  without guard digits, but we know of none.
031: *
032: *  Arguments
033: *  =========
034: *
035: *  JOBZ    (input) CHARACTER*1
036: *          = 'N':  Compute eigenvalues only;
037: *          = 'V':  Compute eigenvalues and eigenvectors.
038: *
039: *  UPLO    (input) CHARACTER*1
040: *          = 'U':  Upper triangle of A is stored;
041: *          = 'L':  Lower triangle of A is stored.
042: *
043: *  N       (input) INTEGER
044: *          The order of the matrix A.  N >= 0.
045: *
046: *  KD      (input) INTEGER
047: *          The number of superdiagonals of the matrix A if UPLO = 'U',
048: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
049: *
050: *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
051: *          On entry, the upper or lower triangle of the symmetric band
052: *          matrix A, stored in the first KD+1 rows of the array.  The
053: *          j-th column of A is stored in the j-th column of the array AB
054: *          as follows:
055: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
056: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
057: *
058: *          On exit, AB is overwritten by values generated during the
059: *          reduction to tridiagonal form.  If UPLO = 'U', the first
060: *          superdiagonal and the diagonal of the tridiagonal matrix T
061: *          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
062: *          the diagonal and first subdiagonal of T are returned in the
063: *          first two rows of AB.
064: *
065: *  LDAB    (input) INTEGER
066: *          The leading dimension of the array AB.  LDAB >= KD + 1.
067: *
068: *  W       (output) DOUBLE PRECISION array, dimension (N)
069: *          If INFO = 0, the eigenvalues in ascending order.
070: *
071: *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, N)
072: *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
073: *          eigenvectors of the matrix A, with the i-th column of Z
074: *          holding the eigenvector associated with W(i).
075: *          If JOBZ = 'N', then Z is not referenced.
076: *
077: *  LDZ     (input) INTEGER
078: *          The leading dimension of the array Z.  LDZ >= 1, and if
079: *          JOBZ = 'V', LDZ >= max(1,N).
080: *
081: *  WORK    (workspace/output) DOUBLE PRECISION array,
082: *                                         dimension (LWORK)
083: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
084: *
085: *  LWORK   (input) INTEGER
086: *          The dimension of the array WORK.
087: *          IF N <= 1,                LWORK must be at least 1.
088: *          If JOBZ  = 'N' and N > 2, LWORK must be at least 2*N.
089: *          If JOBZ  = 'V' and N > 2, LWORK must be at least
090: *                         ( 1 + 5*N + 2*N**2 ).
091: *
092: *          If LWORK = -1, then a workspace query is assumed; the routine
093: *          only calculates the optimal sizes of the WORK and IWORK
094: *          arrays, returns these values as the first entries of the WORK
095: *          and IWORK arrays, and no error message related to LWORK or
096: *          LIWORK is issued by XERBLA.
097: *
098: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
099: *          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
270: