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