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