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