SRC/zsytrs.f

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