File:  [local] / rpl / lapack / lapack / zgtts2.f
Revision 1.6: download - view: text, annotated - select for diffs - revision graph
Fri Aug 13 21:04:04 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_0_19, rpl-4_0_18, HEAD
Patches pour OS/2

    1:       SUBROUTINE ZGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            ITRANS, LDB, N, NRHS
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            IPIV( * )
   13:       COMPLEX*16         B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  ZGTTS2 solves one of the systems of equations
   20: *     A * X = B,  A**T * X = B,  or  A**H * X = B,
   21: *  with a tridiagonal matrix A using the LU factorization computed
   22: *  by ZGTTRF.
   23: *
   24: *  Arguments
   25: *  =========
   26: *
   27: *  ITRANS  (input) INTEGER
   28: *          Specifies the form of the system of equations.
   29: *          = 0:  A * X = B     (No transpose)
   30: *          = 1:  A**T * X = B  (Transpose)
   31: *          = 2:  A**H * X = B  (Conjugate transpose)
   32: *
   33: *  N       (input) INTEGER
   34: *          The order of the matrix A.
   35: *
   36: *  NRHS    (input) INTEGER
   37: *          The number of right hand sides, i.e., the number of columns
   38: *          of the matrix B.  NRHS >= 0.
   39: *
   40: *  DL      (input) COMPLEX*16 array, dimension (N-1)
   41: *          The (n-1) multipliers that define the matrix L from the
   42: *          LU factorization of A.
   43: *
   44: *  D       (input) COMPLEX*16 array, dimension (N)
   45: *          The n diagonal elements of the upper triangular matrix U from
   46: *          the LU factorization of A.
   47: *
   48: *  DU      (input) COMPLEX*16 array, dimension (N-1)
   49: *          The (n-1) elements of the first super-diagonal of U.
   50: *
   51: *  DU2     (input) COMPLEX*16 array, dimension (N-2)
   52: *          The (n-2) elements of the second super-diagonal of U.
   53: *
   54: *  IPIV    (input) INTEGER array, dimension (N)
   55: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
   56: *          interchanged with row IPIV(i).  IPIV(i) will always be either
   57: *          i or i+1; IPIV(i) = i indicates a row interchange was not
   58: *          required.
   59: *
   60: *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
   61: *          On entry, the matrix of right hand side vectors B.
   62: *          On exit, B is overwritten by the solution vectors X.
   63: *
   64: *  LDB     (input) INTEGER
   65: *          The leading dimension of the array B.  LDB >= max(1,N).
   66: *
   67: *  =====================================================================
   68: *
   69: *     .. Local Scalars ..
   70:       INTEGER            I, J
   71:       COMPLEX*16         TEMP
   72: *     ..
   73: *     .. Intrinsic Functions ..
   74:       INTRINSIC          DCONJG
   75: *     ..
   76: *     .. Executable Statements ..
   77: *
   78: *     Quick return if possible
   79: *
   80:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
   81:      $   RETURN
   82: *
   83:       IF( ITRANS.EQ.0 ) THEN
   84: *
   85: *        Solve A*X = B using the LU factorization of A,
   86: *        overwriting each right hand side vector with its solution.
   87: *
   88:          IF( NRHS.LE.1 ) THEN
   89:             J = 1
   90:    10       CONTINUE
   91: *
   92: *           Solve L*x = b.
   93: *
   94:             DO 20 I = 1, N - 1
   95:                IF( IPIV( I ).EQ.I ) THEN
   96:                   B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
   97:                ELSE
   98:                   TEMP = B( I, J )
   99:                   B( I, J ) = B( I+1, J )
  100:                   B( I+1, J ) = TEMP - DL( I )*B( I, J )
  101:                END IF
  102:    20       CONTINUE
  103: *
  104: *           Solve U*x = b.
  105: *
  106:             B( N, J ) = B( N, J ) / D( N )
  107:             IF( N.GT.1 )
  108:      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
  109:      $                       D( N-1 )
  110:             DO 30 I = N - 2, 1, -1
  111:                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
  112:      $                     B( I+2, J ) ) / D( I )
  113:    30       CONTINUE
  114:             IF( J.LT.NRHS ) THEN
  115:                J = J + 1
  116:                GO TO 10
  117:             END IF
  118:          ELSE
  119:             DO 60 J = 1, NRHS
  120: *
  121: *           Solve L*x = b.
  122: *
  123:                DO 40 I = 1, N - 1
  124:                   IF( IPIV( I ).EQ.I ) THEN
  125:                      B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
  126:                   ELSE
  127:                      TEMP = B( I, J )
  128:                      B( I, J ) = B( I+1, J )
  129:                      B( I+1, J ) = TEMP - DL( I )*B( I, J )
  130:                   END IF
  131:    40          CONTINUE
  132: *
  133: *           Solve U*x = b.
  134: *
  135:                B( N, J ) = B( N, J ) / D( N )
  136:                IF( N.GT.1 )
  137:      $            B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
  138:      $                          D( N-1 )
  139:                DO 50 I = N - 2, 1, -1
  140:                   B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
  141:      $                        B( I+2, J ) ) / D( I )
  142:    50          CONTINUE
  143:    60       CONTINUE
  144:          END IF
  145:       ELSE IF( ITRANS.EQ.1 ) THEN
  146: *
  147: *        Solve A**T * X = B.
  148: *
  149:          IF( NRHS.LE.1 ) THEN
  150:             J = 1
  151:    70       CONTINUE
  152: *
  153: *           Solve U**T * x = b.
  154: *
  155:             B( 1, J ) = B( 1, J ) / D( 1 )
  156:             IF( N.GT.1 )
  157:      $         B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
  158:             DO 80 I = 3, N
  159:                B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
  160:      $                     B( I-2, J ) ) / D( I )
  161:    80       CONTINUE
  162: *
  163: *           Solve L**T * x = b.
  164: *
  165:             DO 90 I = N - 1, 1, -1
  166:                IF( IPIV( I ).EQ.I ) THEN
  167:                   B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
  168:                ELSE
  169:                   TEMP = B( I+1, J )
  170:                   B( I+1, J ) = B( I, J ) - DL( I )*TEMP
  171:                   B( I, J ) = TEMP
  172:                END IF
  173:    90       CONTINUE
  174:             IF( J.LT.NRHS ) THEN
  175:                J = J + 1
  176:                GO TO 70
  177:             END IF
  178:          ELSE
  179:             DO 120 J = 1, NRHS
  180: *
  181: *           Solve U**T * x = b.
  182: *
  183:                B( 1, J ) = B( 1, J ) / D( 1 )
  184:                IF( N.GT.1 )
  185:      $            B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
  186:                DO 100 I = 3, N
  187:                   B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
  188:      $                        DU2( I-2 )*B( I-2, J ) ) / D( I )
  189:   100          CONTINUE
  190: *
  191: *           Solve L**T * x = b.
  192: *
  193:                DO 110 I = N - 1, 1, -1
  194:                   IF( IPIV( I ).EQ.I ) THEN
  195:                      B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
  196:                   ELSE
  197:                      TEMP = B( I+1, J )
  198:                      B( I+1, J ) = B( I, J ) - DL( I )*TEMP
  199:                      B( I, J ) = TEMP
  200:                   END IF
  201:   110          CONTINUE
  202:   120       CONTINUE
  203:          END IF
  204:       ELSE
  205: *
  206: *        Solve A**H * X = B.
  207: *
  208:          IF( NRHS.LE.1 ) THEN
  209:             J = 1
  210:   130       CONTINUE
  211: *
  212: *           Solve U**H * x = b.
  213: *
  214:             B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) )
  215:             IF( N.GT.1 )
  216:      $         B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) ) /
  217:      $                     DCONJG( D( 2 ) )
  218:             DO 140 I = 3, N
  219:                B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )*B( I-1, J )-
  220:      $                     DCONJG( DU2( I-2 ) )*B( I-2, J ) ) /
  221:      $                     DCONJG( D( I ) )
  222:   140       CONTINUE
  223: *
  224: *           Solve L**H * x = b.
  225: *
  226:             DO 150 I = N - 1, 1, -1
  227:                IF( IPIV( I ).EQ.I ) THEN
  228:                   B( I, J ) = B( I, J ) - DCONJG( DL( I ) )*B( I+1, J )
  229:                ELSE
  230:                   TEMP = B( I+1, J )
  231:                   B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP
  232:                   B( I, J ) = TEMP
  233:                END IF
  234:   150       CONTINUE
  235:             IF( J.LT.NRHS ) THEN
  236:                J = J + 1
  237:                GO TO 130
  238:             END IF
  239:          ELSE
  240:             DO 180 J = 1, NRHS
  241: *
  242: *           Solve U**H * x = b.
  243: *
  244:                B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) )
  245:                IF( N.GT.1 )
  246:      $            B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) )
  247:      $                         / DCONJG( D( 2 ) )
  248:                DO 160 I = 3, N
  249:                   B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )*
  250:      $                        B( I-1, J )-DCONJG( DU2( I-2 ) )*
  251:      $                        B( I-2, J ) ) / DCONJG( D( I ) )
  252:   160          CONTINUE
  253: *
  254: *           Solve L**H * x = b.
  255: *
  256:                DO 170 I = N - 1, 1, -1
  257:                   IF( IPIV( I ).EQ.I ) THEN
  258:                      B( I, J ) = B( I, J ) - DCONJG( DL( I ) )*
  259:      $                           B( I+1, J )
  260:                   ELSE
  261:                      TEMP = B( I+1, J )
  262:                      B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP
  263:                      B( I, J ) = TEMP
  264:                   END IF
  265:   170          CONTINUE
  266:   180       CONTINUE
  267:          END IF
  268:       END IF
  269: *
  270: *     End of ZGTTS2
  271: *
  272:       END

CVSweb interface <joel.bertrand@systella.fr>