1*> \brief \b SSYGST
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSYGST + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssygst.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssygst.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssygst.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SSYGST( 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*       REAL               A( LDA, * ), B( LDB, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> SSYGST 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 SPOTRF.
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*>          = 'U':  Upper triangle of A is stored and B is factored as
63*>                  U**T*U;
64*>          = 'L':  Lower triangle of A is stored and B is factored as
65*>                  L*L**T.
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 REAL 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 REAL array, dimension (LDB,N)
98*>          The triangular factor from the Cholesky factorization of B,
99*>          as returned by SPOTRF.
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*> \ingroup realSYcomputational
124*
125*  =====================================================================
126      SUBROUTINE SSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
127*
128*  -- LAPACK computational routine --
129*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
130*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131*
132*     .. Scalar Arguments ..
133      CHARACTER          UPLO
134      INTEGER            INFO, ITYPE, LDA, LDB, N
135*     ..
136*     .. Array Arguments ..
137      REAL               A( LDA, * ), B( LDB, * )
138*     ..
139*
140*  =====================================================================
141*
142*     .. Parameters ..
143      REAL               ONE, HALF
144      PARAMETER          ( ONE = 1.0, HALF = 0.5 )
145*     ..
146*     .. Local Scalars ..
147      LOGICAL            UPPER
148      INTEGER            K, KB, NB
149*     ..
150*     .. External Subroutines ..
151      EXTERNAL           SSYGS2, SSYMM, SSYR2K, STRMM, STRSM, XERBLA
152*     ..
153*     .. Intrinsic Functions ..
154      INTRINSIC          MAX, MIN
155*     ..
156*     .. External Functions ..
157      LOGICAL            LSAME
158      INTEGER            ILAENV
159      EXTERNAL           LSAME, ILAENV
160*     ..
161*     .. Executable Statements ..
162*
163*     Test the input parameters.
164*
165      INFO = 0
166      UPPER = LSAME( UPLO, 'U' )
167      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
168         INFO = -1
169      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
170         INFO = -2
171      ELSE IF( N.LT.0 ) THEN
172         INFO = -3
173      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
174         INFO = -5
175      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
176         INFO = -7
177      END IF
178      IF( INFO.NE.0 ) THEN
179         CALL XERBLA( 'SSYGST', -INFO )
180         RETURN
181      END IF
182*
183*     Quick return if possible
184*
185      IF( N.EQ.0 )
186     $   RETURN
187*
188*     Determine the block size for this environment.
189*
190      NB = ILAENV( 1, 'SSYGST', UPLO, N, -1, -1, -1 )
191*
192      IF( NB.LE.1 .OR. NB.GE.N ) THEN
193*
194*        Use unblocked code
195*
196         CALL SSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
197      ELSE
198*
199*        Use blocked code
200*
201         IF( ITYPE.EQ.1 ) THEN
202            IF( UPPER ) THEN
203*
204*              Compute inv(U**T)*A*inv(U)
205*
206               DO 10 K = 1, N, NB
207                  KB = MIN( N-K+1, NB )
208*
209*                 Update the upper triangle of A(k:n,k:n)
210*
211                  CALL SSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
212     $                         B( K, K ), LDB, INFO )
213                  IF( K+KB.LE.N ) THEN
214                     CALL STRSM( 'Left', UPLO, 'Transpose', 'Non-unit',
215     $                           KB, N-K-KB+1, ONE, B( K, K ), LDB,
216     $                           A( K, K+KB ), LDA )
217                     CALL SSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
218     $                           A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
219     $                           A( K, K+KB ), LDA )
220                     CALL SSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE,
221     $                            A( K, K+KB ), LDA, B( K, K+KB ), LDB,
222     $                            ONE, A( K+KB, K+KB ), LDA )
223                     CALL SSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
224     $                           A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
225     $                           A( K, K+KB ), LDA )
226                     CALL STRSM( 'Right', UPLO, 'No transpose',
227     $                           'Non-unit', KB, N-K-KB+1, ONE,
228     $                           B( K+KB, K+KB ), LDB, A( K, K+KB ),
229     $                           LDA )
230                  END IF
231   10          CONTINUE
232            ELSE
233*
234*              Compute inv(L)*A*inv(L**T)
235*
236               DO 20 K = 1, N, NB
237                  KB = MIN( N-K+1, NB )
238*
239*                 Update the lower triangle of A(k:n,k:n)
240*
241                  CALL SSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
242     $                         B( K, K ), LDB, INFO )
243                  IF( K+KB.LE.N ) THEN
244                     CALL STRSM( 'Right', UPLO, 'Transpose', 'Non-unit',
245     $                           N-K-KB+1, KB, ONE, B( K, K ), LDB,
246     $                           A( K+KB, K ), LDA )
247                     CALL SSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
248     $                           A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
249     $                           A( K+KB, K ), LDA )
250                     CALL SSYR2K( UPLO, 'No transpose', N-K-KB+1, KB,
251     $                            -ONE, A( K+KB, K ), LDA, B( K+KB, K ),
252     $                            LDB, ONE, A( K+KB, K+KB ), LDA )
253                     CALL SSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
254     $                           A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
255     $                           A( K+KB, K ), LDA )
256                     CALL STRSM( 'Left', UPLO, 'No transpose',
257     $                           'Non-unit', N-K-KB+1, KB, ONE,
258     $                           B( K+KB, K+KB ), LDB, A( K+KB, K ),
259     $                           LDA )
260                  END IF
261   20          CONTINUE
262            END IF
263         ELSE
264            IF( UPPER ) THEN
265*
266*              Compute U*A*U**T
267*
268               DO 30 K = 1, N, NB
269                  KB = MIN( N-K+1, NB )
270*
271*                 Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
272*
273                  CALL STRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
274     $                        K-1, KB, ONE, B, LDB, A( 1, K ), LDA )
275                  CALL SSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
276     $                        LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
277                  CALL SSYR2K( UPLO, 'No transpose', K-1, KB, ONE,
278     $                         A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
279     $                         LDA )
280                  CALL SSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
281     $                        LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
282                  CALL STRMM( 'Right', UPLO, 'Transpose', 'Non-unit',
283     $                        K-1, KB, ONE, B( K, K ), LDB, A( 1, K ),
284     $                        LDA )
285                  CALL SSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
286     $                         B( K, K ), LDB, INFO )
287   30          CONTINUE
288            ELSE
289*
290*              Compute L**T*A*L
291*
292               DO 40 K = 1, N, NB
293                  KB = MIN( N-K+1, NB )
294*
295*                 Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
296*
297                  CALL STRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
298     $                        KB, K-1, ONE, B, LDB, A( K, 1 ), LDA )
299                  CALL SSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
300     $                        LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
301                  CALL SSYR2K( UPLO, 'Transpose', K-1, KB, ONE,
302     $                         A( K, 1 ), LDA, B( K, 1 ), LDB, ONE, A,
303     $                         LDA )
304                  CALL SSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
305     $                        LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
306                  CALL STRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB,
307     $                        K-1, ONE, B( K, K ), LDB, A( K, 1 ), LDA )
308                  CALL SSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
309     $                         B( K, K ), LDB, INFO )
310   40          CONTINUE
311            END IF
312         END IF
313      END IF
314      RETURN
315*
316*     End of SSYGST
317*
318      END
319