1       SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
  2      $                   LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
  3 *
  4 *  -- LAPACK driver routine (version 3.3.1) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *  -- April 2011                                                      --
  8 *
  9 *     .. Scalar Arguments ..
 10       CHARACTER          JOBZ, UPLO
 11       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
 12 *     ..
 13 *     .. Array Arguments ..
 14       INTEGER            IWORK( * )
 15       DOUBLE PRECISION   RWORK( * ), W( * )
 16       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 17 *     ..
 18 *
 19 *  Purpose
 20 *  =======
 21 *
 22 *  ZHEGVD computes all the eigenvalues, and optionally, the eigenvectors
 23 *  of a complex generalized Hermitian-definite eigenproblem, of the form
 24 *  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
 25 *  B are assumed to be Hermitian and B is also positive definite.
 26 *  If eigenvectors are desired, it uses a divide and conquer algorithm.
 27 *
 28 *  The divide and conquer algorithm makes very mild assumptions about
 29 *  floating point arithmetic. It will work on machines with a guard
 30 *  digit in add/subtract, or on those binary machines without guard
 31 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
 32 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
 33 *  without guard digits, but we know of none.
 34 *
 35 *  Arguments
 36 *  =========
 37 *
 38 *  ITYPE   (input) INTEGER
 39 *          Specifies the problem type to be solved:
 40 *          = 1:  A*x = (lambda)*B*x
 41 *          = 2:  A*B*x = (lambda)*x
 42 *          = 3:  B*A*x = (lambda)*x
 43 *
 44 *  JOBZ    (input) CHARACTER*1
 45 *          = 'N':  Compute eigenvalues only;
 46 *          = 'V':  Compute eigenvalues and eigenvectors.
 47 *
 48 *  UPLO    (input) CHARACTER*1
 49 *          = 'U':  Upper triangles of A and B are stored;
 50 *          = 'L':  Lower triangles of A and B are stored.
 51 *
 52 *  N       (input) INTEGER
 53 *          The order of the matrices A and B.  N >= 0.
 54 *
 55 *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
 56 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the
 57 *          leading N-by-N upper triangular part of A contains the
 58 *          upper triangular part of the matrix A.  If UPLO = 'L',
 59 *          the leading N-by-N lower triangular part of A contains
 60 *          the lower triangular part of the matrix A.
 61 *
 62 *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
 63 *          matrix Z of eigenvectors.  The eigenvectors are normalized
 64 *          as follows:
 65 *          if ITYPE = 1 or 2, Z**H*B*Z = I;
 66 *          if ITYPE = 3, Z**H*inv(B)*Z = I.
 67 *          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
 68 *          or the lower triangle (if UPLO='L') of A, including the
 69 *          diagonal, is destroyed.
 70 *
 71 *  LDA     (input) INTEGER
 72 *          The leading dimension of the array A.  LDA >= max(1,N).
 73 *
 74 *  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
 75 *          On entry, the Hermitian matrix B.  If UPLO = 'U', the
 76 *          leading N-by-N upper triangular part of B contains the
 77 *          upper triangular part of the matrix B.  If UPLO = 'L',
 78 *          the leading N-by-N lower triangular part of B contains
 79 *          the lower triangular part of the matrix B.
 80 *
 81 *          On exit, if INFO <= N, the part of B containing the matrix is
 82 *          overwritten by the triangular factor U or L from the Cholesky
 83 *          factorization B = U**H*U or B = L*L**H.
 84 *
 85 *  LDB     (input) INTEGER
 86 *          The leading dimension of the array B.  LDB >= max(1,N).
 87 *
 88 *  W       (output) DOUBLE PRECISION array, dimension (N)
 89 *          If INFO = 0, the eigenvalues in ascending order.
 90 *
 91 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
 92 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 93 *
 94 *  LWORK   (input) INTEGER
 95 *          The length of the array WORK.
 96 *          If N <= 1,                LWORK >= 1.
 97 *          If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
 98 *          If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.
 99 *
