Eigenvalues and Eigenvectors of a Symmetric Matrix in Banded Storage (Expert Driver)

Calling Sequence

CALL DSBEVX	(JOBZ, RANGE, UPLO, N, NDIAG, DA, LDA, DQ, LDQ, DVL, DVU, IL, IU, DABTOL, NFOUND, DW, DZ, LDZ, DWORK, IWORK2, IFAIL, INFO)
CALL SSBEVX	(JOBZ, RANGE, UPLO, N, NDIAG, SA, LDA, SQ, LDQ, SVL, SVU, IL, IU, SABTOL, NFOUND, SW, SZ, LDZ, SWORK, IWORK2, IFAIL, INFO)
void dsbevx	(char jobz, char range, char uplo, int n, int ndiag, double *da, int lda, double *dq, int ldq, double dvl, double dvu, int il, int iu, double dabtol, int *nfound, double *dw, double *dz, int ldz, int *ifail, int *info)
void ssbevx	(char jobz, char range, char uplo, int n, int ndiag, float *sa, int lda, float *sq, int ldq, float svl, float svu, int il, int iu, float sabtol, int *nfound, float *sw, float *sz, int ldz, int *ifail, int *info)

Arguments

JOBZ	Indicates whether to compute eigenvalues only or to compute both eigenvalues and eigenvectors. The legal values for JOBZ are listed below. Any values not listed below are illegal.
	'N' or 'n'	Compute eigenvalues only.
	'V' or 'v'	Compute both eigenvalues and eigenvectors.
RANGE	Indicates whether to compute all eigenvalues or some subset. Legal values for RANGE are shown below. Any values not shown below are illegal.
	'A' or 'a'	All eigenvalues will be computed.
	'I' or 'i'	The ILth through IUth eigenvalues will be computed; see the descriptions of the IL and IU arguments below.
	'V' or 'v'	All eigenvalues in the half-open interval (VL,VU] will be computed; see the descriptions of the VL and VU arguments below.
UPLO	Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.
	'U' or 'u'	xA contains the upper triangle.
	'L' or 'l'	xA contains the lower triangle.
N	Order of the matrix A. N 0.
NDIAG	Number of superdiagonals or subdiagonals of the matrix A. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.
	If UPLO = 'U' or 'u', NDIAG is the number of superdiagonals.
	If UPLO = 'L' or 'l', NDIAG is the number of subdiagonals.
xA	On entry, the upper or lower triangle of the matrix A. On exit, A has been overwritten.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA NDIAG + 1.
xQ	On exit, if JOBZ = 'V' or 'v' then Q contains the N×N orthogonal matrix used in the reduction to tridiagonal form. If JOBZ = 'N' or 'n' then Q is not used.
LDQ	Leading dimension of the array Q as specified in a dimension or type statement. LDQ max(1, N).
xVL, xVU	If RANGE = 'V' or 'v' then only eigenvalues in the half-open interval (VL,VU] will be computed. Otherwise not used.
IL, IU	If RANGE = 'I' or 'i' then only the ILth through IUth eigenvalues will be computed. Otherwise not used. 1 IL IU N.
xABTOL	The absolute error tolerance for the eigenvalues. An approximate eigenvalue is considered to have converged when it lies in an interval [a,b] of width less than or equal to ABTOL + × max(|a|, |b|) where is machine precision. If ABTOL 0 then × ||T||1 will be used in its place, where ||T||1 is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
NFOUND	On exit, the total number of eigenvalues found. 0 NFOUND N. If RANGE = 'A' or 'a' then NFOUND = N. If RANGE = 'I' or 'i' then NFOUND = IU - IL + 1.
xW	On exit, W(1:NFOUND) contains the computed eigenvalues in ascending order.
xZ	On exit, if JOBZ = 'V' or 'v' then the first NFOUND columns of Z contain the orthonormal eigenvectors corresponding to the selected eigenvalues. If an eigenvector fails to converge then the corresponding column of Z contains an approximation to that eigenvector and its index is returned in IFAIL. If JOBZ = 'N' or 'n' then Z is not used.
	Note that Z must have at least NFOUND columns. If RANGE = 'V' or 'v' then NFOUND may not be known in advance and an upper bound must be used. For A(LDZ,K) where NFOUND > K, the eigenvectors will be stored into memory with unpredictable results.
