SRC/ssytrs2.f

   1 *> \brief \b SSYTRS2
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSYTRS2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrs2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrs2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrs2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
  22 *                           WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDA, LDB, N, NRHS
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IPIV( * )
  30 *       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SSYTRS2 solves a system of linear equations A*X = B with a real
  40 *> symmetric matrix A using the factorization A = U*D*U**T or
  41 *> A = L*D*L**T computed by SSYTRF and converted by SSYCONV.
  42 *> \endverbatim
  43 *
  44 *  Arguments:
  45 *  ==========
  46 *
  47 *> \param[in] UPLO
  48 *> \verbatim
  49 *>          UPLO is CHARACTER*1
  50 *>          Specifies whether the details of the factorization are stored
  51 *>          as an upper or lower triangular matrix.
  52 *>          = 'U':  Upper triangular, form is A = U*D*U**T;
  53 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] N
  57 *> \verbatim
  58 *>          N is INTEGER
  59 *>          The order of the matrix A.  N >= 0.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] NRHS
  63 *> \verbatim
  64 *>          NRHS is INTEGER
  65 *>          The number of right hand sides, i.e., the number of columns
  66 *>          of the matrix B.  NRHS >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is REAL array, dimension (LDA,N)
  72 *>          The block diagonal matrix D and the multipliers used to
  73 *>          obtain the factor U or L as computed by SSYTRF.
  74 *>          Note that A is input / output. This might be counter-intuitive,
  75 *>          and one may think that A is input only. A is input / output. This
  76 *>          is because, at the start of the subroutine, we permute A in a
  77 *>          "better" form and then we permute A back to its original form at
  78 *>          the end.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] LDA
  82 *> \verbatim
  83 *>          LDA is INTEGER
  84 *>          The leading dimension of the array A.  LDA >= max(1,N).
  85 *> \endverbatim
  86 *>
  87 *> \param[in] IPIV
  88 *> \verbatim
  89 *>          IPIV is INTEGER array, dimension (N)
  90 *>          Details of the interchanges and the block structure of D
  91 *>          as determined by SSYTRF.
  92 *> \endverbatim
  93 *>
  94 *> \param[in,out] B
  95 *> \verbatim
  96 *>          B is REAL array, dimension (LDB,NRHS)
  97 *>          On entry, the right hand side matrix B.
  98 *>          On exit, the solution matrix X.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDB
 102 *> \verbatim
 103 *>          LDB is INTEGER
 104 *>          The leading dimension of the array B.  LDB >= max(1,N).
 105 *> \endverbatim
 106 *>
 107 *> \param[out] WORK
 108 *> \verbatim
 109 *>          WORK is REAL array, dimension (N)
 110 *> \endverbatim
 111 *>
 112 *> \param[out] INFO
 113 *> \verbatim
 114 *>          INFO is INTEGER
 115 *>          = 0:  successful exit
 116 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 117 *> \endverbatim
 118 *
 119 *  Authors:
 120 *  ========
 121 *
 122 *> \author Univ. of Tennessee
 123 *> \author Univ. of California Berkeley
 124 *> \author Univ. of Colorado Denver
 125 *> \author NAG Ltd.
 126 *
 127 *> \date November 2015
 128 *
 129 *> \ingroup realSYcomputational
 130 *
 131 *  =====================================================================
 132       SUBROUTINE SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
 133      $                    WORK, INFO )
 134 *
 135 *  -- LAPACK computational routine (version 3.6.0) --
 136 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 137 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 138 *     November 2015
 139 *
 140 *     .. Scalar Arguments ..
 141       CHARACTER          UPLO
 142       INTEGER            INFO, LDA, LDB, N, NRHS
 143 *     ..
 144 *     .. Array Arguments ..
 145       INTEGER            IPIV( * )
 146       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
 147 *     ..
 148 *
 149 *  =====================================================================
 150 *
 151 *     .. Parameters ..
 152       REAL               ONE
 153       PARAMETER          ( ONE = 1.0E+0 )
 154 *     ..
 155 *     .. Local Scalars ..
 156       LOGICAL            UPPER
 157       INTEGER            I, IINFO, J, K, KP
 158       REAL               AK, AKM1, AKM1K, BK, BKM1, DENOM
 159 *     ..
 160 *     .. External Functions ..
 161       LOGICAL            LSAME
 162       EXTERNAL           LSAME
 163 *     ..
 164 *     .. External Subroutines ..
 165       EXTERNAL           SSCAL, SSYCONV, SSWAP, STRSM, XERBLA
 166 *     ..
 167 *     .. Intrinsic Functions ..
 168       INTRINSIC          MAX
 169 *     ..
 170 *     .. Executable Statements ..
 171 *
 172       INFO = 0
 173       UPPER = LSAME( UPLO, 'U' )
 174       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 175          INFO = -1
 176       ELSE IF( N.LT.0 ) THEN
 177          INFO = -2
 178       ELSE IF( NRHS.LT.0 ) THEN
 179          INFO = -3
 180       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 181          INFO = -5
 182       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 183          INFO = -8
 184       END IF
 185       IF( INFO.NE.0 ) THEN
 186          CALL XERBLA( 'SSYTRS2', -INFO )
 187          RETURN
 188       END IF
 189 *
 190 *     Quick return if possible
 191 *
 192       IF( N.EQ.0 .OR. NRHS.EQ.0 )
 193      $   RETURN
 194 *
 195 *     Convert A
 196 *
 197       CALL SSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
 198 *
 199       IF( UPPER ) THEN
 200 *
 201 *        Solve A*X = B, where A = U*D*U**T.
 202 *
 203 *       P**T * B
 204         K=N
 205         DO WHILE ( K .GE. 1 )
 206          IF( IPIV( K ).GT.0 ) THEN
 207 *           1 x 1 diagonal block
 208 *           Interchange rows K and IPIV(K).
 209             KP = IPIV( K )
 210             IF( KP.NE.K )
 211      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 212             K=K-1
 213          ELSE
 214 *           2 x 2 diagonal block
 215 *           Interchange rows K-1 and -IPIV(K).
 216             KP = -IPIV( K )
 217             IF( KP.EQ.-IPIV( K-1 ) )
 218      $         CALL SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
 219             K=K-2
 220          END IF
 221         END DO
 222 *
 223 *  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
 224 *
 225         CALL STRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
 226 *
 227 *  Compute D \ B -> B   [ D \ (U \P**T * B) ]
 228 *
 229          I=N
 230          DO WHILE ( I .GE. 1 )
 231             IF( IPIV(I) .GT. 0 ) THEN
 232               CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
 233             ELSEIF ( I .GT. 1) THEN
 234                IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
 235                   AKM1K = WORK(I)
 236                   AKM1 = A( I-1, I-1 ) / AKM1K
 237                   AK = A( I, I ) / AKM1K
 238                   DENOM = AKM1*AK - ONE
 239                   DO 15 J = 1, NRHS
 240                      BKM1 = B( I-1, J ) / AKM1K
 241                      BK = B( I, J ) / AKM1K
 242                      B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
 243                      B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
 244  15              CONTINUE
 245                I = I - 1
 246                ENDIF
 247             ENDIF
 248             I = I - 1
 249          END DO
 250 *
 251 *      Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
 252 *
 253          CALL STRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
 254 *
 255 *       P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
 256 *
 257         K=1
 258         DO WHILE ( K .LE. N )
 259          IF( IPIV( K ).GT.0 ) THEN
 260 *           1 x 1 diagonal block
 261 *           Interchange rows K and IPIV(K).
 262             KP = IPIV( K )
 263             IF( KP.NE.K )
 264      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 265             K=K+1
 266          ELSE
 267 *           2 x 2 diagonal block
 268 *           Interchange rows K-1 and -IPIV(K).
 269             KP = -IPIV( K )
 270             IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
 271      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 272             K=K+2
 273          ENDIF
 274         END DO
 275 *
 276       ELSE
 277 *
 278 *        Solve A*X = B, where A = L*D*L**T.
 279 *
 280 *       P**T * B
 281         K=1
 282         DO WHILE ( K .LE. N )
 283          IF( IPIV( K ).GT.0 ) THEN
 284 *           1 x 1 diagonal block
 285 *           Interchange rows K and IPIV(K).
 286             KP = IPIV( K )
 287             IF( KP.NE.K )
 288      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 289             K=K+1
 290          ELSE
 291 *           2 x 2 diagonal block
 292 *           Interchange rows K and -IPIV(K+1).
 293             KP = -IPIV( K+1 )
 294             IF( KP.EQ.-IPIV( K ) )
 295      $         CALL SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
 296             K=K+2
 297          ENDIF
 298         END DO
 299 *
 300 *  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
 301 *
 302         CALL STRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
 303 *
 304 *  Compute D \ B -> B   [ D \ (L \P**T * B) ]
 305 *
 306          I=1
 307          DO WHILE ( I .LE. N )
 308             IF( IPIV(I) .GT. 0 ) THEN
 309               CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
 310             ELSE
 311                   AKM1K = WORK(I)
 312                   AKM1 = A( I, I ) / AKM1K
 313                   AK = A( I+1, I+1 ) / AKM1K
 314                   DENOM = AKM1*AK - ONE
 315                   DO 25 J = 1, NRHS
 316                      BKM1 = B( I, J ) / AKM1K
 317                      BK = B( I+1, J ) / AKM1K
 318                      B( I, J ) = ( AK*BKM1-BK ) / DENOM
 319                      B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
 320  25              CONTINUE
 321                   I = I + 1
 322             ENDIF
 323             I = I + 1
 324          END DO
 325 *
 326 *  Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
 327 *
 328         CALL STRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
 329 *
 330 *       P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
 331 *
 332         K=N
 333         DO WHILE ( K .GE. 1 )
 334          IF( IPIV( K ).GT.0 ) THEN
 335 *           1 x 1 diagonal block
 336 *           Interchange rows K and IPIV(K).
 337             KP = IPIV( K )
 338             IF( KP.NE.K )
 339      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 340             K=K-1
 341          ELSE
 342 *           2 x 2 diagonal block
 343 *           Interchange rows K-1 and -IPIV(K).
 344             KP = -IPIV( K )
 345             IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
 346      $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 347             K=K-2
 348          ENDIF
 349         END DO
 350 *
 351       END IF
 352 *
 353 *     Revert A
 354 *
 355       CALL SSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 356 *
 357       RETURN
 358 *
 359 *     End of SSYTRS2
 360 *
 361       END