1*> \brief \b DSYGS2 reduces a symmetric definite generalized eigenproblem to standard form, using the factorization results obtained from spotrf (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 DSYGS2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygs2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygs2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygs2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSYGS2( 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*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> DSYGS2 reduces a real symmetric-definite generalized eigenproblem
38*> to standard form.
39*>
40*> If ITYPE = 1, the problem is A*x = lambda*B*x,
41*> and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
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**T or L**T *A*L.
45*>
46*> B must have been previously factorized as U**T *U or L*L**T by DPOTRF.
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] ITYPE
53*> \verbatim
54*>          ITYPE is INTEGER
55*>          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
56*>          = 2 or 3: compute U*A*U**T or L**T *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*>          symmetric 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 DOUBLE PRECISION array, dimension (LDA,N)
77*>          On entry, the symmetric 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] B
96*> \verbatim
97*>          B is DOUBLE PRECISION array, dimension (LDB,N)
98*>          The triangular factor from the Cholesky factorization of B,
99*>          as returned by DPOTRF.
100*> \endverbatim
101*>
102*> \param[in] LDB
103*> \verbatim
104*>          LDB is INTEGER
105*>          The leading dimension of the array B.  LDB >= max(1,N).
106*> \endverbatim
107*>
108*> \param[out] INFO
109*> \verbatim
110*>          INFO is INTEGER
111*>          = 0:  successful exit.
112*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
113*> \endverbatim
114*
115*  Authors:
116*  ========
117*
118*> \author Univ. of Tennessee
119*> \author Univ. of California Berkeley
120*> \author Univ. of Colorado Denver
121*> \author NAG Ltd.
122*
123*> \date September 2012
124*
125*> \ingroup doubleSYcomputational
126*
127*  =====================================================================
128      SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
129*
130*  -- LAPACK computational routine (version 3.4.2) --
131*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
132*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133*     September 2012
134*
135*     .. Scalar Arguments ..
136      CHARACTER          UPLO
137      INTEGER            INFO, ITYPE, LDA, LDB, N
138*     ..
139*     .. Array Arguments ..
140      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
141*     ..
142*
143*  =====================================================================
144*
145*     .. Parameters ..
146      DOUBLE PRECISION   ONE, HALF
147      PARAMETER          ( ONE = 1.0D0, HALF = 0.5D0 )
148*     ..
149*     .. Local Scalars ..
150      LOGICAL            UPPER
151      INTEGER            K
152      DOUBLE PRECISION   AKK, BKK, CT
153*     ..
154*     .. External Subroutines ..
155      EXTERNAL           DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA
156*     ..
157*     .. Intrinsic Functions ..
158      INTRINSIC          MAX
159*     ..
160*     .. External Functions ..
161      LOGICAL            LSAME
162      EXTERNAL           LSAME
163*     ..
164*     .. Executable Statements ..
165*
166*     Test the input parameters.
167*
168      INFO = 0
169      UPPER = LSAME( UPLO, 'U' )
170      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
171         INFO = -1
172      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
173         INFO = -2
174      ELSE IF( N.LT.0 ) THEN
175         INFO = -3
176      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
177         INFO = -5
178      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
179         INFO = -7
180      END IF
181      IF( INFO.NE.0 ) THEN
182         CALL XERBLA( 'DSYGS2', -INFO )
183         RETURN
184      END IF
185*
186      IF( ITYPE.EQ.1 ) THEN
187         IF( UPPER ) THEN
188*
189*           Compute inv(U**T)*A*inv(U)
190*
191            DO 10 K = 1, N
192*
193*              Update the upper triangle of A(k:n,k:n)
194*
195               AKK = A( K, K )
196               BKK = B( K, K )
197               AKK = AKK / BKK**2
198               A( K, K ) = AKK
199               IF( K.LT.N ) THEN
200                  CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
201                  CT = -HALF*AKK
202                  CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
203     $                        LDA )
204                  CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA,
205     $                        B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
206                  CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
207     $                        LDA )
208                  CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K,
209     $                        B( K+1, K+1 ), LDB, A( K, K+1 ), LDA )
210               END IF
211   10       CONTINUE
212         ELSE
213*
214*           Compute inv(L)*A*inv(L**T)
215*
216            DO 20 K = 1, N
217*
218*              Update the lower triangle of A(k:n,k:n)
219*
220               AKK = A( K, K )
221               BKK = B( K, K )
222               AKK = AKK / BKK**2
223               A( K, K ) = AKK
224               IF( K.LT.N ) THEN
225                  CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
226                  CT = -HALF*AKK
227                  CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
228                  CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1,
229     $                        B( K+1, K ), 1, A( K+1, K+1 ), LDA )
230                  CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
231                  CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
232     $                        B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
233               END IF
234   20       CONTINUE
235         END IF
236      ELSE
237         IF( UPPER ) THEN
238*
239*           Compute U*A*U**T
240*
241            DO 30 K = 1, N
242*
243*              Update the upper triangle of A(1:k,1:k)
244*
245               AKK = A( K, K )
246               BKK = B( K, K )
247               CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
248     $                     LDB, A( 1, K ), 1 )
249               CT = HALF*AKK
250               CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
251               CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1,
252     $                     A, LDA )
253               CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
254               CALL DSCAL( K-1, BKK, A( 1, K ), 1 )
255               A( K, K ) = AKK*BKK**2
256   30       CONTINUE
257         ELSE
258*
259*           Compute L**T *A*L
260*
261            DO 40 K = 1, N
262*
263*              Update the lower triangle of A(1:k,1:k)
264*
265               AKK = A( K, K )
266               BKK = B( K, K )
267               CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB,
268     $                     A( K, 1 ), LDA )
269               CT = HALF*AKK
270               CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
271               CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ),
272     $                     LDB, A, LDA )
273               CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
274               CALL DSCAL( K-1, BKK, A( K, 1 ), LDA )
275               A( K, K ) = AKK*BKK**2
276   40       CONTINUE
277         END IF
278      END IF
279      RETURN
280*
281*     End of DSYGS2
282*
283      END
284