1*> \brief \b ZSYTRS_ROOK 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZSYTRS_ROOK + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrs_rook.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrs_rook.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrs_rook.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZSYTRS_ROOK( 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* COMPLEX*16 A( LDA, * ), B( LDB, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZSYTRS_ROOK solves a system of linear equations A*X = B with 39*> a complex symmetric matrix A using the factorization A = U*D*U**T or 40*> A = L*D*L**T computed by ZSYTRF_ROOK. 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 COMPLEX*16 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 ZSYTRF_ROOK. 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 ZSYTRF_ROOK. 86*> \endverbatim 87*> 88*> \param[in,out] B 89*> \verbatim 90*> B is COMPLEX*16 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*> \ingroup complex16SYcomputational 117* 118*> \par Contributors: 119* ================== 120*> 121*> \verbatim 122*> 123*> December 2016, Igor Kozachenko, 124*> Computer Science Division, 125*> University of California, Berkeley 126*> 127*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 128*> School of Mathematics, 129*> University of Manchester 130*> 131*> \endverbatim 132* 133* ===================================================================== 134 SUBROUTINE ZSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 135 $ INFO ) 136* 137* -- LAPACK computational routine -- 138* -- LAPACK is a software package provided by Univ. of Tennessee, -- 139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 140* 141* .. Scalar Arguments .. 142 CHARACTER UPLO 143 INTEGER INFO, LDA, LDB, N, NRHS 144* .. 145* .. Array Arguments .. 146 INTEGER IPIV( * ) 147 COMPLEX*16 A( LDA, * ), B( LDB, * ) 148* .. 149* 150* ===================================================================== 151* 152* .. Parameters .. 153 COMPLEX*16 CONE 154 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 155* .. 156* .. Local Scalars .. 157 LOGICAL UPPER 158 INTEGER J, K, KP 159 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM 160* .. 161* .. External Functions .. 162 LOGICAL LSAME 163 EXTERNAL LSAME 164* .. 165* .. External Subroutines .. 166 EXTERNAL ZGEMV, ZGERU, ZSCAL, ZSWAP, XERBLA 167* .. 168* .. Intrinsic Functions .. 169 INTRINSIC MAX 170* .. 171* .. Executable Statements .. 172* 173 INFO = 0 174 UPPER = LSAME( UPLO, 'U' ) 175 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 176 INFO = -1 177 ELSE IF( N.LT.0 ) THEN 178 INFO = -2 179 ELSE IF( NRHS.LT.0 ) THEN 180 INFO = -3 181 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 182 INFO = -5 183 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 184 INFO = -8 185 END IF 186 IF( INFO.NE.0 ) THEN 187 CALL XERBLA( 'ZSYTRS_ROOK', -INFO ) 188 RETURN 189 END IF 190* 191* Quick return if possible 192* 193 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 194 $ RETURN 195* 196 IF( UPPER ) THEN 197* 198* Solve A*X = B, where A = U*D*U**T. 199* 200* First solve U*D*X = B, overwriting B with X. 201* 202* K is the main loop index, decreasing from N to 1 in steps of 203* 1 or 2, depending on the size of the diagonal blocks. 204* 205 K = N 206 10 CONTINUE 207* 208* If K < 1, exit from loop. 209* 210 IF( K.LT.1 ) 211 $ GO TO 30 212* 213 IF( IPIV( K ).GT.0 ) THEN 214* 215* 1 x 1 diagonal block 216* 217* Interchange rows K and IPIV(K). 218* 219 KP = IPIV( K ) 220 IF( KP.NE.K ) 221 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 222* 223* Multiply by inv(U(K)), where U(K) is the transformation 224* stored in column K of A. 225* 226 CALL ZGERU( K-1, NRHS, -CONE, A( 1, K ), 1, B( K, 1 ), LDB, 227 $ B( 1, 1 ), LDB ) 228* 229* Multiply by the inverse of the diagonal block. 230* 231 CALL ZSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB ) 232 K = K - 1 233 ELSE 234* 235* 2 x 2 diagonal block 236* 237* Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1) 238* 239 KP = -IPIV( K ) 240 IF( KP.NE.K ) 241 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 242* 243 KP = -IPIV( K-1 ) 244 IF( KP.NE.K-1 ) 245 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 246* 247* Multiply by inv(U(K)), where U(K) is the transformation 248* stored in columns K-1 and K of A. 249* 250 IF( K.GT.2 ) THEN 251 CALL ZGERU( K-2, NRHS,-CONE, A( 1, K ), 1, B( K, 1 ), 252 $ LDB, B( 1, 1 ), LDB ) 253 CALL ZGERU( K-2, NRHS,-CONE, A( 1, K-1 ), 1, B( K-1, 1 ), 254 $ LDB, B( 1, 1 ), LDB ) 255 END IF 256* 257* Multiply by the inverse of the diagonal block. 258* 259 AKM1K = A( K-1, K ) 260 AKM1 = A( K-1, K-1 ) / AKM1K 261 AK = A( K, K ) / AKM1K 262 DENOM = AKM1*AK - CONE 263 DO 20 J = 1, NRHS 264 BKM1 = B( K-1, J ) / AKM1K 265 BK = B( K, J ) / AKM1K 266 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM 267 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM 268 20 CONTINUE 269 K = K - 2 270 END IF 271* 272 GO TO 10 273 30 CONTINUE 274* 275* Next solve U**T *X = B, overwriting B with X. 276* 277* K is the main loop index, increasing from 1 to N in steps of 278* 1 or 2, depending on the size of the diagonal blocks. 279* 280 K = 1 281 40 CONTINUE 282* 283* If K > N, exit from loop. 284* 285 IF( K.GT.N ) 286 $ GO TO 50 287* 288 IF( IPIV( K ).GT.0 ) THEN 289* 290* 1 x 1 diagonal block 291* 292* Multiply by inv(U**T(K)), where U(K) is the transformation 293* stored in column K of A. 294* 295 IF( K.GT.1 ) 296 $ CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B, 297 $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB ) 298* 299* Interchange rows K and IPIV(K). 300* 301 KP = IPIV( K ) 302 IF( KP.NE.K ) 303 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 304 K = K + 1 305 ELSE 306* 307* 2 x 2 diagonal block 308* 309* Multiply by inv(U**T(K+1)), where U(K+1) is the transformation 310* stored in columns K and K+1 of A. 311* 312 IF( K.GT.1 ) THEN 313 CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B, 314 $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB ) 315 CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B, 316 $ LDB, A( 1, K+1 ), 1, CONE, B( K+1, 1 ), LDB ) 317 END IF 318* 319* Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1). 320* 321 KP = -IPIV( K ) 322 IF( KP.NE.K ) 323 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 324* 325 KP = -IPIV( K+1 ) 326 IF( KP.NE.K+1 ) 327 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 328* 329 K = K + 2 330 END IF 331* 332 GO TO 40 333 50 CONTINUE 334* 335 ELSE 336* 337* Solve A*X = B, where A = L*D*L**T. 338* 339* First solve L*D*X = B, overwriting B with X. 340* 341* K is the main loop index, increasing from 1 to N in steps of 342* 1 or 2, depending on the size of the diagonal blocks. 343* 344 K = 1 345 60 CONTINUE 346* 347* If K > N, exit from loop. 348* 349 IF( K.GT.N ) 350 $ GO TO 80 351* 352 IF( IPIV( K ).GT.0 ) THEN 353* 354* 1 x 1 diagonal block 355* 356* Interchange rows K and IPIV(K). 357* 358 KP = IPIV( K ) 359 IF( KP.NE.K ) 360 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 361* 362* Multiply by inv(L(K)), where L(K) is the transformation 363* stored in column K of A. 364* 365 IF( K.LT.N ) 366 $ CALL ZGERU( N-K, NRHS, -CONE, A( K+1, K ), 1, B( K, 1 ), 367 $ LDB, B( K+1, 1 ), LDB ) 368* 369* Multiply by the inverse of the diagonal block. 370* 371 CALL ZSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB ) 372 K = K + 1 373 ELSE 374* 375* 2 x 2 diagonal block 376* 377* Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1) 378* 379 KP = -IPIV( K ) 380 IF( KP.NE.K ) 381 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 382* 383 KP = -IPIV( K+1 ) 384 IF( KP.NE.K+1 ) 385 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 386* 387* Multiply by inv(L(K)), where L(K) is the transformation 388* stored in columns K and K+1 of A. 389* 390 IF( K.LT.N-1 ) THEN 391 CALL ZGERU( N-K-1, NRHS,-CONE, A( K+2, K ), 1, B( K, 1 ), 392 $ LDB, B( K+2, 1 ), LDB ) 393 CALL ZGERU( N-K-1, NRHS,-CONE, A( K+2, K+1 ), 1, 394 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) 395 END IF 396* 397* Multiply by the inverse of the diagonal block. 398* 399 AKM1K = A( K+1, K ) 400 AKM1 = A( K, K ) / AKM1K 401 AK = A( K+1, K+1 ) / AKM1K 402 DENOM = AKM1*AK - CONE 403 DO 70 J = 1, NRHS 404 BKM1 = B( K, J ) / AKM1K 405 BK = B( K+1, J ) / AKM1K 406 B( K, J ) = ( AK*BKM1-BK ) / DENOM 407 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 408 70 CONTINUE 409 K = K + 2 410 END IF 411* 412 GO TO 60 413 80 CONTINUE 414* 415* Next solve L**T *X = B, overwriting B with X. 416* 417* K is the main loop index, decreasing from N to 1 in steps of 418* 1 or 2, depending on the size of the diagonal blocks. 419* 420 K = N 421 90 CONTINUE 422* 423* If K < 1, exit from loop. 424* 425 IF( K.LT.1 ) 426 $ GO TO 100 427* 428 IF( IPIV( K ).GT.0 ) THEN 429* 430* 1 x 1 diagonal block 431* 432* Multiply by inv(L**T(K)), where L(K) is the transformation 433* stored in column K of A. 434* 435 IF( K.LT.N ) 436 $ CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), 437 $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB ) 438* 439* Interchange rows K and IPIV(K). 440* 441 KP = IPIV( K ) 442 IF( KP.NE.K ) 443 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 444 K = K - 1 445 ELSE 446* 447* 2 x 2 diagonal block 448* 449* Multiply by inv(L**T(K-1)), where L(K-1) is the transformation 450* stored in columns K-1 and K of A. 451* 452 IF( K.LT.N ) THEN 453 CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), 454 $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB ) 455 CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), 456 $ LDB, A( K+1, K-1 ), 1, CONE, B( K-1, 1 ), 457 $ LDB ) 458 END IF 459* 460* Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1) 461* 462 KP = -IPIV( K ) 463 IF( KP.NE.K ) 464 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 465* 466 KP = -IPIV( K-1 ) 467 IF( KP.NE.K-1 ) 468 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 469* 470 K = K - 2 471 END IF 472* 473 GO TO 90 474 100 CONTINUE 475 END IF 476* 477 RETURN 478* 479* End of ZSYTRS_ROOK 480* 481 END 482