100 *          If LWORK = -1, then a workspace query is assumed; the routine
101 *          only calculates the optimal sizes of the WORK, RWORK and
102 *          IWORK arrays, returns these values as the first entries of
103 *          the WORK, RWORK and IWORK arrays, and no error message
104 *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
105 *
106 *  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
107 *          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
108 *
109 *  LRWORK  (input) INTEGER
110 *          The dimension of the array RWORK.
111 *          If N <= 1,                LRWORK >= 1.
112 *          If JOBZ  = 'N' and N > 1, LRWORK >= N.
113 *          If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
114 *
115 *          If LRWORK = -1, then a workspace query is assumed; the
116 *          routine only calculates the optimal sizes of the WORK, RWORK
117 *          and IWORK arrays, returns these values as the first entries
118 *          of the WORK, RWORK and IWORK arrays, and no error message
119 *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
120 *
121 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
122 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
123 *
124 *  LIWORK  (input) INTEGER
125 *          The dimension of the array IWORK.
126 *          If N <= 1,                LIWORK >= 1.
127 *          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
128 *          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
129 *
130 *          If LIWORK = -1, then a workspace query is assumed; the
131 *          routine only calculates the optimal sizes of the WORK, RWORK
132 *          and IWORK arrays, returns these values as the first entries
133 *          of the WORK, RWORK and IWORK arrays, and no error message
134 *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
135 *
136 *  INFO    (output) INTEGER
137 *          = 0:  successful exit
138 *          < 0:  if INFO = -i, the i-th argument had an illegal value
139 *          > 0:  ZPOTRF or ZHEEVD returned an error code:
140 *             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
141 *                    failed to converge; i off-diagonal elements of an
142 *                    intermediate tridiagonal form did not converge to
143 *                    zero;
144 *                    if INFO = i and JOBZ = 'V', then the algorithm
145 *                    failed to compute an eigenvalue while working on
146 *                    the submatrix lying in rows and columns INFO/(N+1)
147 *                    through mod(INFO,N+1);
148 *             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
149 *                    minor of order i of B is not positive definite.
150 *                    The factorization of B could not be completed and
151 *                    no eigenvalues or eigenvectors were computed.
152 *
153 *  Further Details
154 *  ===============
155 *
156 *  Based on contributions by
157 *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
158 *
159 *  Modified so that no backsubstitution is performed if ZHEEVD fails to
160 *  converge (NEIG in old code could be greater than N causing out of
161 *  bounds reference to A - reported by Ralf Meyer).  Also corrected the
162 *  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
163 *  =====================================================================
164 *
165 *     .. Parameters ..
166       COMPLEX*16         CONE
167       PARAMETER          ( CONE = ( 1.0D+00.0D+0 ) )
168 *     ..
169 *     .. Local Scalars ..
170       LOGICAL            LQUERY, UPPER, WANTZ
171       CHARACTER          TRANS
172       INTEGER            LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
173 *     ..
174 *     .. External Functions ..
175       LOGICAL            LSAME
176       EXTERNAL           LSAME
177 *     ..
178 *     .. External Subroutines ..
179       EXTERNAL           XERBLA, ZHEEVD, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
180 *     ..
181 *     .. Intrinsic Functions ..
182       INTRINSIC          DBLEMAX
183 *     ..
184 *     .. Executable Statements ..
185 *
186 *     Test the input parameters.
187 *
188       WANTZ = LSAME( JOBZ, 'V' )
189       UPPER = LSAME( UPLO, 'U' )
190       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
191 *
192       INFO = 0
193       IF( N.LE.1 ) THEN
194          LWMIN = 1
195          LRWMIN = 1
196          LIWMIN = 1
197       ELSE IF( WANTZ ) THEN
198          LWMIN = 2*+ N*N
199          LRWMIN = 1 + 5*+ 2*N*N
200          LIWMIN = 3 + 5*N
201       ELSE
202          LWMIN = N + 1
203          LRWMIN = N
204          LIWMIN = 1
205       END IF
206       LOPT = LWMIN
207       LROPT = LRWMIN
208       LIOPT = LIWMIN
209       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
210          INFO = -1
211       ELSE IF.NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
212          INFO = -2
213       ELSE IF.NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
214          INFO = -3
215       ELSE IF( N.LT.0 ) THEN
216          INFO = -4
217       ELSE IF( LDA.LT.MAX1, N ) ) THEN
218          INFO = -6
219       ELSE IF( LDB.LT.MAX1, N ) ) THEN
220          INFO = -8
221       END IF
222 *
223       IF( INFO.EQ.0 ) THEN
224          WORK( 1 ) = LOPT
225          RWORK( 1 ) = LROPT
226          IWORK( 1 ) = LIOPT
227 *
228          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
229             INFO = -11
230          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
231             INFO = -13
232          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
233             INFO = -15
234          END IF
235       END IF
236 *
237       IF( INFO.NE.0 ) THEN
238          CALL XERBLA( 'ZHEGVD'-INFO )
239          RETURN
240       ELSE IF( LQUERY ) THEN
241          RETURN
242       END IF
243 *
244 *     Quick return if possible
245 *
246       IF( N.EQ.0 )
247      $   RETURN
248 *
249 *     Form a Cholesky factorization of B.
250 *
251       CALL ZPOTRF( UPLO, N, B, LDB, INFO )
252       IF( INFO.NE.0 ) THEN
253          INFO = N + INFO
254          RETURN
255       END IF
256 *
257 *     Transform problem to standard eigenvalue problem and solve.
258 *
259       CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
260       CALL ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK,
261      $             IWORK, LIWORK, INFO )
262       LOPT = MAXDBLE( LOPT ), DBLE( WORK( 1 ) ) )
263       LROPT = MAXDBLE( LROPT ), DBLE( RWORK( 1 ) ) )
264       LIOPT = MAXDBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
265 *
266       IF( WANTZ .AND. INFO.EQ.0 ) THEN
267 *
268 *        Backtransform eigenvectors to the original problem.
269 *
270          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
271 *
272 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
273 *           backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
274 *
275             IF( UPPER ) THEN
276                TRANS = 'N'
277             ELSE
278                TRANS = 'C'
279             END IF
280 *
281             CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
282      $                  B, LDB, A, LDA )
283 *
284          ELSE IF( ITYPE.EQ.3 ) THEN
285 *
286 *           For B*A*x=(lambda)*x;
287 *           backtransform eigenvectors: x = L*y or U**H *y
288 *
289             IF( UPPER ) THEN
290                TRANS = 'C'
291             ELSE
292                TRANS = 'N'
293             END IF
294 *
295             CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
296      $                  B, LDB, A, LDA )
297          END IF
298       END IF
299 *
300       WORK( 1 ) = LOPT
301       RWORK( 1 ) = LROPT
302       IWORK( 1 ) = LIOPT
303 *
304       RETURN
305 *
306 *     End of ZHEGVD
307 *
308       END