Return to zgtrfs.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

    1:       SUBROUTINE ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2,
    2:      $                   IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK,
    3:      $                   INFO )
    4: *
    5: *  -- LAPACK routine (version 3.2) --
    6: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    7: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    8: *     November 2006
    9: *
   10: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
   11: *
   12: *     .. Scalar Arguments ..
   13:       CHARACTER          TRANS
   14:       INTEGER            INFO, LDB, LDX, N, NRHS
   15: *     ..
   16: *     .. Array Arguments ..
   17:       INTEGER            IPIV( * )
   18:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
   19:       COMPLEX*16         B( LDB, * ), D( * ), DF( * ), DL( * ),
   20:      $                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
   21:      $                   WORK( * ), X( LDX, * )
   22: *     ..
   23: *
   24: *  Purpose
   25: *  =======
   26: *
   27: *  ZGTRFS improves the computed solution to a system of linear
   28: *  equations when the coefficient matrix is tridiagonal, and provides
   29: *  error bounds and backward error estimates for the solution.
   30: *
   31: *  Arguments
   32: *  =========
   33: *
   34: *  TRANS   (input) CHARACTER*1
   35: *          Specifies the form of the system of equations:
   36: *          = 'N':  A * X = B     (No transpose)
   37: *          = 'T':  A**T * X = B  (Transpose)
   38: *          = 'C':  A**H * X = B  (Conjugate transpose)
   39: *
   40: *  N       (input) INTEGER
   41: *          The order of the matrix A.  N >= 0.
   42: *
   43: *  NRHS    (input) INTEGER
   44: *          The number of right hand sides, i.e., the number of columns
   45: *          of the matrix B.  NRHS >= 0.
   46: *
   47: *  DL      (input) COMPLEX*16 array, dimension (N-1)
   48: *          The (n-1) subdiagonal elements of A.
   49: *
   50: *  D       (input) COMPLEX*16 array, dimension (N)
   51: *          The diagonal elements of A.
   52: *
   53: *  DU      (input) COMPLEX*16 array, dimension (N-1)
   54: *          The (n-1) superdiagonal elements of A.
   55: *
   56: *  DLF     (input) COMPLEX*16 array, dimension (N-1)
   57: *          The (n-1) multipliers that define the matrix L from the
   58: *          LU factorization of A as computed by ZGTTRF.
   59: *
   60: *  DF      (input) COMPLEX*16 array, dimension (N)
   61: *          The n diagonal elements of the upper triangular matrix U from
   62: *          the LU factorization of A.
   63: *
   64: *  DUF     (input) COMPLEX*16 array, dimension (N-1)
   65: *          The (n-1) elements of the first superdiagonal of U.
   66: *
   67: *  DU2     (input) COMPLEX*16 array, dimension (N-2)
   68: *          The (n-2) elements of the second superdiagonal of U.
   69: *
   70: *  IPIV    (input) INTEGER array, dimension (N)
   71: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
   72: *          interchanged with row IPIV(i).  IPIV(i) will always be either
   73: *          i or i+1; IPIV(i) = i indicates a row interchange was not
   74: *          required.
   75: *
   76: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
   77: *          The right hand side matrix B.
   78: *
   79: *  LDB     (input) INTEGER
   80: *          The leading dimension of the array B.  LDB >= max(1,N).
   81: *
   82: *  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
   83: *          On entry, the solution matrix X, as computed by ZGTTRS.
   84: *          On exit, the improved solution matrix X.
   85: *
   86: *  LDX     (input) INTEGER
   87: *          The leading dimension of the array X.  LDX >= max(1,N).
   88: *
   89: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
   90: *          The estimated forward error bound for each solution vector
   91: *          X(j) (the j-th column of the solution matrix X).
   92: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
   93: *          is an estimated upper bound for the magnitude of the largest
   94: *          element in (X(j) - XTRUE) divided by the magnitude of the
   95: *          largest element in X(j).  The estimate is as reliable as
   96: *          the estimate for RCOND, and is almost always a slight
   97: *          overestimate of the true error.
   98: *
   99: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  100: *          The componentwise relative backward error of each solution
  101: *          vector X(j) (i.e., the smallest relative change in
  102: *          any element of A or B that makes X(j) an exact solution).
  103: *
  104: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
  105: *
  106: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
  107: *
  108: *  INFO    (output) INTEGER
  109: *          = 0:  successful exit
  110: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  111: *
  112: *  Internal Parameters
  113: *  ===================
  114: *
  115: *  ITMAX is the maximum number of steps of iterative refinement.
  116: *
  117: *  =====================================================================
  118: *
  119: *     .. Parameters ..
  120:       INTEGER            ITMAX
  121:       PARAMETER          ( ITMAX = 5 )
  122:       DOUBLE PRECISION   ZERO, ONE
  123:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  124:       DOUBLE PRECISION   TWO
  125:       PARAMETER          ( TWO = 2.0D+0 )
  126:       DOUBLE PRECISION   THREE
  127:       PARAMETER          ( THREE = 3.0D+0 )
  128: *     ..
  129: *     .. Local Scalars ..
  130:       LOGICAL            NOTRAN
  131:       CHARACTER          TRANSN, TRANST
  132:       INTEGER            COUNT, I, J, KASE, NZ
  133:       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
  134:       COMPLEX*16         ZDUM
  135: *     ..
  136: *     .. Local Arrays ..
  137:       INTEGER            ISAVE( 3 )
  138: *     ..
  139: *     .. External Subroutines ..
  140:       EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZGTTRS, ZLACN2, ZLAGTM
  141: *     ..
  142: *     .. Intrinsic Functions ..
  143:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, MAX
  144: *     ..
  145: *     .. External Functions ..
  146:       LOGICAL            LSAME
  147:       DOUBLE PRECISION   DLAMCH
  148:       EXTERNAL           LSAME, DLAMCH
  149: *     ..
  150: *     .. Statement Functions ..
  151:       DOUBLE PRECISION   CABS1
  152: *     ..
  153: *     .. Statement Function definitions ..
  154:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  155: *     ..
  156: *     .. Executable Statements ..
  157: *
  158: *     Test the input parameters.
  159: *
  160:       INFO = 0
  161:       NOTRAN = LSAME( TRANS, 'N' )
  162:       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
  163:      $    LSAME( TRANS, 'C' ) ) THEN
  164:          INFO = -1
  165:       ELSE IF( N.LT.0 ) THEN
  166:          INFO = -2
  167:       ELSE IF( NRHS.LT.0 ) THEN
  168:          INFO = -3
  169:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  170:          INFO = -13
  171:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  172:          INFO = -15
  173:       END IF
  174:       IF( INFO.NE.0 ) THEN
  175:          CALL XERBLA( 'ZGTRFS', -INFO )
  176:          RETURN
  177:       END IF
  178: *
  179: *     Quick return if possible
  180: *
  181:       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
  182:          DO 10 J = 1, NRHS
  183:             FERR( J ) = ZERO
  184:             BERR( J ) = ZERO
  185:    10    CONTINUE
  186:          RETURN
  187:       END IF
  188: *
  189:       IF( NOTRAN ) THEN
  190:          TRANSN = 'N'
  191:          TRANST = 'C'
  192:       ELSE
  193:          TRANSN = 'C'
  194:          TRANST = 'N'
  195:       END IF
  196: *
  197: *     NZ = maximum number of nonzero elements in each row of A, plus 1
  198: *
  199:       NZ = 4
  200:       EPS = DLAMCH( 'Epsilon' )
  201:       SAFMIN = DLAMCH( 'Safe minimum' )
  202:       SAFE1 = NZ*SAFMIN
  203:       SAFE2 = SAFE1 / EPS
  204: *
  205: *     Do for each right hand side
  206: *
  207:       DO 110 J = 1, NRHS
  208: *
  209:          COUNT = 1
  210:          LSTRES = THREE
  211:    20    CONTINUE
  212: *
  213: *        Loop until stopping criterion is satisfied.
  214: *
  215: *        Compute residual R = B - op(A) * X,
  216: *        where op(A) = A, A**T, or A**H, depending on TRANS.
  217: *
  218:          CALL ZCOPY( N, B( 1, J ), 1, WORK, 1 )
  219:          CALL ZLAGTM( TRANS, N, 1, -ONE, DL, D, DU, X( 1, J ), LDX, ONE,
  220:      $                WORK, N )
  221: *
  222: *        Compute abs(op(A))*abs(x) + abs(b) for use in the backward
  223: *        error bound.
  224: *
  225:          IF( NOTRAN ) THEN
  226:             IF( N.EQ.1 ) THEN
  227:                RWORK( 1 ) = CABS1( B( 1, J ) ) +
  228:      $                      CABS1( D( 1 ) )*CABS1( X( 1, J ) )
  229:             ELSE
  230:                RWORK( 1 ) = CABS1( B( 1, J ) ) +
  231:      $                      CABS1( D( 1 ) )*CABS1( X( 1, J ) ) +
  232:      $                      CABS1( DU( 1 ) )*CABS1( X( 2, J ) )
  233:                DO 30 I = 2, N - 1
  234:                   RWORK( I ) = CABS1( B( I, J ) ) +
  235:      $                         CABS1( DL( I-1 ) )*CABS1( X( I-1, J ) ) +
  236:      $                         CABS1( D( I ) )*CABS1( X( I, J ) ) +
  237:      $                         CABS1( DU( I ) )*CABS1( X( I+1, J ) )
  238:    30          CONTINUE
  239:                RWORK( N ) = CABS1( B( N, J ) ) +
  240:      $                      CABS1( DL( N-1 ) )*CABS1( X( N-1, J ) ) +
  241:      $                      CABS1( D( N ) )*CABS1( X( N, J ) )
  242:             END IF
  243:          ELSE
  244:             IF( N.EQ.1 ) THEN
  245:                RWORK( 1 ) = CABS1( B( 1, J ) ) +
  246:      $                      CABS1( D( 1 ) )*CABS1( X( 1, J ) )
  247:             ELSE
  248:                RWORK( 1 ) = CABS1( B( 1, J ) ) +
  249:      $                      CABS1( D( 1 ) )*CABS1( X( 1, J ) ) +
  250:      $                      CABS1( DL( 1 ) )*CABS1( X( 2, J ) )
  251:                DO 40 I = 2, N - 1
  252:                   RWORK( I ) = CABS1( B( I, J ) ) +
  253:      $                         CABS1( DU( I-1 ) )*CABS1( X( I-1, J ) ) +
  254:      $                         CABS1( D( I ) )*CABS1( X( I, J ) ) +
  255:      $                         CABS1( DL( I ) )*CABS1( X( I+1, J ) )
  256:    40          CONTINUE
  257:                RWORK( N ) = CABS1( B( N, J ) ) +
  258:      $                      CABS1( DU( N-1 ) )*CABS1( X( N-1, J ) ) +
  259:      $                      CABS1( D( N ) )*CABS1( X( N, J ) )
  260:             END IF
  261:          END IF
  262: *
  263: *        Compute componentwise relative backward error from formula
  264: *
  265: *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
  266: *
  267: *        where abs(Z) is the componentwise absolute value of the matrix
  268: *        or vector Z.  If the i-th component of the denominator is less
  269: *        than SAFE2, then SAFE1 is added to the i-th components of the
  270: *        numerator and denominator before dividing.
  271: *
  272:          S = ZERO
  273:          DO 50 I = 1, N
  274:             IF( RWORK( I ).GT.SAFE2 ) THEN
  275:                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
  276:             ELSE
  277:                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
  278:      $             ( RWORK( I )+SAFE1 ) )
  279:             END IF
  280:    50    CONTINUE
  281:          BERR( J ) = S
  282: *
  283: *        Test stopping criterion. Continue iterating if
  284: *           1) The residual BERR(J) is larger than machine epsilon, and
  285: *           2) BERR(J) decreased by at least a factor of 2 during the
  286: *              last iteration, and
  287: *           3) At most ITMAX iterations tried.
  288: *
  289:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
  290:      $       COUNT.LE.ITMAX ) THEN
  291: *
  292: *           Update solution and try again.
  293: *
  294:             CALL ZGTTRS( TRANS, N, 1, DLF, DF, DUF, DU2, IPIV, WORK, N,
  295:      $                   INFO )
  296:             CALL ZAXPY( N, DCMPLX( ONE ), WORK, 1, X( 1, J ), 1 )
  297:             LSTRES = BERR( J )
  298:             COUNT = COUNT + 1
  299:             GO TO 20
  300:          END IF
  301: *
  302: *        Bound error from formula
  303: *
  304: *        norm(X - XTRUE) / norm(X) .le. FERR =
  305: *        norm( abs(inv(op(A)))*
  306: *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
  307: *
  308: *        where
  309: *          norm(Z) is the magnitude of the largest component of Z
  310: *          inv(op(A)) is the inverse of op(A)
  311: *          abs(Z) is the componentwise absolute value of the matrix or
  312: *             vector Z
  313: *          NZ is the maximum number of nonzeros in any row of A, plus 1
  314: *          EPS is machine epsilon
  315: *
  316: *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
  317: *        is incremented by SAFE1 if the i-th component of
  318: *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
  319: *
  320: *        Use ZLACN2 to estimate the infinity-norm of the matrix
  321: *           inv(op(A)) * diag(W),
  322: *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
  323: *
  324:          DO 60 I = 1, N
  325:             IF( RWORK( I ).GT.SAFE2 ) THEN
  326:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
  327:             ELSE
  328:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
  329:      $                      SAFE1
  330:             END IF
  331:    60    CONTINUE
  332: *
  333:          KASE = 0
  334:    70    CONTINUE
  335:          CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
  336:          IF( KASE.NE.0 ) THEN
  337:             IF( KASE.EQ.1 ) THEN
  338: *
  339: *              Multiply by diag(W)*inv(op(A)**H).
  340: *
  341:                CALL ZGTTRS( TRANST, N, 1, DLF, DF, DUF, DU2, IPIV, WORK,
  342:      $                      N, INFO )
  343:                DO 80 I = 1, N
  344:                   WORK( I ) = RWORK( I )*WORK( I )
  345:    80          CONTINUE
  346:             ELSE
  347: *
  348: *              Multiply by inv(op(A))*diag(W).
  349: *
  350:                DO 90 I = 1, N
  351:                   WORK( I ) = RWORK( I )*WORK( I )
  352:    90          CONTINUE
  353:                CALL ZGTTRS( TRANSN, N, 1, DLF, DF, DUF, DU2, IPIV, WORK,
  354:      $                      N, INFO )
  355:             END IF
  356:             GO TO 70
  357:          END IF
  358: *
  359: *        Normalize error.
  360: *
  361:          LSTRES = ZERO
  362:          DO 100 I = 1, N
  363:             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
  364:   100    CONTINUE
  365:          IF( LSTRES.NE.ZERO )
  366:      $      FERR( J ) = FERR( J ) / LSTRES
  367: *
  368:   110 CONTINUE
  369: *
  370:       RETURN
  371: *
  372: *     End of ZGTRFS
  373: *
  374:       END