1*> \brief \b CHEGVD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHEGVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chegvd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chegvd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chegvd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
22*                          LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBZ, UPLO
26*       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       REAL               RWORK( * ), W( * )
31*       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CHEGVD computes all the eigenvalues, and optionally, the eigenvectors
41*> of a complex generalized Hermitian-definite eigenproblem, of the form
42*> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
43*> B are assumed to be Hermitian and B is also positive definite.
44*> If eigenvectors are desired, it uses a divide and conquer algorithm.
45*>
46*> The divide and conquer algorithm makes very mild assumptions about
47*> floating point arithmetic. It will work on machines with a guard
48*> digit in add/subtract, or on those binary machines without guard
49*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
50*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
51*> without guard digits, but we know of none.
52*> \endverbatim
53*
54*  Arguments:
55*  ==========
56*
57*> \param[in] ITYPE
58*> \verbatim
59*>          ITYPE is INTEGER
60*>          Specifies the problem type to be solved:
61*>          = 1:  A*x = (lambda)*B*x
62*>          = 2:  A*B*x = (lambda)*x
63*>          = 3:  B*A*x = (lambda)*x
64*> \endverbatim
65*>
66*> \param[in] JOBZ
67*> \verbatim
68*>          JOBZ is CHARACTER*1
69*>          = 'N':  Compute eigenvalues only;
70*>          = 'V':  Compute eigenvalues and eigenvectors.
71*> \endverbatim
72*>
73*> \param[in] UPLO
74*> \verbatim
75*>          UPLO is CHARACTER*1
76*>          = 'U':  Upper triangles of A and B are stored;
77*>          = 'L':  Lower triangles of A and B are stored.
78*> \endverbatim
79*>
80*> \param[in] N
81*> \verbatim
82*>          N is INTEGER
83*>          The order of the matrices A and B.  N >= 0.
84*> \endverbatim
85*>
86*> \param[in,out] A
87*> \verbatim
88*>          A is COMPLEX array, dimension (LDA, N)
89*>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
90*>          leading N-by-N upper triangular part of A contains the
91*>          upper triangular part of the matrix A.  If UPLO = 'L',
92*>          the leading N-by-N lower triangular part of A contains
93*>          the lower triangular part of the matrix A.
94*>
95*>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
96*>          matrix Z of eigenvectors.  The eigenvectors are normalized
97*>          as follows:
98*>          if ITYPE = 1 or 2, Z**H*B*Z = I;
99*>          if ITYPE = 3, Z**H*inv(B)*Z = I.
100*>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
101*>          or the lower triangle (if UPLO='L') of A, including the
102*>          diagonal, is destroyed.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*>          LDA is INTEGER
108*>          The leading dimension of the array A.  LDA >= max(1,N).
109*> \endverbatim
110*>
111*> \param[in,out] B
112*> \verbatim
113*>          B is COMPLEX array, dimension (LDB, N)
114*>          On entry, the Hermitian matrix B.  If UPLO = 'U', the
115*>          leading N-by-N upper triangular part of B contains the
116*>          upper triangular part of the matrix B.  If UPLO = 'L',
117*>          the leading N-by-N lower triangular part of B contains
118*>          the lower triangular part of the matrix B.
119*>
120*>          On exit, if INFO <= N, the part of B containing the matrix is
121*>          overwritten by the triangular factor U or L from the Cholesky
122*>          factorization B = U**H*U or B = L*L**H.
123*> \endverbatim
124*>
125*> \param[in] LDB
126*> \verbatim
127*>          LDB is INTEGER
128*>          The leading dimension of the array B.  LDB >= max(1,N).
129*> \endverbatim
130*>
131*> \param[out] W
132*> \verbatim
133*>          W is REAL array, dimension (N)
134*>          If INFO = 0, the eigenvalues in ascending order.
135*> \endverbatim
136*>
137*> \param[out] WORK
138*> \verbatim
139*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
140*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
141*> \endverbatim
142*>
143*> \param[in] LWORK
144*> \verbatim
145*>          LWORK is INTEGER
146*>          The length of the array WORK.
147*>          If N <= 1,                LWORK >= 1.
148*>          If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
149*>          If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.
150*>
151*>          If LWORK = -1, then a workspace query is assumed; the routine
152*>          only calculates the optimal sizes of the WORK, RWORK and
153*>          IWORK arrays, returns these values as the first entries of
154*>          the WORK, RWORK and IWORK arrays, and no error message
155*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
156*> \endverbatim
157*>
158*> \param[out] RWORK
159*> \verbatim
160*>          RWORK is REAL array, dimension (MAX(1,LRWORK))
161*>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
162*> \endverbatim
163*>
164*> \param[in] LRWORK
165*> \verbatim
166*>          LRWORK is INTEGER
167*>          The dimension of the array RWORK.
168*>          If N <= 1,                LRWORK >= 1.
169*>          If JOBZ  = 'N' and N > 1, LRWORK >= N.
170*>          If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
171*>
172*>          If LRWORK = -1, then a workspace query is assumed; the
173*>          routine only calculates the optimal sizes of the WORK, RWORK
174*>          and IWORK arrays, returns these values as the first entries
175*>          of the WORK, RWORK and IWORK arrays, and no error message
176*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
177*> \endverbatim
178*>
179*> \param[out] IWORK
180*> \verbatim
181*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
182*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
183*> \endverbatim
184*>
185*> \param[in] LIWORK
186*> \verbatim
187*>          LIWORK is INTEGER
188*>          The dimension of the array IWORK.
189*>          If N <= 1,                LIWORK >= 1.
190*>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
191*>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
192*>
193*>          If LIWORK = -1, then a workspace query is assumed; the
194*>          routine only calculates the optimal sizes of the WORK, RWORK
195*>          and IWORK arrays, returns these values as the first entries
196*>          of the WORK, RWORK and IWORK arrays, and no error message
197*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
198*> \endverbatim
199*>
200*> \param[out] INFO
201*> \verbatim
202*>          INFO is INTEGER
203*>          = 0:  successful exit
204*>          < 0:  if INFO = -i, the i-th argument had an illegal value
205*>          > 0:  CPOTRF or CHEEVD returned an error code:
206*>             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
207*>                    failed to converge; i off-diagonal elements of an
208*>                    intermediate tridiagonal form did not converge to
209*>                    zero;
210*>                    if INFO = i and JOBZ = 'V', then the algorithm
211*>                    failed to compute an eigenvalue while working on
212*>                    the submatrix lying in rows and columns INFO/(N+1)
213*>                    through mod(INFO,N+1);
214*>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
215*>                    minor of order i of B is not positive definite.
216*>                    The factorization of B could not be completed and
217*>                    no eigenvalues or eigenvectors were computed.
218*> \endverbatim
219*
220*  Authors:
221*  ========
222*
223*> \author Univ. of Tennessee
224*> \author Univ. of California Berkeley
225*> \author Univ. of Colorado Denver
226*> \author NAG Ltd.
227*
228*> \ingroup complexHEeigen
229*
230*> \par Further Details:
231*  =====================
232*>
233*> \verbatim
234*>
235*>  Modified so that no backsubstitution is performed if CHEEVD fails to
236*>  converge (NEIG in old code could be greater than N causing out of
237*>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
238*>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
239*> \endverbatim
240*
241*> \par Contributors:
242*  ==================
243*>
244*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
245*>
246*  =====================================================================
247      SUBROUTINE CHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
248     $                   LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
249*
250*  -- LAPACK driver routine --
251*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
252*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253*
254*     .. Scalar Arguments ..
255      CHARACTER          JOBZ, UPLO
256      INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
257*     ..
258*     .. Array Arguments ..
259      INTEGER            IWORK( * )
260      REAL               RWORK( * ), W( * )
261      COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
262*     ..
263*
264*  =====================================================================
265*
266*     .. Parameters ..
267      COMPLEX            CONE
268      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
269*     ..
270*     .. Local Scalars ..
271      LOGICAL            LQUERY, UPPER, WANTZ
272      CHARACTER          TRANS
273      INTEGER            LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
274*     ..
275*     .. External Functions ..
276      LOGICAL            LSAME
277      EXTERNAL           LSAME
278*     ..
279*     .. External Subroutines ..
280      EXTERNAL           CHEEVD, CHEGST, CPOTRF, CTRMM, CTRSM, XERBLA
281*     ..
282*     .. Intrinsic Functions ..
283      INTRINSIC          MAX, REAL
284*     ..
285*     .. Executable Statements ..
286*
287*     Test the input parameters.
288*
289      WANTZ = LSAME( JOBZ, 'V' )
290      UPPER = LSAME( UPLO, 'U' )
291      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
292*
293      INFO = 0
294      IF( N.LE.1 ) THEN
295         LWMIN = 1
296         LRWMIN = 1
297         LIWMIN = 1
298      ELSE IF( WANTZ ) THEN
299         LWMIN = 2*N + N*N
300         LRWMIN = 1 + 5*N + 2*N*N
301         LIWMIN = 3 + 5*N
302      ELSE
303         LWMIN = N + 1
304         LRWMIN = N
305         LIWMIN = 1
306      END IF
307      LOPT = LWMIN
308      LROPT = LRWMIN
309      LIOPT = LIWMIN
310      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
311         INFO = -1
312      ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
313         INFO = -2
314      ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
315         INFO = -3
316      ELSE IF( N.LT.0 ) THEN
317         INFO = -4
318      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
319         INFO = -6
320      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
321         INFO = -8
322      END IF
323*
324      IF( INFO.EQ.0 ) THEN
325         WORK( 1 ) = LOPT
326         RWORK( 1 ) = LROPT
327         IWORK( 1 ) = LIOPT
328*
329         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
330            INFO = -11
331         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
332            INFO = -13
333         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
334            INFO = -15
335         END IF
336      END IF
337*
338      IF( INFO.NE.0 ) THEN
339         CALL XERBLA( 'CHEGVD', -INFO )
340         RETURN
341      ELSE IF( LQUERY ) THEN
342         RETURN
343      END IF
344*
345*     Quick return if possible
346*
347      IF( N.EQ.0 )
348     $   RETURN
349*
350*     Form a Cholesky factorization of B.
351*
352      CALL CPOTRF( UPLO, N, B, LDB, INFO )
353      IF( INFO.NE.0 ) THEN
354         INFO = N + INFO
355         RETURN
356      END IF
357*
358*     Transform problem to standard eigenvalue problem and solve.
359*
360      CALL CHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
361      CALL CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK,
362     $             IWORK, LIWORK, INFO )
363      LOPT = MAX( REAL( LOPT ), REAL( WORK( 1 ) ) )
364      LROPT = MAX( REAL( LROPT ), REAL( RWORK( 1 ) ) )
365      LIOPT = MAX( REAL( LIOPT ), REAL( IWORK( 1 ) ) )
366*
367      IF( WANTZ .AND. INFO.EQ.0 ) THEN
368*
369*        Backtransform eigenvectors to the original problem.
370*
371         IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
372*
373*           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
374*           backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
375*
376            IF( UPPER ) THEN
377               TRANS = 'N'
378            ELSE
379               TRANS = 'C'
380            END IF
381*
382            CALL CTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
383     $                  B, LDB, A, LDA )
384*
385         ELSE IF( ITYPE.EQ.3 ) THEN
386*
387*           For B*A*x=(lambda)*x;
388*           backtransform eigenvectors: x = L*y or U**H *y
389*
390            IF( UPPER ) THEN
391               TRANS = 'C'
392            ELSE
393               TRANS = 'N'
394            END IF
395*
396            CALL CTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
397     $                  B, LDB, A, LDA )
398         END IF
399      END IF
400*
401      WORK( 1 ) = LOPT
402      RWORK( 1 ) = LROPT
403      IWORK( 1 ) = LIOPT
404*
405      RETURN
406*
407*     End of CHEGVD
408*
409      END
410