SRC/dsptrs.f

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