1*> \brief \b CHEGS2 reduces a Hermitian definite generalized eigenproblem to standard form, using the factorization results obtained from cpotrf (unblocked algorithm).
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHEGS2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chegs2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chegs2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chegs2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, ITYPE, LDA, LDB, N
26*       ..
27*       .. Array Arguments ..
28*       COMPLEX            A( LDA, * ), B( LDB, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> CHEGS2 reduces a complex Hermitian-definite generalized
38*> eigenproblem to standard form.
39*>
40*> If ITYPE = 1, the problem is A*x = lambda*B*x,
41*> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
42*>
43*> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
44*> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H *A*L.
45*>
46*> B must have been previously factorized as U**H *U or L*L**H by ZPOTRF.
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] ITYPE
53*> \verbatim
54*>          ITYPE is INTEGER
55*>          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
56*>          = 2 or 3: compute U*A*U**H or L**H *A*L.
57*> \endverbatim
58*>
59*> \param[in] UPLO
60*> \verbatim
61*>          UPLO is CHARACTER*1
62*>          Specifies whether the upper or lower triangular part of the
63*>          Hermitian matrix A is stored, and how B has been factorized.
64*>          = 'U':  Upper triangular
65*>          = 'L':  Lower triangular
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*>          N is INTEGER
71*>          The order of the matrices A and B.  N >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] A
75*> \verbatim
76*>          A is COMPLEX array, dimension (LDA,N)
77*>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
78*>          n by n upper triangular part of A contains the upper
79*>          triangular part of the matrix A, and the strictly lower
80*>          triangular part of A is not referenced.  If UPLO = 'L', the
81*>          leading n by n lower triangular part of A contains the lower
82*>          triangular part of the matrix A, and the strictly upper
83*>          triangular part of A is not referenced.
84*>
85*>          On exit, if INFO = 0, the transformed matrix, stored in the
86*>          same format as A.
87*> \endverbatim
88*>
89*> \param[in] LDA
90*> \verbatim
91*>          LDA is INTEGER
92*>          The leading dimension of the array A.  LDA >= max(1,N).
93*> \endverbatim
94*>
95*> \param[in,out] B
96*> \verbatim
97*>          B is COMPLEX array, dimension (LDB,N)
98*>          The triangular factor from the Cholesky factorization of B,
99*>          as returned by CPOTRF.
100*>          B is modified by the routine but restored on exit.
101*> \endverbatim
102*>
103*> \param[in] LDB
104*> \verbatim
105*>          LDB is INTEGER
106*>          The leading dimension of the array B.  LDB >= max(1,N).
107*> \endverbatim
108*>
109*> \param[out] INFO
110*> \verbatim
111*>          INFO is INTEGER
112*>          = 0:  successful exit.
113*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
114*> \endverbatim
115*
116*  Authors:
117*  ========
118*
119*> \author Univ. of Tennessee
120*> \author Univ. of California Berkeley
121*> \author Univ. of Colorado Denver
122*> \author NAG Ltd.
123*
124*> \date December 2016
125*
126*> \ingroup complexHEcomputational
127*
128*  =====================================================================
129      SUBROUTINE CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
130*
131*  -- LAPACK computational routine (version 3.7.0) --
132*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
133*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134*     December 2016
135*
136*     .. Scalar Arguments ..
137      CHARACTER          UPLO
138      INTEGER            INFO, ITYPE, LDA, LDB, N
139*     ..
140*     .. Array Arguments ..
141      COMPLEX            A( LDA, * ), B( LDB, * )
142*     ..
143*
144*  =====================================================================
145*
146*     .. Parameters ..
147      REAL               ONE, HALF
148      PARAMETER          ( ONE = 1.0E+0, HALF = 0.5E+0 )
149      COMPLEX            CONE
150      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
151*     ..
152*     .. Local Scalars ..
153      LOGICAL            UPPER
154      INTEGER            K
155      REAL               AKK, BKK
156      COMPLEX            CT
157*     ..
158*     .. External Subroutines ..
159      EXTERNAL           CAXPY, CHER2, CLACGV, CSSCAL, CTRMV, CTRSV,
160     $                   XERBLA
161*     ..
162*     .. Intrinsic Functions ..
163      INTRINSIC          MAX
164*     ..
165*     .. External Functions ..
166      LOGICAL            LSAME
167      EXTERNAL           LSAME
168*     ..
169*     .. Executable Statements ..
170*
171*     Test the input parameters.
172*
173      INFO = 0
174      UPPER = LSAME( UPLO, 'U' )
175      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
176         INFO = -1
177      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
178         INFO = -2
179      ELSE IF( N.LT.0 ) THEN
180         INFO = -3
181      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
182         INFO = -5
183      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
184         INFO = -7
185      END IF
186      IF( INFO.NE.0 ) THEN
187         CALL XERBLA( 'CHEGS2', -INFO )
188         RETURN
189      END IF
190*
191      IF( ITYPE.EQ.1 ) THEN
192         IF( UPPER ) THEN
193*
194*           Compute inv(U**H)*A*inv(U)
195*
196            DO 10 K = 1, N
197*
198*              Update the upper triangle of A(k:n,k:n)
199*
200               AKK = A( K, K )
201               BKK = B( K, K )
202               AKK = AKK / BKK**2
203               A( K, K ) = AKK
204               IF( K.LT.N ) THEN
205                  CALL CSSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
206                  CT = -HALF*AKK
207                  CALL CLACGV( N-K, A( K, K+1 ), LDA )
208                  CALL CLACGV( N-K, B( K, K+1 ), LDB )
209                  CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
210     $                        LDA )
211                  CALL CHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
212     $                        B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
213                  CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
214     $                        LDA )
215                  CALL CLACGV( N-K, B( K, K+1 ), LDB )
216                  CALL CTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
217     $                        N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
218     $                        LDA )
219                  CALL CLACGV( N-K, A( K, K+1 ), LDA )
220               END IF
221   10       CONTINUE
222         ELSE
223*
224*           Compute inv(L)*A*inv(L**H)
225*
226            DO 20 K = 1, N
227*
228*              Update the lower triangle of A(k:n,k:n)
229*
230               AKK = A( K, K )
231               BKK = B( K, K )
232               AKK = AKK / BKK**2
233               A( K, K ) = AKK
234               IF( K.LT.N ) THEN
235                  CALL CSSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
236                  CT = -HALF*AKK
237                  CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
238                  CALL CHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
239     $                        B( K+1, K ), 1, A( K+1, K+1 ), LDA )
240                  CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
241                  CALL CTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
242     $                        B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
243               END IF
244   20       CONTINUE
245         END IF
246      ELSE
247         IF( UPPER ) THEN
248*
249*           Compute U*A*U**H
250*
251            DO 30 K = 1, N
252*
253*              Update the upper triangle of A(1:k,1:k)
254*
255               AKK = A( K, K )
256               BKK = B( K, K )
257               CALL CTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
258     $                     LDB, A( 1, K ), 1 )
259               CT = HALF*AKK
260               CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
261               CALL CHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
262     $                     A, LDA )
263               CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
264               CALL CSSCAL( K-1, BKK, A( 1, K ), 1 )
265               A( K, K ) = AKK*BKK**2
266   30       CONTINUE
267         ELSE
268*
269*           Compute L**H *A*L
270*
271            DO 40 K = 1, N
272*
273*              Update the lower triangle of A(1:k,1:k)
274*
275               AKK = A( K, K )
276               BKK = B( K, K )
277               CALL CLACGV( K-1, A( K, 1 ), LDA )
278               CALL CTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1,
279     $                     B, LDB, A( K, 1 ), LDA )
280               CT = HALF*AKK
281               CALL CLACGV( K-1, B( K, 1 ), LDB )
282               CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
283               CALL CHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
284     $                     LDB, A, LDA )
285               CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
286               CALL CLACGV( K-1, B( K, 1 ), LDB )
287               CALL CSSCAL( K-1, BKK, A( K, 1 ), LDA )
288               CALL CLACGV( K-1, A( K, 1 ), LDA )
289               A( K, K ) = AKK*BKK**2
290   40       CONTINUE
291         END IF
292      END IF
293      RETURN
294*
295*     End of CHEGS2
296*
297      END
298