LDZ	Leading dimension of the array Z as specified in a dimension or type statement. LDZ 1. If JOBZ = 'V' or 'v' then LDZ max(1, N).
xWORK	Scratch array with a dimension of 7 × N.
IWORK2	Scratch array with a dimension of 5 × N.
IFAIL	On exit, an N-element array that indicates which eigenvectors that failed to converge. If JOBZ = 'V' or 'v' and INFO = 0 then all N elements contain zero. If INFO > 0 then IFAIL contains the indices of the eigenvectors that failed to converge. If JOBZ = 'N' or 'n' then IFAIL is not used.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	INFO > 0	There were i eigenvectors, where i = INFO, that failed to converge. Their indices are stored in IFAIL.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDEVEC, LDIWRK, LDQ, LDWORK, N, NDIAG

PARAMETER (N = 4)

PARAMETER (NDIAG = 1)

PARAMETER (LDA = NDIAG + 1)

PARAMETER (LDEVEC = N)

PARAMETER (LDIWRK = 5 * N)

PARAMETER (LDWORK = 7 * N)

PARAMETER (LDQ = N)

C

DOUBLE PRECISION A(LDA,N), EVALS(N), EVECS(LDEVEC,N), TEMP

DOUBLE PRECISION Q(LDQ,N), WORK(LDWORK)

INTEGER ICOL, IFAIL(N), INFO, IROW, ITEMP

INTEGER IWORK(LDIWRK), NFOUND

C

EXTERNAL DSBEVX

INTRINSIC ABS

C

C Initialize the array A to store the symmetric tridiagonal

C matrix A shown below.

C

C 1 0 0 0

C A = 0 2 1 0

C 0 1 2 0

C 0 0 0 1

C

DATA A / 8D8, 1.0D0, 0.0D0, 2.0D0, 1.0D0, 2.0D0, 0.0D0,

$ 1.0D0 /

C

PRINT 1000

PRINT 1010, A(2,1), A(1,2), 0.0D0, 0.0D0

PRINT 1010, A(1,2), A(2,2), A(1,3), 0.0D0

PRINT 1010, 0.0D0, A(1,3), A(2,3), A(1,4)

PRINT 1010, 0.0D0, 0.0D0, A(1,4), A(2,4)

PRINT 1020

DO 100, IROW = 1, LDA

PRINT 1010, (A(IROW,ICOL), ICOL = 1, N)

100 CONTINUE

C

C Compute the eigenvalues and right eigenvectors of A.

C

CALL DSBEVX ('VECTORS AND EIGENVALUES', 'ALL EIGENVALUES',

$ 'UPPER TRIANGLE A STORED', N, NDIAG, A, LDA,

$ Q, LDQ, TEMP, TEMP, ITEMP, ITEMP, 0.0D0,

$ NFOUND, EVALS, EVECS, LDEVEC, WORK, IWORK,

$ IFAIL, INFO)

IF (INFO .NE. 0) THEN

IF (INFO .LT. 0) THEN

PRINT 1030, ABS(INFO)

STOP 1

ELSE

PRINT 1040, INFO

STOP 2

END IF

END IF

C

C Print the eigenvalues and eigenvectors.

C

PRINT 1050

DO 120, ICOL = 1, N

PRINT 1060, EVALS(ICOL), (EVECS(IROW,ICOL), IROW = 1, N)

120 CONTINUE

C

1000 FORMAT (1X, 'A in full form:')

1010 FORMAT (4(3X, F9.6))

1020 FORMAT (/1X,

$ 'A in symmetric banded form: '

$ ' (* in unused elements)')

1030 FORMAT (/1X, 'Illegal argument to DSBEVX, argument #', I2)

1040 FORMAT (/1X, 'Convergence failure, INFO = ', I2)

1050 FORMAT (/1X, 'Eigenvalue', 10X, 'Eigenvector**T')

1060 FORMAT (1X, F7.2, 6X, '[', 3(F5.3, ', '), F5.3, ']')

C

END

Sample Output

A in full form:

1.000000 0.000000 0.000000 0.000000

0.000000 2.000000 1.000000 0.000000

0.000000 1.000000 2.000000 0.000000

0.000000 0.000000 0.000000 1.000000

A in symmetric banded form: (* in unused elements)

********* 0.000000 1.000000 0.000000

1.000000 2.000000 2.000000 1.000000

Eigenvalue Eigenvector**T

1.00 [1.000, 0.000, 0.000, 0.000]

1.00 [0.000, -.707, 0.707, 0.000]

1.00 [0.000, 0.000, 0.000, 1.000]

3.00 [0.000, 0.707, 0.707, 0.000]