001:       SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INFO, ITYPE, LDA, LDB, N
010: *     ..
011: *     .. Array Arguments ..
012:       COMPLEX*16         A( LDA, * ), B( LDB, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  ZHEGST reduces a complex Hermitian-definite generalized
019: *  eigenproblem to standard form.
020: *
021: *  If ITYPE = 1, the problem is A*x = lambda*B*x,
022: *  and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
023: *
024: *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
025: *  B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
026: *
027: *  B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  ITYPE   (input) INTEGER
033: *          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
034: *          = 2 or 3: compute U*A*U**H or L**H*A*L.
035: *
036: *  UPLO    (input) CHARACTER*1
037: *          = 'U':  Upper triangle of A is stored and B is factored as
038: *                  U**H*U;
039: *          = 'L':  Lower triangle of A is stored and B is factored as
040: *                  L*L**H.
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrices A and B.  N >= 0.
044: *
045: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
046: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
047: *          N-by-N upper triangular part of A contains the upper
048: *          triangular part of the matrix A, and the strictly lower
049: *          triangular part of A is not referenced.  If UPLO = 'L', the
050: *          leading N-by-N lower triangular part of A contains the lower
051: *          triangular part of the matrix A, and the strictly upper
052: *          triangular part of A is not referenced.
053: *
054: *          On exit, if INFO = 0, the transformed matrix, stored in the
055: *          same format as A.
056: *
057: *  LDA     (input) INTEGER
058: *          The leading dimension of the array A.  LDA >= max(1,N).
059: *
060: *  B       (input) COMPLEX*16 array, dimension (LDB,N)
061: *          The triangular factor from the Cholesky factorization of B,
062: *          as returned by ZPOTRF.
063: *
064: *  LDB     (input) INTEGER
065: *          The leading dimension of the array B.  LDB >= max(1,N).
066: *
067: *  INFO    (output) INTEGER
068: *          = 0:  successful exit
069: *          < 0:  if INFO = -i, the i-th argument had an illegal value
070: *
071: *  =====================================================================
072: *
073: *     .. Parameters ..
074:       DOUBLE PRECISION   ONE
075:       PARAMETER          ( ONE = 1.0D+0 )
076:       COMPLEX*16         CONE, HALF
077:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
078:      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
079: *     ..
080: *     .. Local Scalars ..
081:       LOGICAL            UPPER
082:       INTEGER            K, KB, NB
083: *     ..
084: *     .. External Subroutines ..
085:       EXTERNAL           XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
086: *     ..
087: *     .. Intrinsic Functions ..
088:       INTRINSIC          MAX, MIN
089: *     ..
090: *     .. External Functions ..
091:       LOGICAL            LSAME
092:       INTEGER            ILAENV
093:       EXTERNAL           LSAME, ILAENV
094: *     ..
095: *     .. Executable Statements ..
096: *
097: *     Test the input parameters.
098: *
099:       INFO = 0
100:       UPPER = LSAME( UPLO, 'U' )
101:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
102:          INFO = -1
103:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
104:          INFO = -2
105:       ELSE IF( N.LT.0 ) THEN
106:          INFO = -3
107:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
108:          INFO = -5
109:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
110:          INFO = -7
111:       END IF
112:       IF( INFO.NE.0 ) THEN
113:          CALL XERBLA( 'ZHEGST', -INFO )
114:          RETURN
115:       END IF
116: *
117: *     Quick return if possible
118: *
119:       IF( N.EQ.0 )
120:      $   RETURN
121: *
122: *     Determine the block size for this environment.
123: *
124:       NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
125: *
126:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
127: *
128: *        Use unblocked code
129: *
130:          CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
131:       ELSE
132: *
133: *        Use blocked code
134: *
135:          IF( ITYPE.EQ.1 ) THEN
136:             IF( UPPER ) THEN
137: *
138: *              Compute inv(U')*A*inv(U)
139: *
140:                DO 10 K = 1, N, NB
141:                   KB = MIN( N-K+1, NB )
142: *
143: *                 Update the upper triangle of A(k:n,k:n)
144: *
145:                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
146:      $                         B( K, K ), LDB, INFO )
147:                   IF( K+KB.LE.N ) THEN
148:                      CALL ZTRSM( 'Left', UPLO, 'Conjugate transpose',
149:      $                           'Non-unit', KB, N-K-KB+1, CONE,
150:      $                           B( K, K ), LDB, A( K, K+KB ), LDA )
151:                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
152:      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
153:      $                           CONE, A( K, K+KB ), LDA )
154:                      CALL ZHER2K( UPLO, 'Conjugate transpose', N-K-KB+1,
155:      $                            KB, -CONE, A( K, K+KB ), LDA,
156:      $                            B( K, K+KB ), LDB, ONE,
157:      $                            A( K+KB, K+KB ), LDA )
158:                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
159:      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
160:      $                           CONE, A( K, K+KB ), LDA )
161:                      CALL ZTRSM( 'Right', UPLO, 'No transpose',
162:      $                           'Non-unit', KB, N-K-KB+1, CONE,
163:      $                           B( K+KB, K+KB ), LDB, A( K, K+KB ),
164:      $                           LDA )
165:                   END IF
166:    10          CONTINUE
167:             ELSE
168: *
169: *              Compute inv(L)*A*inv(L')
170: *
171:                DO 20 K = 1, N, NB
172:                   KB = MIN( N-K+1, NB )
173: *
174: *                 Update the lower triangle of A(k:n,k:n)
175: *
176:                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
177:      $                         B( K, K ), LDB, INFO )
178:                   IF( K+KB.LE.N ) THEN
179:                      CALL ZTRSM( 'Right', UPLO, 'Conjugate transpose',
180:      $                           'Non-unit', N-K-KB+1, KB, CONE,
181:      $                           B( K, K ), LDB, A( K+KB, K ), LDA )
182:                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
183:      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
184:      $                           CONE, A( K+KB, K ), LDA )
185:                      CALL ZHER2K( UPLO, 'No transpose', N-K-KB+1, KB,
186:      $                            -CONE, A( K+KB, K ), LDA,
187:      $                            B( K+KB, K ), LDB, ONE,
188:      $                            A( K+KB, K+KB ), LDA )
189:                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
190:      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
191:      $                           CONE, A( K+KB, K ), LDA )
192:                      CALL ZTRSM( 'Left', UPLO, 'No transpose',
193:      $                           'Non-unit', N-K-KB+1, KB, CONE,
194:      $                           B( K+KB, K+KB ), LDB, A( K+KB, K ),
195:      $                           LDA )
196:                   END IF
197:    20          CONTINUE
198:             END IF
199:          ELSE
200:             IF( UPPER ) THEN
201: *
202: *              Compute U*A*U'
203: *
204:                DO 30 K = 1, N, NB
205:                   KB = MIN( N-K+1, NB )
206: *
207: *                 Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
208: *
209:                   CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
210:      $                        K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
211:                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
212:      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
213:      $                        LDA )
214:                   CALL ZHER2K( UPLO, 'No transpose', K-1, KB, CONE,
215:      $                         A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
216:      $                         LDA )
217:                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
218:      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
219:      $                        LDA )
220:                   CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose',
221:      $                        'Non-unit', K-1, KB, CONE, B( K, K ), LDB,
222:      $                        A( 1, K ), LDA )
223:                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
224:      $                         B( K, K ), LDB, INFO )
225:    30          CONTINUE
226:             ELSE
227: *
228: *              Compute L'*A*L
229: *
230:                DO 40 K = 1, N, NB
231:                   KB = MIN( N-K+1, NB )
232: *
233: *                 Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
234: *
235:                   CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
236:      $                        KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
237:                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
238:      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
239:      $                        LDA )
240:                   CALL ZHER2K( UPLO, 'Conjugate transpose', K-1, KB,
241:      $                         CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
242:      $                         ONE, A, LDA )
243:                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
244:      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
245:      $                        LDA )
246:                   CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose',
247:      $                        'Non-unit', KB, K-1, CONE, B( K, K ), LDB,
248:      $                        A( K, 1 ), LDA )
249:                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
250:      $                         B( K, K ), LDB, INFO )
251:    40          CONTINUE
252:             END IF
253:          END IF
254:       END IF
255:       RETURN
256: *
257: *     End of ZHEGST
258: *
259:       END
260: