src/LAPACK/dsytrs.f

   1       SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          UPLO
   9       INTEGER            INFO, LDA, LDB, N, NRHS
  10 *     ..
  11 *     .. Array Arguments ..
  12       INTEGER            IPIV( * )
  13       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  DSYTRS solves a system of linear equations A*X = B with a real
  20 *  symmetric matrix A using the factorization A = U*D*U**T or
  21 *  A = L*D*L**T computed by DSYTRF.
  22 *
  23 *  Arguments
  24 *  =========
  25 *
  26 *  UPLO    (input) CHARACTER*1
  27 *          Specifies whether the details of the factorization are stored
  28 *          as an upper or lower triangular matrix.
  29 *          = 'U':  Upper triangular, form is A = U*D*U**T;
  30 *          = 'L':  Lower triangular, form is A = L*D*L**T.
  31 *
  32 *  N       (input) INTEGER
  33 *          The order of the matrix A.  N >= 0.
  34 *
  35 *  NRHS    (input) INTEGER
  36 *          The number of right hand sides, i.e., the number of columns
  37 *          of the matrix B.  NRHS >= 0.
  38 *
  39 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
  40 *          The block diagonal matrix D and the multipliers used to
  41 *          obtain the factor U or L as computed by DSYTRF.
  42 *
  43 *  LDA     (input) INTEGER
  44 *          The leading dimension of the array A.  LDA >= max(1,N).
  45 *
  46 *  IPIV    (input) INTEGER array, dimension (N)
  47 *          Details of the interchanges and the block structure of D
  48 *          as determined by DSYTRF.
  49 *
  50 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  51 *          On entry, the right hand side matrix B.
  52 *          On exit, the solution matrix X.
  53 *
  54 *  LDB     (input) INTEGER
  55 *          The leading dimension of the array B.  LDB >= max(1,N).
  56 *
  57 *  INFO    (output) INTEGER
  58 *          = 0:  successful exit
  59 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  60 *
  61 *  =====================================================================
  62 *
  63 *     .. Parameters ..
  64       DOUBLE PRECISION   ONE
  65       PARAMETER          ( ONE = 1.0D+0 )
  66 *     ..
  67 *     .. Local Scalars ..
  68       LOGICAL            UPPER
  69       INTEGER            J, K, KP
  70       DOUBLE PRECISION   AK, AKM1, AKM1K, BK, BKM1, DENOM
  71 *     ..
  72 *     .. External Functions ..
  73       LOGICAL            LSAME
  74       EXTERNAL           LSAME
  75 *     ..
  76 *     .. External Subroutines ..
  77       EXTERNAL           DGEMV, DGER, DSCAL, DSWAP, XERBLA
  78 *     ..
  79 *     .. Intrinsic Functions ..
  80       INTRINSIC          MAX
  81 *     ..
  82 *     .. Executable Statements ..
  83 *
  84       INFO = 0
  85       UPPER = LSAME( UPLO, 'U' )
  86       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  87          INFO = -1
  88       ELSE IF( N.LT.0 ) THEN
  89          INFO = -2
  90       ELSE IF( NRHS.LT.0 ) THEN
  91          INFO = -3
  92       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  93          INFO = -5
  94       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  95          INFO = -8
  96       END IF
  97       IF( INFO.NE.0 ) THEN
  98          CALL XERBLA( 'DSYTRS', -INFO )
  99          RETURN
 100       END IF
 101 *
 102 *     Quick return if possible
 103 *
 104       IF( N.EQ.0 .OR. NRHS.EQ.0 )
 105      $   RETURN
 106 *
 107       IF( UPPER ) THEN
 108 *
 109 *        Solve A*X = B, where A = U*D*U'.
 110 *
 111 *        First solve U*D*X = B, overwriting B with X.
 112 *
 113 *        K is the main loop index, decreasing from N to 1 in steps of
 114 *        1 or 2, depending on the size of the diagonal blocks.
 115 *
 116          K = N
 117    10    CONTINUE
 118 *
 119 *        If K < 1, exit from loop.
 120 *
 121          IF( K.LT.1 )
 122      $      GO TO 30
 123 *
 124          IF( IPIV( K ).GT.0 ) THEN
 125 *
 126 *           1 x 1 diagonal block
 127 *
 128 *           Interchange rows K and IPIV(K).
 129 *
 130             KP = IPIV( K )
 131             IF( KP.NE.K )
 132      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 133 *
 134 *           Multiply by inv(U(K)), where U(K) is the transformation
 135 *           stored in column K of A.
 136 *
 137             CALL DGER( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
 138      $                 B( 1, 1 ), LDB )
 139 *
 140 *           Multiply by the inverse of the diagonal block.
 141 *
 142             CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
 143             K = K - 1
 144          ELSE
 145 *
 146 *           2 x 2 diagonal block
 147 *
 148 *           Interchange rows K-1 and -IPIV(K).
 149 *
 150             KP = -IPIV( K )
 151             IF( KP.NE.K-1 )
 152      $         CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
 153 *
 154 *           Multiply by inv(U(K)), where U(K) is the transformation
 155 *           stored in columns K-1 and K of A.
 156 *
 157             CALL DGER( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
 158      $                 B( 1, 1 ), LDB )
 159             CALL DGER( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
 160      $                 LDB, B( 1, 1 ), LDB )
 161 *
 162 *           Multiply by the inverse of the diagonal block.
 163 *
 164             AKM1K = A( K-1, K )
 165             AKM1 = A( K-1, K-1 ) / AKM1K
 166             AK = A( K, K ) / AKM1K
 167             DENOM = AKM1*AK - ONE
 168             DO 20 J = 1, NRHS
 169                BKM1 = B( K-1, J ) / AKM1K
 170                BK = B( K, J ) / AKM1K
 171                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
 172                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
 173    20       CONTINUE
 174             K = K - 2
 175          END IF
 176 *
 177          GO TO 10
 178    30    CONTINUE
 179 *
 180 *        Next solve U'*X = B, overwriting B with X.
 181 *
 182 *        K is the main loop index, increasing from 1 to N in steps of
 183 *        1 or 2, depending on the size of the diagonal blocks.
 184 *
 185          K = 1
 186    40    CONTINUE
 187 *
 188 *        If K > N, exit from loop.
 189 *
 190          IF( K.GT.N )
 191      $      GO TO 50
 192 *
 193          IF( IPIV( K ).GT.0 ) THEN
 194 *
 195 *           1 x 1 diagonal block
 196 *
 197 *           Multiply by inv(U'(K)), where U(K) is the transformation
 198 *           stored in column K of A.
 199 *
 200             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ),
 201      $                  1, ONE, B( K, 1 ), LDB )
 202 *
 203 *           Interchange rows K and IPIV(K).
 204 *
 205             KP = IPIV( K )
 206             IF( KP.NE.K )
 207      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 208             K = K + 1
 209          ELSE
 210 *
 211 *           2 x 2 diagonal block
 212 *
 213 *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
 214 *           stored in columns K and K+1 of A.
 215 *
 216             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ),
 217      $                  1, ONE, B( K, 1 ), LDB )
 218             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
 219      $                  A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
 220 *
 221 *           Interchange rows K and -IPIV(K).
 222 *
 223             KP = -IPIV( K )
 224             IF( KP.NE.K )
 225      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 226             K = K + 2
 227          END IF
 228 *
 229          GO TO 40
 230    50    CONTINUE
 231 *
 232       ELSE
 233 *
 234 *        Solve A*X = B, where A = L*D*L'.
 235 *
 236 *        First solve L*D*X = B, overwriting B with X.
 237 *
 238 *        K is the main loop index, increasing from 1 to N in steps of
 239 *        1 or 2, depending on the size of the diagonal blocks.
 240 *
 241          K = 1
 242    60    CONTINUE
 243 *
 244 *        If K > N, exit from loop.
 245 *
 246          IF( K.GT.N )
 247      $      GO TO 80
 248 *
 249          IF( IPIV( K ).GT.0 ) THEN
 250 *
 251 *           1 x 1 diagonal block
 252 *
 253 *           Interchange rows K and IPIV(K).
 254 *
 255             KP = IPIV( K )
 256             IF( KP.NE.K )
 257      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 258 *
 259 *           Multiply by inv(L(K)), where L(K) is the transformation
 260 *           stored in column K of A.
 261 *
 262             IF( K.LT.N )
 263      $         CALL DGER( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
 264      $                    LDB, B( K+1, 1 ), LDB )
 265 *
 266 *           Multiply by the inverse of the diagonal block.
 267 *
 268             CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
 269             K = K + 1
 270          ELSE
 271 *
 272 *           2 x 2 diagonal block
 273 *
 274 *           Interchange rows K+1 and -IPIV(K).
 275 *
 276             KP = -IPIV( K )
 277             IF( KP.NE.K+1 )
 278      $         CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
 279 *
 280 *           Multiply by inv(L(K)), where L(K) is the transformation
 281 *           stored in columns K and K+1 of A.
 282 *
 283             IF( K.LT.N-1 ) THEN
 284                CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
 285      $                    LDB, B( K+2, 1 ), LDB )
 286                CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
 287      $                    B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
 288             END IF
 289 *
 290 *           Multiply by the inverse of the diagonal block.
 291 *
 292             AKM1K = A( K+1, K )
 293             AKM1 = A( K, K ) / AKM1K
 294             AK = A( K+1, K+1 ) / AKM1K
 295             DENOM = AKM1*AK - ONE
 296             DO 70 J = 1, NRHS
 297                BKM1 = B( K, J ) / AKM1K
 298                BK = B( K+1, J ) / AKM1K
 299                B( K, J ) = ( AK*BKM1-BK ) / DENOM
 300                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
 301    70       CONTINUE
 302             K = K + 2
 303          END IF
 304 *
 305          GO TO 60
 306    80    CONTINUE
 307 *
 308 *        Next solve L'*X = B, overwriting B with X.
 309 *
 310 *        K is the main loop index, decreasing from N to 1 in steps of
 311 *        1 or 2, depending on the size of the diagonal blocks.
 312 *
 313          K = N
 314    90    CONTINUE
 315 *
 316 *        If K < 1, exit from loop.
 317 *
 318          IF( K.LT.1 )
 319      $      GO TO 100
 320 *
 321          IF( IPIV( K ).GT.0 ) THEN
 322 *
 323 *           1 x 1 diagonal block
 324 *
 325 *           Multiply by inv(L'(K)), where L(K) is the transformation
 326 *           stored in column K of A.
 327 *
 328             IF( K.LT.N )
 329      $         CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
 330      $                     LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
 331 *
 332 *           Interchange rows K and IPIV(K).
 333 *
 334             KP = IPIV( K )
 335             IF( KP.NE.K )
 336      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 337             K = K - 1
 338          ELSE
 339 *
 340 *           2 x 2 diagonal block
 341 *
 342 *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
 343 *           stored in columns K-1 and K of A.
 344 *
 345             IF( K.LT.N ) THEN
 346                CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
 347      $                     LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
 348                CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
 349      $                     LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ),
 350      $                     LDB )
 351             END IF
 352 *
 353 *           Interchange rows K and -IPIV(K).
 354 *
 355             KP = -IPIV( K )
 356             IF( KP.NE.K )
 357      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
 358             K = K - 2
 359          END IF
 360 *
 361          GO TO 90
 362   100    CONTINUE
 363       END IF
 364 *
 365       RETURN
 366 *
 367 *     End of DSYTRS
 368 *
 369       END