Return to zptsvx.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE ZPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
    2:      $                   RCOND, FERR, BERR, WORK, RWORK, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       CHARACTER          FACT
   11:       INTEGER            INFO, LDB, LDX, N, NRHS
   12:       DOUBLE PRECISION   RCOND
   13: *     ..
   14: *     .. Array Arguments ..
   15:       DOUBLE PRECISION   BERR( * ), D( * ), DF( * ), FERR( * ),
   16:      $                   RWORK( * )
   17:       COMPLEX*16         B( LDB, * ), E( * ), EF( * ), WORK( * ),
   18:      $                   X( LDX, * )
   19: *     ..
   20: *
   21: *  Purpose
   22: *  =======
   23: *
   24: *  ZPTSVX uses the factorization A = L*D*L**H to compute the solution
   25: *  to a complex system of linear equations A*X = B, where A is an
   26: *  N-by-N Hermitian positive definite tridiagonal matrix and X and B
   27: *  are N-by-NRHS matrices.
   28: *
   29: *  Error bounds on the solution and a condition estimate are also
   30: *  provided.
   31: *
   32: *  Description
   33: *  ===========
   34: *
   35: *  The following steps are performed:
   36: *
   37: *  1. If FACT = 'N', the matrix A is factored as A = L*D*L**H, where L
   38: *     is a unit lower bidiagonal matrix and D is diagonal.  The
   39: *     factorization can also be regarded as having the form
   40: *     A = U**H*D*U.
   41: *
   42: *  2. If the leading i-by-i principal minor is not positive definite,
   43: *     then the routine returns with INFO = i. Otherwise, the factored
   44: *     form of A is used to estimate the condition number of the matrix
   45: *     A.  If the reciprocal of the condition number is less than machine
   46: *     precision, INFO = N+1 is returned as a warning, but the routine
   47: *     still goes on to solve for X and compute error bounds as
   48: *     described below.
   49: *
   50: *  3. The system of equations is solved for X using the factored form
   51: *     of A.
   52: *
   53: *  4. Iterative refinement is applied to improve the computed solution
   54: *     matrix and calculate error bounds and backward error estimates
   55: *     for it.
   56: *
   57: *  Arguments
   58: *  =========
   59: *
   60: *  FACT    (input) CHARACTER*1
   61: *          Specifies whether or not the factored form of the matrix
   62: *          A is supplied on entry.
   63: *          = 'F':  On entry, DF and EF contain the factored form of A.
   64: *                  D, E, DF, and EF will not be modified.
   65: *          = 'N':  The matrix A will be copied to DF and EF and
   66: *                  factored.
   67: *
   68: *  N       (input) INTEGER
   69: *          The order of the matrix A.  N >= 0.
   70: *
   71: *  NRHS    (input) INTEGER
   72: *          The number of right hand sides, i.e., the number of columns
   73: *          of the matrices B and X.  NRHS >= 0.
   74: *
   75: *  D       (input) DOUBLE PRECISION array, dimension (N)
   76: *          The n diagonal elements of the tridiagonal matrix A.
   77: *
   78: *  E       (input) COMPLEX*16 array, dimension (N-1)
   79: *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
   80: *
   81: *  DF      (input or output) DOUBLE PRECISION array, dimension (N)
   82: *          If FACT = 'F', then DF is an input argument and on entry
   83: *          contains the n diagonal elements of the diagonal matrix D
   84: *          from the L*D*L**H factorization of A.
   85: *          If FACT = 'N', then DF is an output argument and on exit
   86: *          contains the n diagonal elements of the diagonal matrix D
   87: *          from the L*D*L**H factorization of A.
   88: *
   89: *  EF      (input or output) COMPLEX*16 array, dimension (N-1)
   90: *          If FACT = 'F', then EF is an input argument and on entry
   91: *          contains the (n-1) subdiagonal elements of the unit
   92: *          bidiagonal factor L from the L*D*L**H factorization of A.
   93: *          If FACT = 'N', then EF is an output argument and on exit
   94: *          contains the (n-1) subdiagonal elements of the unit
   95: *          bidiagonal factor L from the L*D*L**H factorization of A.
   96: *
   97: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
   98: *          The N-by-NRHS right hand side matrix B.
   99: *
  100: *  LDB     (input) INTEGER
  101: *          The leading dimension of the array B.  LDB >= max(1,N).
  102: *
  103: *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
  104: *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
  105: *
  106: *  LDX     (input) INTEGER
  107: *          The leading dimension of the array X.  LDX >= max(1,N).
  108: *
  109: *  RCOND   (output) DOUBLE PRECISION
  110: *          The reciprocal condition number of the matrix A.  If RCOND
  111: *          is less than the machine precision (in particular, if
  112: *          RCOND = 0), the matrix is singular to working precision.
  113: *          This condition is indicated by a return code of INFO > 0.
  114: *
  115: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  116: *          The forward error bound for each solution vector
  117: *          X(j) (the j-th column of the solution matrix X).
  118: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
  119: *          is an estimated upper bound for the magnitude of the largest
  120: *          element in (X(j) - XTRUE) divided by the magnitude of the
  121: *          largest element in X(j).
  122: *
  123: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  124: *          The componentwise relative backward error of each solution
  125: *          vector X(j) (i.e., the smallest relative change in any
  126: *          element of A or B that makes X(j) an exact solution).
  127: *
  128: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
  129: *
  130: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
  131: *
  132: *  INFO    (output) INTEGER
  133: *          = 0:  successful exit
  134: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  135: *          > 0:  if INFO = i, and i is
  136: *                <= N:  the leading minor of order i of A is
  137: *                       not positive definite, so the factorization
  138: *                       could not be completed, and the solution has not
  139: *                       been computed. RCOND = 0 is returned.
  140: *                = N+1: U is nonsingular, but RCOND is less than machine
  141: *                       precision, meaning that the matrix is singular
  142: *                       to working precision.  Nevertheless, the
  143: *                       solution and error bounds are computed because
  144: *                       there are a number of situations where the
  145: *                       computed solution can be more accurate than the
  146: *                       value of RCOND would suggest.
  147: *
  148: *  =====================================================================
  149: *
  150: *     .. Parameters ..
  151:       DOUBLE PRECISION   ZERO
  152:       PARAMETER          ( ZERO = 0.0D+0 )
  153: *     ..
  154: *     .. Local Scalars ..
  155:       LOGICAL            NOFACT
  156:       DOUBLE PRECISION   ANORM
  157: *     ..
  158: *     .. External Functions ..
  159:       LOGICAL            LSAME
  160:       DOUBLE PRECISION   DLAMCH, ZLANHT
  161:       EXTERNAL           LSAME, DLAMCH, ZLANHT
  162: *     ..
  163: *     .. External Subroutines ..
  164:       EXTERNAL           DCOPY, XERBLA, ZCOPY, ZLACPY, ZPTCON, ZPTRFS,
  165:      $                   ZPTTRF, ZPTTRS
  166: *     ..
  167: *     .. Intrinsic Functions ..
  168:       INTRINSIC          MAX
  169: *     ..
  170: *     .. Executable Statements ..
  171: *
  172: *     Test the input parameters.
  173: *
  174:       INFO = 0
  175:       NOFACT = LSAME( FACT, 'N' )
  176:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
  177:          INFO = -1
  178:       ELSE IF( N.LT.0 ) THEN
  179:          INFO = -2
  180:       ELSE IF( NRHS.LT.0 ) THEN
  181:          INFO = -3
  182:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  183:          INFO = -9
  184:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  185:          INFO = -11
  186:       END IF
  187:       IF( INFO.NE.0 ) THEN
  188:          CALL XERBLA( 'ZPTSVX', -INFO )
  189:          RETURN
  190:       END IF
  191: *
  192:       IF( NOFACT ) THEN
  193: *
  194: *        Compute the L*D*L' (or U'*D*U) factorization of A.
  195: *
  196:          CALL DCOPY( N, D, 1, DF, 1 )
  197:          IF( N.GT.1 )
  198:      $      CALL ZCOPY( N-1, E, 1, EF, 1 )
  199:          CALL ZPTTRF( N, DF, EF, INFO )
  200: *
  201: *        Return if INFO is non-zero.
  202: *
  203:          IF( INFO.GT.0 )THEN
  204:             RCOND = ZERO
  205:             RETURN
  206:          END IF
  207:       END IF
  208: *
  209: *     Compute the norm of the matrix A.
  210: *
  211:       ANORM = ZLANHT( '1', N, D, E )
  212: *
  213: *     Compute the reciprocal of the condition number of A.
  214: *
  215:       CALL ZPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
  216: *
  217: *     Compute the solution vectors X.
  218: *
  219:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
  220:       CALL ZPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO )
  221: *
  222: *     Use iterative refinement to improve the computed solutions and
  223: *     compute error bounds and backward error estimates for them.
  224: *
  225:       CALL ZPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
  226:      $             BERR, WORK, RWORK, INFO )
  227: *
  228: *     Set INFO = N+1 if the matrix is singular to working precision.
  229: *
  230:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
  231:      $   INFO = N + 1
  232: *
  233:       RETURN
  234: *
  235: *     End of ZPTSVX
  236: *
  237:       END