SRC/ssytrs_rook.f

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