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: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOBZ, UPLO
010:       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IWORK( * )
014:       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
021: *  a real symmetric band matrix A. If eigenvectors are desired, it uses
022: *  a divide and conquer algorithm.
023: *
024: *  The divide and conquer algorithm makes very mild assumptions about
025: *  floating point arithmetic. It will work on machines with a guard
026: *  digit in add/subtract, or on those binary machines without guard
027: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
028: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
029: *  without guard digits, but we know of none.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  JOBZ    (input) CHARACTER*1
035: *          = 'N':  Compute eigenvalues only;
036: *          = 'V':  Compute eigenvalues and eigenvectors.
037: *
038: *  UPLO    (input) CHARACTER*1
039: *          = 'U':  Upper triangle of A is stored;
040: *          = 'L':  Lower triangle of A is stored.
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrix A.  N >= 0.
044: *
045: *  KD      (input) INTEGER
046: *          The number of superdiagonals of the matrix A if UPLO = 'U',
047: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
048: *
049: *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
050: *          On entry, the upper or lower triangle of the symmetric band
051: *          matrix A, stored in the first KD+1 rows of the array.  The
052: *          j-th column of A is stored in the j-th column of the array AB
053: *          as follows:
054: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
055: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
056: *
057: *          On exit, AB is overwritten by values generated during the
058: *          reduction to tridiagonal form.  If UPLO = 'U', the first
059: *          superdiagonal and the diagonal of the tridiagonal matrix T
060: *          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
061: *          the diagonal and first subdiagonal of T are returned in the
062: *          first two rows of AB.
063: *
064: *  LDAB    (input) INTEGER
065: *          The leading dimension of the array AB.  LDAB >= KD + 1.
066: *
067: *  W       (output) DOUBLE PRECISION array, dimension (N)
068: *          If INFO = 0, the eigenvalues in ascending order.
069: *
070: *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, N)
071: *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
072: *          eigenvectors of the matrix A, with the i-th column of Z
073: *          holding the eigenvector associated with W(i).
074: *          If JOBZ = 'N', then Z is not referenced.
075: *
076: *  LDZ     (input) INTEGER
077: *          The leading dimension of the array Z.  LDZ >= 1, and if
078: *          JOBZ = 'V', LDZ >= max(1,N).
079: *
080: *  WORK    (workspace/output) DOUBLE PRECISION array,
081: *                                         dimension (LWORK)
082: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
083: *
084: *  LWORK   (input) INTEGER
085: *          The dimension of the array WORK.
086: *          IF N <= 1,                LWORK must be at least 1.
087: *          If JOBZ  = 'N' and N > 2, LWORK must be at least 2*N.
088: *          If JOBZ  = 'V' and N > 2, LWORK must be at least
089: *                         ( 1 + 5*N + 2*N**2 ).
090: *
091: *          If LWORK = -1, then a workspace query is assumed; the routine
092: *          only calculates the optimal sizes of the WORK and IWORK
093: *          arrays, returns these values as the first entries of the WORK
094: *          and IWORK arrays, and no error message related to LWORK or
095: *          LIWORK is issued by XERBLA.
096: *
097: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
098: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
099: *
100: *  LIWORK  (input) INTEGER
101: *          The dimension of the array LIWORK.
102: *          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
103: *          If JOBZ  = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
104: *
105: *          If LIWORK = -1, then a workspace query is assumed; the
106: *          routine only calculates the optimal sizes of the WORK and
107: *          IWORK arrays, returns these values as the first entries of
108: *          the WORK and IWORK arrays, and no error message related to
109: *          LWORK or LIWORK is issued by XERBLA.
110: *
111: *  INFO    (output) INTEGER
112: *          = 0:  successful exit
113: *          < 0:  if INFO = -i, the i-th argument had an illegal value
114: *          > 0:  if INFO = i, the algorithm failed to converge; i
115: *                off-diagonal elements of an intermediate tridiagonal
116: *                form did not converge to zero.
117: *
118: *  =====================================================================
119: *
120: *     .. Parameters ..
121:       DOUBLE PRECISION   ZERO, ONE
122:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
123: *     ..
124: *     .. Local Scalars ..
125:       LOGICAL            LOWER, LQUERY, WANTZ
126:       INTEGER            IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
127:      $                   LLWRK2, LWMIN
128:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
129:      $                   SMLNUM
130: *     ..
131: *     .. External Functions ..
132:       LOGICAL            LSAME
133:       DOUBLE PRECISION   DLAMCH, DLANSB
134:       EXTERNAL           LSAME, DLAMCH, DLANSB
135: *     ..
136: *     .. External Subroutines ..
137:       EXTERNAL           DGEMM, DLACPY, DLASCL, DSBTRD, DSCAL, DSTEDC,
138:      $                   DSTERF, XERBLA
139: *     ..
140: *     .. Intrinsic Functions ..
141:       INTRINSIC          SQRT
142: *     ..
143: *     .. Executable Statements ..
144: *
145: *     Test the input parameters.
146: *
147:       WANTZ = LSAME( JOBZ, 'V' )
148:       LOWER = LSAME( UPLO, 'L' )
149:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
150: *
151:       INFO = 0
152:       IF( N.LE.1 ) THEN
153:          LIWMIN = 1
154:          LWMIN = 1
155:       ELSE
156:          IF( WANTZ ) THEN
157:             LIWMIN = 3 + 5*N
158:             LWMIN = 1 + 5*N + 2*N**2
159:          ELSE
160:             LIWMIN = 1
161:             LWMIN = 2*N
162:          END IF
163:       END IF
164:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
165:          INFO = -1
166:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
167:          INFO = -2
168:       ELSE IF( N.LT.0 ) THEN
169:          INFO = -3
170:       ELSE IF( KD.LT.0 ) THEN
171:          INFO = -4
172:       ELSE IF( LDAB.LT.KD+1 ) THEN
173:          INFO = -6
174:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
175:          INFO = -9
176:       END IF
177: *
178:       IF( INFO.EQ.0 ) THEN
179:          WORK( 1 ) = LWMIN
180:          IWORK( 1 ) = LIWMIN
181: *
182:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
183:             INFO = -11
184:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
185:             INFO = -13
186:          END IF
187:       END IF
188: *
189:       IF( INFO.NE.0 ) THEN
190:          CALL XERBLA( 'DSBEVD', -INFO )
191:          RETURN
192:       ELSE IF( LQUERY ) THEN
193:          RETURN
194:       END IF
195: *
196: *     Quick return if possible
197: *
198:       IF( N.EQ.0 )
199:      $   RETURN
200: *
201:       IF( N.EQ.1 ) THEN
202:          W( 1 ) = AB( 1, 1 )
203:          IF( WANTZ )
204:      $      Z( 1, 1 ) = ONE
205:          RETURN
206:       END IF
207: *
208: *     Get machine constants.
209: *
210:       SAFMIN = DLAMCH( 'Safe minimum' )
211:       EPS = DLAMCH( 'Precision' )
212:       SMLNUM = SAFMIN / EPS
213:       BIGNUM = ONE / SMLNUM
214:       RMIN = SQRT( SMLNUM )
215:       RMAX = SQRT( BIGNUM )
216: *
217: *     Scale matrix to allowable range, if necessary.
218: *
219:       ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
220:       ISCALE = 0
221:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
222:          ISCALE = 1
223:          SIGMA = RMIN / ANRM
224:       ELSE IF( ANRM.GT.RMAX ) THEN
225:          ISCALE = 1
226:          SIGMA = RMAX / ANRM
227:       END IF
228:       IF( ISCALE.EQ.1 ) THEN
229:          IF( LOWER ) THEN
230:             CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
231:          ELSE
232:             CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
233:          END IF
234:       END IF
235: *
236: *     Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
237: *
238:       INDE = 1
239:       INDWRK = INDE + N
240:       INDWK2 = INDWRK + N*N
241:       LLWRK2 = LWORK - INDWK2 + 1
242:       CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
243:      $             WORK( INDWRK ), IINFO )
244: *
245: *     For eigenvalues only, call DSTERF.  For eigenvectors, call SSTEDC.
246: *
247:       IF( .NOT.WANTZ ) THEN
248:          CALL DSTERF( N, W, WORK( INDE ), INFO )
249:       ELSE
250:          CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
251:      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
252:          CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
253:      $               ZERO, WORK( INDWK2 ), N )
254:          CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
255:       END IF
256: *
257: *     If matrix was scaled, then rescale eigenvalues appropriately.
258: *
259:       IF( ISCALE.EQ.1 )
260:      $   CALL DSCAL( N, ONE / SIGMA, W, 1 )
261: *
262:       WORK( 1 ) = LWMIN
263:       IWORK( 1 ) = LIWMIN
264:       RETURN
265: *
266: *     End of DSBEVD
267: *
268:       END
269: