SRC/zhetrs_rook.f

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