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