001:       SUBROUTINE CHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
002:      $                   LWORK, RWORK, LRWORK, 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, LRWORK, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IWORK( * )
014:       REAL               RWORK( * ), W( * )
015:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CHEGVD computes all the eigenvalues, and optionally, the eigenvectors
022: *  of a complex generalized Hermitian-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 Hermitian 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) COMPLEX array, dimension (LDA, N)
055: *          On entry, the Hermitian 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**H*B*Z = I;
065: *          if ITYPE = 3, Z**H*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) COMPLEX array, dimension (LDB, N)
074: *          On entry, the Hermitian 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**H*U or B = L*L**H.
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) COMPLEX array, dimension (MAX(1,LWORK))
091: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
092: *
093: *  LWORK   (input) INTEGER
094: *          The length of the array WORK.
095: *          If N <= 1,                LWORK >= 1.
096: *          If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
097: *          If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.
098: *
099: *          If LWORK = -1, then a workspace query is assumed; the routine
100: *          only calculates the optimal sizes of the WORK, RWORK and
101: *          IWORK arrays, returns these values as the first entries of
102: *          the WORK, RWORK and IWORK arrays, and no error message
103: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
104: *
105: *  RWORK   (workspace/output) REAL array, dimension (MAX(1,LRWORK))
106: *          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
107: *
108: *  LRWORK  (input) INTEGER
109: *          The dimension of the array RWORK.
110: *          If N <= 1,                LRWORK >= 1.
111: *          If JOBZ  = 'N' and N > 1, LRWORK >= N.
112: *          If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
113: *
114: *          If LRWORK = -1, then a workspace query is assumed; the
115: *          routine only calculates the optimal sizes of the WORK, RWORK
116: *          and IWORK arrays, returns these values as the first entries
117: *          of the WORK, RWORK and IWORK arrays, and no error message
118: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
119: *
120: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
121: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
122: *
123: *  LIWORK  (input) INTEGER
124: *          The dimension of the array IWORK.
125: *          If N <= 1,                LIWORK >= 1.
126: *          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
127: *          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
128: *
129: *          If LIWORK = -1, then a workspace query is assumed; the
130: *          routine only calculates the optimal sizes of the WORK, RWORK
131: *          and IWORK arrays, returns these values as the first entries
132: *          of the WORK, RWORK and IWORK arrays, and no error message
133: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
134: *
135: *  INFO    (output) INTEGER
136: *          = 0:  successful exit
137: *          < 0:  if INFO = -i, the i-th argument had an illegal value
138: *          > 0:  CPOTRF or CHEEVD returned an error code:
139: *             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
140: *                    failed to converge; i off-diagonal elements of an
141: *                    intermediate tridiagonal form did not converge to
142: *                    zero;
143: *                    if INFO = i and JOBZ = 'V', then the algorithm
144: *                    failed to compute an eigenvalue while working on
145: *                    the submatrix lying in rows and columns INFO/(N+1)
146: *                    through mod(INFO,N+1);
147: *             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
148: *                    minor of order i of B is not positive definite.
149: *                    The factorization of B could not be completed and
150: *                    no eigenvalues or eigenvectors were computed.
151: *
152: *  Further Details
153: *  ===============
154: *
155: *  Based on contributions by
156: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
157: *
158: *  Modified so that no backsubstitution is performed if CHEEVD fails to
159: *  converge (NEIG in old code could be greater than N causing out of
160: *  bounds reference to A - reported by Ralf Meyer).  Also corrected the
161: *  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
162: *  =====================================================================
163: *
164: *     .. Parameters ..
165:       COMPLEX            CONE
166:       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
167: *     ..
168: *     .. Local Scalars ..
169:       LOGICAL            LQUERY, UPPER, WANTZ
170:       CHARACTER          TRANS
171:       INTEGER            LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
172: *     ..
173: *     .. External Functions ..
174:       LOGICAL            LSAME
175:       EXTERNAL           LSAME
176: *     ..
177: *     .. External Subroutines ..
178:       EXTERNAL           CHEEVD, CHEGST, CPOTRF, CTRMM, CTRSM, XERBLA
179: *     ..
180: *     .. Intrinsic Functions ..
181:       INTRINSIC          MAX, REAL
182: *     ..
183: *     .. Executable Statements ..
184: *
185: *     Test the input parameters.
186: *
187:       WANTZ = LSAME( JOBZ, 'V' )
188:       UPPER = LSAME( UPLO, 'U' )
189:       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
190: *
191:       INFO = 0
192:       IF( N.LE.1 ) THEN
193:          LWMIN = 1
194:          LRWMIN = 1
195:          LIWMIN = 1
196:       ELSE IF( WANTZ ) THEN
197:          LWMIN = 2*N + N*N
198:          LRWMIN = 1 + 5*N + 2*N*N
199:          LIWMIN = 3 + 5*N
200:       ELSE
201:          LWMIN = N + 1
202:          LRWMIN = N
203:          LIWMIN = 1
204:       END IF
205:       LOPT = LWMIN
206:       LROPT = LRWMIN
207:       LIOPT = LIWMIN
208:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
209:          INFO = -1
210:       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
211:          INFO = -2
212:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
213:          INFO = -3
214:       ELSE IF( N.LT.0 ) THEN
215:          INFO = -4
216:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
217:          INFO = -6
218:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
219:          INFO = -8
220:       END IF
221: *
222:       IF( INFO.EQ.0 ) THEN
223:          WORK( 1 ) = LOPT
224:          RWORK( 1 ) = LROPT
225:          IWORK( 1 ) = LIOPT
226: *
227:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
228:             INFO = -11
229:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
230:             INFO = -13
231:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
232:             INFO = -15
233:          END IF
234:       END IF
235: *
236:       IF( INFO.NE.0 ) THEN
237:          CALL XERBLA( 'CHEGVD', -INFO )
238:          RETURN
239:       ELSE IF( LQUERY ) THEN
240:          RETURN
241:       END IF
242: *
243: *     Quick return if possible
244: *
245:       IF( N.EQ.0 )
246:      $   RETURN
247: *
248: *     Form a Cholesky factorization of B.
249: *
250:       CALL CPOTRF( UPLO, N, B, LDB, INFO )
251:       IF( INFO.NE.0 ) THEN
252:          INFO = N + INFO
253:          RETURN
254:       END IF
255: *
256: *     Transform problem to standard eigenvalue problem and solve.
257: *
258:       CALL CHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
259:       CALL CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK,
260:      $             IWORK, LIWORK, INFO )
261:       LOPT = MAX( REAL( LOPT ), REAL( WORK( 1 ) ) )
262:       LROPT = MAX( REAL( LROPT ), REAL( RWORK( 1 ) ) )
263:       LIOPT = MAX( REAL( LIOPT ), REAL( IWORK( 1 ) ) )
264: *
265:       IF( WANTZ .AND. INFO.EQ.0 ) THEN
266: *
267: *        Backtransform eigenvectors to the original problem.
268: *
269:          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
270: *
271: *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
272: *           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
273: *
274:             IF( UPPER ) THEN
275:                TRANS = 'N'
276:             ELSE
277:                TRANS = 'C'
278:             END IF
279: *
280:             CALL CTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
281:      $                  B, LDB, A, LDA )
282: *
283:          ELSE IF( ITYPE.EQ.3 ) THEN
284: *
285: *           For B*A*x=(lambda)*x;
286: *           backtransform eigenvectors: x = L*y or U'*y
287: *
288:             IF( UPPER ) THEN
289:                TRANS = 'C'
290:             ELSE
291:                TRANS = 'N'
292:             END IF
293: *
294:             CALL CTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
295:      $                  B, LDB, A, LDA )
296:          END IF
297:       END IF
298: *
299:       WORK( 1 ) = LOPT
300:       RWORK( 1 ) = LROPT
301:       IWORK( 1 ) = LIOPT
302: *
303:       RETURN
304: *
305: *     End of CHEGVD
306: *
307:       END
308: