1*> \brief \b DSYTRS 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSYTRS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrs.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrs.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrs.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, LDA, LDB, N, NRHS 26* .. 27* .. Array Arguments .. 28* INTEGER IPIV( * ) 29* DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DSYTRS solves a system of linear equations A*X = B with a real 39*> symmetric matrix A using the factorization A = U*D*U**T or 40*> A = L*D*L**T computed by DSYTRF. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] UPLO 47*> \verbatim 48*> UPLO is CHARACTER*1 49*> Specifies whether the details of the factorization are stored 50*> as an upper or lower triangular matrix. 51*> = 'U': Upper triangular, form is A = U*D*U**T; 52*> = 'L': Lower triangular, form is A = L*D*L**T. 53*> \endverbatim 54*> 55*> \param[in] N 56*> \verbatim 57*> N is INTEGER 58*> The order of the matrix A. N >= 0. 59*> \endverbatim 60*> 61*> \param[in] NRHS 62*> \verbatim 63*> NRHS is INTEGER 64*> The number of right hand sides, i.e., the number of columns 65*> of the matrix B. NRHS >= 0. 66*> \endverbatim 67*> 68*> \param[in] A 69*> \verbatim 70*> A is DOUBLE PRECISION array, dimension (LDA,N) 71*> The block diagonal matrix D and the multipliers used to 72*> obtain the factor U or L as computed by DSYTRF. 73*> \endverbatim 74*> 75*> \param[in] LDA 76*> \verbatim 77*> LDA is INTEGER 78*> The leading dimension of the array A. LDA >= max(1,N). 79*> \endverbatim 80*> 81*> \param[in] IPIV 82*> \verbatim 83*> IPIV is INTEGER array, dimension (N) 84*> Details of the interchanges and the block structure of D 85*> as determined by DSYTRF. 86*> \endverbatim 87*> 88*> \param[in,out] B 89*> \verbatim 90*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 91*> On entry, the right hand side matrix B. 92*> On exit, the solution matrix X. 93*> \endverbatim 94*> 95*> \param[in] LDB 96*> \verbatim 97*> LDB is INTEGER 98*> The leading dimension of the array B. LDB >= max(1,N). 99*> \endverbatim 100*> 101*> \param[out] INFO 102*> \verbatim 103*> INFO is INTEGER 104*> = 0: successful exit 105*> < 0: if INFO = -i, the i-th argument had an illegal value 106*> \endverbatim 107* 108* Authors: 109* ======== 110* 111*> \author Univ. of Tennessee 112*> \author Univ. of California Berkeley 113*> \author Univ. of Colorado Denver 114*> \author NAG Ltd. 115* 116*> \date November 2011 117* 118*> \ingroup doubleSYcomputational 119* 120* ===================================================================== 121 SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) 122* 123* -- LAPACK computational routine (version 3.4.0) -- 124* -- LAPACK is a software package provided by Univ. of Tennessee, -- 125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 126* November 2011 127* 128* .. Scalar Arguments .. 129 CHARACTER UPLO 130 INTEGER INFO, LDA, LDB, N, NRHS 131* .. 132* .. Array Arguments .. 133 INTEGER IPIV( * ) 134 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 135* .. 136* 137* ===================================================================== 138* 139* .. Parameters .. 140 DOUBLE PRECISION ONE 141 PARAMETER ( ONE = 1.0D+0 ) 142* .. 143* .. Local Scalars .. 144 LOGICAL UPPER 145 INTEGER J, K, KP 146 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM 147* .. 148* .. External Functions .. 149 LOGICAL LSAME 150 EXTERNAL LSAME 151* .. 152* .. External Subroutines .. 153 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA 154* .. 155* .. Intrinsic Functions .. 156 INTRINSIC MAX 157* .. 158* .. Executable Statements .. 159* 160 INFO = 0 161 UPPER = LSAME( UPLO, 'U' ) 162 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 163 INFO = -1 164 ELSE IF( N.LT.0 ) THEN 165 INFO = -2 166 ELSE IF( NRHS.LT.0 ) THEN 167 INFO = -3 168 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 169 INFO = -5 170 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 171 INFO = -8 172 END IF 173 IF( INFO.NE.0 ) THEN 174 CALL XERBLA( 'DSYTRS', -INFO ) 175 RETURN 176 END IF 177* 178* Quick return if possible 179* 180 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 181 $ RETURN 182* 183 IF( UPPER ) THEN 184* 185* Solve A*X = B, where A = U*D*U**T. 186* 187* First solve U*D*X = B, overwriting B with X. 188* 189* K is the main loop index, decreasing from N to 1 in steps of 190* 1 or 2, depending on the size of the diagonal blocks. 191* 192 K = N 193 10 CONTINUE 194* 195* If K < 1, exit from loop. 196* 197 IF( K.LT.1 ) 198 $ GO TO 30 199* 200 IF( IPIV( K ).GT.0 ) THEN 201* 202* 1 x 1 diagonal block 203* 204* Interchange rows K and IPIV(K). 205* 206 KP = IPIV( K ) 207 IF( KP.NE.K ) 208 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 209* 210* Multiply by inv(U(K)), where U(K) is the transformation 211* stored in column K of A. 212* 213 CALL DGER( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, 214 $ B( 1, 1 ), LDB ) 215* 216* Multiply by the inverse of the diagonal block. 217* 218 CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) 219 K = K - 1 220 ELSE 221* 222* 2 x 2 diagonal block 223* 224* Interchange rows K-1 and -IPIV(K). 225* 226 KP = -IPIV( K ) 227 IF( KP.NE.K-1 ) 228 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 229* 230* Multiply by inv(U(K)), where U(K) is the transformation 231* stored in columns K-1 and K of A. 232* 233 CALL DGER( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, 234 $ B( 1, 1 ), LDB ) 235 CALL DGER( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ), 236 $ LDB, B( 1, 1 ), LDB ) 237* 238* Multiply by the inverse of the diagonal block. 239* 240 AKM1K = A( K-1, K ) 241 AKM1 = A( K-1, K-1 ) / AKM1K 242 AK = A( K, K ) / AKM1K 243 DENOM = AKM1*AK - ONE 244 DO 20 J = 1, NRHS 245 BKM1 = B( K-1, J ) / AKM1K 246 BK = B( K, J ) / AKM1K 247 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM 248 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM 249 20 CONTINUE 250 K = K - 2 251 END IF 252* 253 GO TO 10 254 30 CONTINUE 255* 256* Next solve U**T *X = B, overwriting B with X. 257* 258* K is the main loop index, increasing from 1 to N in steps of 259* 1 or 2, depending on the size of the diagonal blocks. 260* 261 K = 1 262 40 CONTINUE 263* 264* If K > N, exit from loop. 265* 266 IF( K.GT.N ) 267 $ GO TO 50 268* 269 IF( IPIV( K ).GT.0 ) THEN 270* 271* 1 x 1 diagonal block 272* 273* Multiply by inv(U**T(K)), where U(K) is the transformation 274* stored in column K of A. 275* 276 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), 277 $ 1, ONE, B( K, 1 ), LDB ) 278* 279* Interchange rows K and IPIV(K). 280* 281 KP = IPIV( K ) 282 IF( KP.NE.K ) 283 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 284 K = K + 1 285 ELSE 286* 287* 2 x 2 diagonal block 288* 289* Multiply by inv(U**T(K+1)), where U(K+1) is the transformation 290* stored in columns K and K+1 of A. 291* 292 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), 293 $ 1, ONE, B( K, 1 ), LDB ) 294 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, 295 $ A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB ) 296* 297* Interchange rows K and -IPIV(K). 298* 299 KP = -IPIV( K ) 300 IF( KP.NE.K ) 301 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 302 K = K + 2 303 END IF 304* 305 GO TO 40 306 50 CONTINUE 307* 308 ELSE 309* 310* Solve A*X = B, where A = L*D*L**T. 311* 312* First solve L*D*X = B, overwriting B with X. 313* 314* K is the main loop index, increasing from 1 to N in steps of 315* 1 or 2, depending on the size of the diagonal blocks. 316* 317 K = 1 318 60 CONTINUE 319* 320* If K > N, exit from loop. 321* 322 IF( K.GT.N ) 323 $ GO TO 80 324* 325 IF( IPIV( K ).GT.0 ) THEN 326* 327* 1 x 1 diagonal block 328* 329* Interchange rows K and IPIV(K). 330* 331 KP = IPIV( K ) 332 IF( KP.NE.K ) 333 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 334* 335* Multiply by inv(L(K)), where L(K) is the transformation 336* stored in column K of A. 337* 338 IF( K.LT.N ) 339 $ CALL DGER( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ), 340 $ LDB, B( K+1, 1 ), LDB ) 341* 342* Multiply by the inverse of the diagonal block. 343* 344 CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) 345 K = K + 1 346 ELSE 347* 348* 2 x 2 diagonal block 349* 350* Interchange rows K+1 and -IPIV(K). 351* 352 KP = -IPIV( K ) 353 IF( KP.NE.K+1 ) 354 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 355* 356* Multiply by inv(L(K)), where L(K) is the transformation 357* stored in columns K and K+1 of A. 358* 359 IF( K.LT.N-1 ) THEN 360 CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ), 361 $ LDB, B( K+2, 1 ), LDB ) 362 CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1, 363 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) 364 END IF 365* 366* Multiply by the inverse of the diagonal block. 367* 368 AKM1K = A( K+1, K ) 369 AKM1 = A( K, K ) / AKM1K 370 AK = A( K+1, K+1 ) / AKM1K 371 DENOM = AKM1*AK - ONE 372 DO 70 J = 1, NRHS 373 BKM1 = B( K, J ) / AKM1K 374 BK = B( K+1, J ) / AKM1K 375 B( K, J ) = ( AK*BKM1-BK ) / DENOM 376 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 377 70 CONTINUE 378 K = K + 2 379 END IF 380* 381 GO TO 60 382 80 CONTINUE 383* 384* Next solve L**T *X = B, overwriting B with X. 385* 386* K is the main loop index, decreasing from N to 1 in steps of 387* 1 or 2, depending on the size of the diagonal blocks. 388* 389 K = N 390 90 CONTINUE 391* 392* If K < 1, exit from loop. 393* 394 IF( K.LT.1 ) 395 $ GO TO 100 396* 397 IF( IPIV( K ).GT.0 ) THEN 398* 399* 1 x 1 diagonal block 400* 401* Multiply by inv(L**T(K)), where L(K) is the transformation 402* stored in column K of A. 403* 404 IF( K.LT.N ) 405 $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 406 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 407* 408* Interchange rows K and IPIV(K). 409* 410 KP = IPIV( K ) 411 IF( KP.NE.K ) 412 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 413 K = K - 1 414 ELSE 415* 416* 2 x 2 diagonal block 417* 418* Multiply by inv(L**T(K-1)), where L(K-1) is the transformation 419* stored in columns K-1 and K of A. 420* 421 IF( K.LT.N ) THEN 422 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 423 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 424 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 425 $ LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ), 426 $ LDB ) 427 END IF 428* 429* Interchange rows K and -IPIV(K). 430* 431 KP = -IPIV( K ) 432 IF( KP.NE.K ) 433 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 434 K = K - 2 435 END IF 436* 437 GO TO 90 438 100 CONTINUE 439 END IF 440* 441 RETURN 442* 443* End of DSYTRS 444* 445 END 446