001:       SUBROUTINE SSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, 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, ITYPE, LDA, LDB, LIWORK, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SSYGVD computes all the eigenvalues, and optionally, the eigenvectors
022: *  of a real generalized symmetric-definite eigenproblem, of the form
023: *  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
024: *  B are assumed to be symmetric and B is also positive definite.
025: *  If eigenvectors are desired, it uses a divide and conquer algorithm.
026: *
027: *  The divide and conquer algorithm makes very mild assumptions about
028: *  floating point arithmetic. It will work on machines with a guard
029: *  digit in add/subtract, or on those binary machines without guard
030: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
031: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
032: *  without guard digits, but we know of none.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  ITYPE   (input) INTEGER
038: *          Specifies the problem type to be solved:
039: *          = 1:  A*x = (lambda)*B*x
040: *          = 2:  A*B*x = (lambda)*x
041: *          = 3:  B*A*x = (lambda)*x
042: *
043: *  JOBZ    (input) CHARACTER*1
044: *          = 'N':  Compute eigenvalues only;
045: *          = 'V':  Compute eigenvalues and eigenvectors.
046: *
047: *  UPLO    (input) CHARACTER*1
048: *          = 'U':  Upper triangles of A and B are stored;
049: *          = 'L':  Lower triangles of A and B are stored.
050: *
051: *  N       (input) INTEGER
052: *          The order of the matrices A and B.  N >= 0.
053: *
054: *  A       (input/output) REAL array, dimension (LDA, N)
055: *          On entry, the symmetric matrix A.  If UPLO = 'U', the
056: *          leading N-by-N upper triangular part of A contains the
057: *          upper triangular part of the matrix A.  If UPLO = 'L',
058: *          the leading N-by-N lower triangular part of A contains
059: *          the lower triangular part of the matrix A.
060: *
061: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
062: *          matrix Z of eigenvectors.  The eigenvectors are normalized
063: *          as follows:
064: *          if ITYPE = 1 or 2, Z**T*B*Z = I;
065: *          if ITYPE = 3, Z**T*inv(B)*Z = I.
066: *          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
067: *          or the lower triangle (if UPLO='L') of A, including the
068: *          diagonal, is destroyed.
069: *
070: *  LDA     (input) INTEGER
071: *          The leading dimension of the array A.  LDA >= max(1,N).
072: *
073: *  B       (input/output) REAL array, dimension (LDB, N)
074: *          On entry, the symmetric matrix B.  If UPLO = 'U', the
075: *          leading N-by-N upper triangular part of B contains the
076: *          upper triangular part of the matrix B.  If UPLO = 'L',
077: *          the leading N-by-N lower triangular part of B contains
078: *          the lower triangular part of the matrix B.
079: *
080: *          On exit, if INFO <= N, the part of B containing the matrix is
081: *          overwritten by the triangular factor U or L from the Cholesky
082: *          factorization B = U**T*U or B = L*L**T.
083: *
084: *  LDB     (input) INTEGER
085: *          The leading dimension of the array B.  LDB >= max(1,N).
086: *
087: *  W       (output) REAL array, dimension (N)
088: *          If INFO = 0, the eigenvalues in ascending order.
089: *
090: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
091: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
092: *
093: *  LWORK   (input) INTEGER
094: *          The dimension of the array WORK.
095: *          If N <= 1,               LWORK >= 1.
096: *          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
097: *          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
098: *
099: *          If LWORK = -1, then a workspace query is assumed; the routine
100: *          only calculates the optimal sizes of the WORK and IWORK
101: *          arrays, returns these values as the first entries of the WORK
102: *          and IWORK arrays, and no error message related to LWORK or
103: *          LIWORK is issued by XERBLA.
104: *
105: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
106: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
107: *
108: *  LIWORK  (input) INTEGER
109: *          The dimension of the array IWORK.
110: *          If N <= 1,                LIWORK >= 1.
111: *          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
112: *          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
113: *
114: *          If LIWORK = -1, then a workspace query is assumed; the
115: *          routine only calculates the optimal sizes of the WORK and
116: *          IWORK arrays, returns these values as the first entries of
117: *          the WORK and IWORK arrays, and no error message related to
118: *          LWORK or LIWORK is issued by XERBLA.
119: *
120: *  INFO    (output) INTEGER
121: *          = 0:  successful exit
122: *          < 0:  if INFO = -i, the i-th argument had an illegal value
123: *          > 0:  SPOTRF or SSYEVD returned an error code:
124: *             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
125: *                    failed to converge; i off-diagonal elements of an
126: *                    intermediate tridiagonal form did not converge to
127: *                    zero;
128: *                    if INFO = i and JOBZ = 'V', then the algorithm
129: *                    failed to compute an eigenvalue while working on
130: *                    the submatrix lying in rows and columns INFO/(N+1)
131: *                    through mod(INFO,N+1);
132: *             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
133: *                    minor of order i of B is not positive definite.
134: *                    The factorization of B could not be completed and
135: *                    no eigenvalues or eigenvectors were computed.
136: *
137: *  Further Details
138: *  ===============
139: *
140: *  Based on contributions by
141: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
142: *
143: *  Modified so that no backsubstitution is performed if SSYEVD fails to
144: *  converge (NEIG in old code could be greater than N causing out of
145: *  bounds reference to A - reported by Ralf Meyer).  Also corrected the
146: *  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
147: *  =====================================================================
148: *
149: *     .. Parameters ..
150:       REAL               ONE
151:       PARAMETER          ( ONE = 1.0E+0 )
152: *     ..
153: *     .. Local Scalars ..
154:       LOGICAL            LQUERY, UPPER, WANTZ
155:       CHARACTER          TRANS
156:       INTEGER            LIOPT, LIWMIN, LOPT, LWMIN
157: *     ..
158: *     .. External Functions ..
159:       LOGICAL            LSAME
160:       EXTERNAL           LSAME
161: *     ..
162: *     .. External Subroutines ..
163:       EXTERNAL           SPOTRF, SSYEVD, SSYGST, STRMM, STRSM, XERBLA
164: *     ..
165: *     .. Intrinsic Functions ..
166:       INTRINSIC          MAX, REAL
167: *     ..
168: *     .. Executable Statements ..
169: *
170: *     Test the input parameters.
171: *
172:       WANTZ = LSAME( JOBZ, 'V' )
173:       UPPER = LSAME( UPLO, 'U' )
174:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
175: *
176:       INFO = 0
177:       IF( N.LE.1 ) THEN
178:          LIWMIN = 1
179:          LWMIN = 1
180:       ELSE IF( WANTZ ) THEN
181:          LIWMIN = 3 + 5*N
182:          LWMIN = 1 + 6*N + 2*N**2
183:       ELSE
184:          LIWMIN = 1
185:          LWMIN = 2*N + 1
186:       END IF
187:       LOPT = LWMIN
188:       LIOPT = LIWMIN
189:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
190:          INFO = -1
191:       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
192:          INFO = -2
193:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
194:          INFO = -3
195:       ELSE IF( N.LT.0 ) THEN
196:          INFO = -4
197:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
198:          INFO = -6
199:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
200:          INFO = -8
201:       END IF
202: *
203:       IF( INFO.EQ.0 ) THEN
204:          WORK( 1 ) = LOPT
205:          IWORK( 1 ) = LIOPT
206: *
207:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
208:             INFO = -11
209:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
210:             INFO = -13
211:          END IF
212:       END IF
213: *
214:       IF( INFO.NE.0 ) THEN
215:          CALL XERBLA( 'SSYGVD', -INFO )
216:          RETURN
217:       ELSE IF( LQUERY ) THEN
218:          RETURN
219:       END IF
220: *
221: *     Quick return if possible
222: *
223:       IF( N.EQ.0 )
224:      $   RETURN
225: *
226: *     Form a Cholesky factorization of B.
227: *
228:       CALL SPOTRF( UPLO, N, B, LDB, INFO )
229:       IF( INFO.NE.0 ) THEN
230:          INFO = N + INFO
231:          RETURN
232:       END IF
233: *
234: *     Transform problem to standard eigenvalue problem and solve.
235: *
236:       CALL SSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
237:       CALL SSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
238:      $             INFO )
239:       LOPT = MAX( REAL( LOPT ), REAL( WORK( 1 ) ) )
240:       LIOPT = MAX( REAL( LIOPT ), REAL( IWORK( 1 ) ) )
241: *
242:       IF( WANTZ .AND. INFO.EQ.0 ) THEN
243: *
244: *        Backtransform eigenvectors to the original problem.
245: *
246:          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
247: *
248: *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
249: *           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
250: *
251:             IF( UPPER ) THEN
252:                TRANS = 'N'
253:             ELSE
254:                TRANS = 'T'
255:             END IF
256: *
257:             CALL STRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
258:      $                  B, LDB, A, LDA )
259: *
260:          ELSE IF( ITYPE.EQ.3 ) THEN
261: *
262: *           For B*A*x=(lambda)*x;
263: *           backtransform eigenvectors: x = L*y or U'*y
264: *
265:             IF( UPPER ) THEN
266:                TRANS = 'T'
267:             ELSE
268:                TRANS = 'N'
269:             END IF
270: *
271:             CALL STRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
272:      $                  B, LDB, A, LDA )
273:          END IF
274:       END IF
275: *
276:       WORK( 1 ) = LOPT
277:       IWORK( 1 ) = LIOPT
278: *
279:       RETURN
280: *
281: *     End of SSYGVD
282: *
283:       END
284: