Return to zptsvx.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief <b> ZPTSVX computes the solution to system of linear equations A * X = B for PT matrices</b>
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZPTSVX + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zptsvx.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zptsvx.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zptsvx.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
   22: *                          RCOND, FERR, BERR, WORK, RWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          FACT
   26: *       INTEGER            INFO, LDB, LDX, N, NRHS
   27: *       DOUBLE PRECISION   RCOND
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       DOUBLE PRECISION   BERR( * ), D( * ), DF( * ), FERR( * ),
   31: *      $                   RWORK( * )
   32: *       COMPLEX*16         B( LDB, * ), E( * ), EF( * ), WORK( * ),
   33: *      $                   X( LDX, * )
   34: *       ..
   35: *
   36: *
   37: *> \par Purpose:
   38: *  =============
   39: *>
   40: *> \verbatim
   41: *>
   42: *> ZPTSVX uses the factorization A = L*D*L**H to compute the solution
   43: *> to a complex system of linear equations A*X = B, where A is an
   44: *> N-by-N Hermitian positive definite tridiagonal matrix and X and B
   45: *> are N-by-NRHS matrices.
   46: *>
   47: *> Error bounds on the solution and a condition estimate are also
   48: *> provided.
   49: *> \endverbatim
   50: *
   51: *> \par Description:
   52: *  =================
   53: *>
   54: *> \verbatim
   55: *>
   56: *> The following steps are performed:
   57: *>
   58: *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**H, where L
   59: *>    is a unit lower bidiagonal matrix and D is diagonal.  The
   60: *>    factorization can also be regarded as having the form
   61: *>    A = U**H*D*U.
   62: *>
   63: *> 2. If the leading i-by-i principal minor is not positive definite,
   64: *>    then the routine returns with INFO = i. Otherwise, the factored
   65: *>    form of A is used to estimate the condition number of the matrix
   66: *>    A.  If the reciprocal of the condition number is less than machine
   67: *>    precision, INFO = N+1 is returned as a warning, but the routine
   68: *>    still goes on to solve for X and compute error bounds as
   69: *>    described below.
   70: *>
   71: *> 3. The system of equations is solved for X using the factored form
   72: *>    of A.
   73: *>
   74: *> 4. Iterative refinement is applied to improve the computed solution
   75: *>    matrix and calculate error bounds and backward error estimates
   76: *>    for it.
   77: *> \endverbatim
   78: *
   79: *  Arguments:
   80: *  ==========
   81: *
   82: *> \param[in] FACT
   83: *> \verbatim
   84: *>          FACT is CHARACTER*1
   85: *>          Specifies whether or not the factored form of the matrix
   86: *>          A is supplied on entry.
   87: *>          = 'F':  On entry, DF and EF contain the factored form of A.
   88: *>                  D, E, DF, and EF will not be modified.
   89: *>          = 'N':  The matrix A will be copied to DF and EF and
   90: *>                  factored.
   91: *> \endverbatim
   92: *>
   93: *> \param[in] N
   94: *> \verbatim
   95: *>          N is INTEGER
   96: *>          The order of the matrix A.  N >= 0.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] NRHS
  100: *> \verbatim
  101: *>          NRHS is INTEGER
  102: *>          The number of right hand sides, i.e., the number of columns
  103: *>          of the matrices B and X.  NRHS >= 0.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] D
  107: *> \verbatim
  108: *>          D is DOUBLE PRECISION array, dimension (N)
  109: *>          The n diagonal elements of the tridiagonal matrix A.
  110: *> \endverbatim
  111: *>
  112: *> \param[in] E
  113: *> \verbatim
  114: *>          E is COMPLEX*16 array, dimension (N-1)
  115: *>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
  116: *> \endverbatim
  117: *>
  118: *> \param[in,out] DF
  119: *> \verbatim
  120: *>          DF is DOUBLE PRECISION array, dimension (N)
  121: *>          If FACT = 'F', then DF is an input argument and on entry
  122: *>          contains the n diagonal elements of the diagonal matrix D
  123: *>          from the L*D*L**H factorization of A.
  124: *>          If FACT = 'N', then DF is an output argument and on exit
  125: *>          contains the n diagonal elements of the diagonal matrix D
  126: *>          from the L*D*L**H factorization of A.
  127: *> \endverbatim
  128: *>
  129: *> \param[in,out] EF
  130: *> \verbatim
  131: *>          EF is COMPLEX*16 array, dimension (N-1)
  132: *>          If FACT = 'F', then EF is an input argument and on entry
  133: *>          contains the (n-1) subdiagonal elements of the unit
  134: *>          bidiagonal factor L from the L*D*L**H factorization of A.
  135: *>          If FACT = 'N', then EF is an output argument and on exit
  136: *>          contains the (n-1) subdiagonal elements of the unit
  137: *>          bidiagonal factor L from the L*D*L**H factorization of A.
  138: *> \endverbatim
  139: *>
  140: *> \param[in] B
  141: *> \verbatim
  142: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  143: *>          The N-by-NRHS right hand side matrix B.
  144: *> \endverbatim
  145: *>
  146: *> \param[in] LDB
  147: *> \verbatim
  148: *>          LDB is INTEGER
  149: *>          The leading dimension of the array B.  LDB >= max(1,N).
  150: *> \endverbatim
  151: *>
  152: *> \param[out] X
  153: *> \verbatim
  154: *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
  155: *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
  156: *> \endverbatim
  157: *>
  158: *> \param[in] LDX
  159: *> \verbatim
  160: *>          LDX is INTEGER
  161: *>          The leading dimension of the array X.  LDX >= max(1,N).
  162: *> \endverbatim
  163: *>
  164: *> \param[out] RCOND
  165: *> \verbatim
  166: *>          RCOND is DOUBLE PRECISION
  167: *>          The reciprocal condition number of the matrix A.  If RCOND
  168: *>          is less than the machine precision (in particular, if
  169: *>          RCOND = 0), the matrix is singular to working precision.
  170: *>          This condition is indicated by a return code of INFO > 0.
  171: *> \endverbatim
  172: *>
  173: *> \param[out] FERR
  174: *> \verbatim
  175: *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
  176: *>          The forward error bound for each solution vector
  177: *>          X(j) (the j-th column of the solution matrix X).
  178: *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
  179: *>          is an estimated upper bound for the magnitude of the largest
  180: *>          element in (X(j) - XTRUE) divided by the magnitude of the
  181: *>          largest element in X(j).
  182: *> \endverbatim
  183: *>
  184: *> \param[out] BERR
  185: *> \verbatim
  186: *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
  187: *>          The componentwise relative backward error of each solution
  188: *>          vector X(j) (i.e., the smallest relative change in any
  189: *>          element of A or B that makes X(j) an exact solution).
  190: *> \endverbatim
  191: *>
  192: *> \param[out] WORK
  193: *> \verbatim
  194: *>          WORK is COMPLEX*16 array, dimension (N)
  195: *> \endverbatim
  196: *>
  197: *> \param[out] RWORK
  198: *> \verbatim
  199: *>          RWORK is DOUBLE PRECISION array, dimension (N)
  200: *> \endverbatim
  201: *>
  202: *> \param[out] INFO
  203: *> \verbatim
  204: *>          INFO is INTEGER
  205: *>          = 0:  successful exit
  206: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  207: *>          > 0:  if INFO = i, and i is
  208: *>                <= N:  the leading minor of order i of A is
  209: *>                       not positive definite, so the factorization
  210: *>                       could not be completed, and the solution has not
  211: *>                       been computed. RCOND = 0 is returned.
  212: *>                = N+1: U is nonsingular, but RCOND is less than machine
  213: *>                       precision, meaning that the matrix is singular
  214: *>                       to working precision.  Nevertheless, the
  215: *>                       solution and error bounds are computed because
  216: *>                       there are a number of situations where the
  217: *>                       computed solution can be more accurate than the
  218: *>                       value of RCOND would suggest.
  219: *> \endverbatim
  220: *
  221: *  Authors:
  222: *  ========
  223: *
  224: *> \author Univ. of Tennessee
  225: *> \author Univ. of California Berkeley
  226: *> \author Univ. of Colorado Denver
  227: *> \author NAG Ltd.
  228: *
  229: *> \date December 2016
  230: *
  231: *> \ingroup complex16PTsolve
  232: *
  233: *  =====================================================================
  234:       SUBROUTINE ZPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
  235:      $                   RCOND, FERR, BERR, WORK, RWORK, INFO )
  236: *
  237: *  -- LAPACK driver routine (version 3.7.0) --
  238: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  239: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  240: *     December 2016
  241: *
  242: *     .. Scalar Arguments ..
  243:       CHARACTER          FACT
  244:       INTEGER            INFO, LDB, LDX, N, NRHS
  245:       DOUBLE PRECISION   RCOND
  246: *     ..
  247: *     .. Array Arguments ..
  248:       DOUBLE PRECISION   BERR( * ), D( * ), DF( * ), FERR( * ),
  249:      $                   RWORK( * )
  250:       COMPLEX*16         B( LDB, * ), E( * ), EF( * ), WORK( * ),
  251:      $                   X( LDX, * )
  252: *     ..
  253: *
  254: *  =====================================================================
  255: *
  256: *     .. Parameters ..
  257:       DOUBLE PRECISION   ZERO
  258:       PARAMETER          ( ZERO = 0.0D+0 )
  259: *     ..
  260: *     .. Local Scalars ..
  261:       LOGICAL            NOFACT
  262:       DOUBLE PRECISION   ANORM
  263: *     ..
  264: *     .. External Functions ..
  265:       LOGICAL            LSAME
  266:       DOUBLE PRECISION   DLAMCH, ZLANHT
  267:       EXTERNAL           LSAME, DLAMCH, ZLANHT
  268: *     ..
  269: *     .. External Subroutines ..
  270:       EXTERNAL           DCOPY, XERBLA, ZCOPY, ZLACPY, ZPTCON, ZPTRFS,
  271:      $                   ZPTTRF, ZPTTRS
  272: *     ..
  273: *     .. Intrinsic Functions ..
  274:       INTRINSIC          MAX
  275: *     ..
  276: *     .. Executable Statements ..
  277: *
  278: *     Test the input parameters.
  279: *
  280:       INFO = 0
  281:       NOFACT = LSAME( FACT, 'N' )
  282:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
  283:          INFO = -1
  284:       ELSE IF( N.LT.0 ) THEN
  285:          INFO = -2
  286:       ELSE IF( NRHS.LT.0 ) THEN
  287:          INFO = -3
  288:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  289:          INFO = -9
  290:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  291:          INFO = -11
  292:       END IF
  293:       IF( INFO.NE.0 ) THEN
  294:          CALL XERBLA( 'ZPTSVX', -INFO )
  295:          RETURN
  296:       END IF
  297: *
  298:       IF( NOFACT ) THEN
  299: *
  300: *        Compute the L*D*L**H (or U**H*D*U) factorization of A.
  301: *
  302:          CALL DCOPY( N, D, 1, DF, 1 )
  303:          IF( N.GT.1 )
  304:      $      CALL ZCOPY( N-1, E, 1, EF, 1 )
  305:          CALL ZPTTRF( N, DF, EF, INFO )
  306: *
  307: *        Return if INFO is non-zero.
  308: *
  309:          IF( INFO.GT.0 )THEN
  310:             RCOND = ZERO
  311:             RETURN
  312:          END IF
  313:       END IF
  314: *
  315: *     Compute the norm of the matrix A.
  316: *
  317:       ANORM = ZLANHT( '1', N, D, E )
  318: *
  319: *     Compute the reciprocal of the condition number of A.
  320: *
  321:       CALL ZPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
  322: *
  323: *     Compute the solution vectors X.
  324: *
  325:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
  326:       CALL ZPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO )
  327: *
  328: *     Use iterative refinement to improve the computed solutions and
  329: *     compute error bounds and backward error estimates for them.
  330: *
  331:       CALL ZPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
  332:      $             BERR, WORK, RWORK, INFO )
  333: *
  334: *     Set INFO = N+1 if the matrix is singular to working precision.
  335: *
  336:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
  337:      $   INFO = N + 1
  338: *
  339:       RETURN
  340: *
  341: *     End of ZPTSVX
  342: *
  343:       END