SRC/zhetrs.f

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