SRC/zhptrs.f

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