1*> \brief \b ZHPGST
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHPGST + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgst.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgst.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgst.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, ITYPE, N
26*       ..
27*       .. Array Arguments ..
28*       COMPLEX*16         AP( * ), BP( * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> ZHPGST reduces a complex Hermitian-definite generalized
38*> eigenproblem to standard form, using packed storage.
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 ZPPTRF.
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*>          = 'U':  Upper triangle of A is stored and B is factored as
63*>                  U**H*U;
64*>          = 'L':  Lower triangle of A is stored and B is factored as
65*>                  L*L**H.
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] AP
75*> \verbatim
76*>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
77*>          On entry, the upper or lower triangle of the Hermitian matrix
78*>          A, packed columnwise in a linear array.  The j-th column of A
79*>          is stored in the array AP as follows:
80*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
81*>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
82*>
83*>          On exit, if INFO = 0, the transformed matrix, stored in the
84*>          same format as A.
85*> \endverbatim
86*>
87*> \param[in] BP
88*> \verbatim
89*>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
90*>          The triangular factor from the Cholesky factorization of B,
91*>          stored in the same format as A, as returned by ZPPTRF.
92*> \endverbatim
93*>
94*> \param[out] INFO
95*> \verbatim
96*>          INFO is INTEGER
97*>          = 0:  successful exit
98*>          < 0:  if INFO = -i, the i-th argument had an illegal value
99*> \endverbatim
100*
101*  Authors:
102*  ========
103*
104*> \author Univ. of Tennessee
105*> \author Univ. of California Berkeley
106*> \author Univ. of Colorado Denver
107*> \author NAG Ltd.
108*
109*> \date December 2016
110*
111*> \ingroup complex16OTHERcomputational
112*
113*  =====================================================================
114      SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
115*
116*  -- LAPACK computational routine (version 3.7.0) --
117*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
118*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
119*     December 2016
120*
121*     .. Scalar Arguments ..
122      CHARACTER          UPLO
123      INTEGER            INFO, ITYPE, N
124*     ..
125*     .. Array Arguments ..
126      COMPLEX*16         AP( * ), BP( * )
127*     ..
128*
129*  =====================================================================
130*
131*     .. Parameters ..
132      DOUBLE PRECISION   ONE, HALF
133      PARAMETER          ( ONE = 1.0D+0, HALF = 0.5D+0 )
134      COMPLEX*16         CONE
135      PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
136*     ..
137*     .. Local Scalars ..
138      LOGICAL            UPPER
139      INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
140      DOUBLE PRECISION   AJJ, AKK, BJJ, BKK
141      COMPLEX*16         CT
142*     ..
143*     .. External Subroutines ..
144      EXTERNAL           XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV,
145     $                   ZTPSV
146*     ..
147*     .. Intrinsic Functions ..
148      INTRINSIC          DBLE
149*     ..
150*     .. External Functions ..
151      LOGICAL            LSAME
152      COMPLEX*16         ZDOTC
153      EXTERNAL           LSAME, ZDOTC
154*     ..
155*     .. Executable Statements ..
156*
157*     Test the input parameters.
158*
159      INFO = 0
160      UPPER = LSAME( UPLO, 'U' )
161      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
162         INFO = -1
163      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
164         INFO = -2
165      ELSE IF( N.LT.0 ) THEN
166         INFO = -3
167      END IF
168      IF( INFO.NE.0 ) THEN
169         CALL XERBLA( 'ZHPGST', -INFO )
170         RETURN
171      END IF
172*
173      IF( ITYPE.EQ.1 ) THEN
174         IF( UPPER ) THEN
175*
176*           Compute inv(U**H)*A*inv(U)
177*
178*           J1 and JJ are the indices of A(1,j) and A(j,j)
179*
180            JJ = 0
181            DO 10 J = 1, N
182               J1 = JJ + 1
183               JJ = JJ + J
184*
185*              Compute the j-th column of the upper triangle of A
186*
187               AP( JJ ) = DBLE( AP( JJ ) )
188               BJJ = BP( JJ )
189               CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J,
190     $                     BP, AP( J1 ), 1 )
191               CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE,
192     $                     AP( J1 ), 1 )
193               CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
194               AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ),
195     $                    1 ) ) / BJJ
196   10       CONTINUE
197         ELSE
198*
199*           Compute inv(L)*A*inv(L**H)
200*
201*           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
202*
203            KK = 1
204            DO 20 K = 1, N
205               K1K1 = KK + N - K + 1
206*
207*              Update the lower triangle of A(k:n,k:n)
208*
209               AKK = AP( KK )
210               BKK = BP( KK )
211               AKK = AKK / BKK**2
212               AP( KK ) = AKK
213               IF( K.LT.N ) THEN
214                  CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
215                  CT = -HALF*AKK
216                  CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
217                  CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1,
218     $                        BP( KK+1 ), 1, AP( K1K1 ) )
219                  CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
220                  CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
221     $                        BP( K1K1 ), AP( KK+1 ), 1 )
222               END IF
223               KK = K1K1
224   20       CONTINUE
225         END IF
226      ELSE
227         IF( UPPER ) THEN
228*
229*           Compute U*A*U**H
230*
231*           K1 and KK are the indices of A(1,k) and A(k,k)
232*
233            KK = 0
234            DO 30 K = 1, N
235               K1 = KK + 1
236               KK = KK + K
237*
238*              Update the upper triangle of A(1:k,1:k)
239*
240               AKK = AP( KK )
241               BKK = BP( KK )
242               CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
243     $                     AP( K1 ), 1 )
244               CT = HALF*AKK
245               CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
246               CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1,
247     $                     AP )
248               CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
249               CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 )
250               AP( KK ) = AKK*BKK**2
251   30       CONTINUE
252         ELSE
253*
254*           Compute L**H *A*L
255*
256*           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
257*
258            JJ = 1
259            DO 40 J = 1, N
260               J1J1 = JJ + N - J + 1
261*
262*              Compute the j-th column of the lower triangle of A
263*
264               AJJ = AP( JJ )
265               BJJ = BP( JJ )
266               AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1,
267     $                    BP( JJ+1 ), 1 )
268               CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
269               CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1,
270     $                     CONE, AP( JJ+1 ), 1 )
271               CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit',
272     $                     N-J+1, BP( JJ ), AP( JJ ), 1 )
273               JJ = J1J1
274   40       CONTINUE
275         END IF
276      END IF
277      RETURN
278*
279*     End of ZHPGST
280*
281      END
282