Return to zsysvx.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief <b> ZSYSVX computes the solution to system of linear equations A * X = B for SY 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 ZSYSVX + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsysvx.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsysvx.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsysvx.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
   22: *                          LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
   23: *                          RWORK, INFO )
   24: *
   25: *       .. Scalar Arguments ..
   26: *       CHARACTER          FACT, UPLO
   27: *       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
   28: *       DOUBLE PRECISION   RCOND
   29: *       ..
   30: *       .. Array Arguments ..
   31: *       INTEGER            IPIV( * )
   32: *       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
   33: *       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
   34: *      $                   WORK( * ), X( LDX, * )
   35: *       ..
   36: *
   37: *
   38: *> \par Purpose:
   39: *  =============
   40: *>
   41: *> \verbatim
   42: *>
   43: *> ZSYSVX uses the diagonal pivoting factorization to compute the
   44: *> solution to a complex system of linear equations A * X = B,
   45: *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
   46: *> matrices.
   47: *>
   48: *> Error bounds on the solution and a condition estimate are also
   49: *> provided.
   50: *> \endverbatim
   51: *
   52: *> \par Description:
   53: *  =================
   54: *>
   55: *> \verbatim
   56: *>
   57: *> The following steps are performed:
   58: *>
   59: *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
   60: *>    The form of the factorization is
   61: *>       A = U * D * U**T,  if UPLO = 'U', or
   62: *>       A = L * D * L**T,  if UPLO = 'L',
   63: *>    where U (or L) is a product of permutation and unit upper (lower)
   64: *>    triangular matrices, and D is symmetric and block diagonal with
   65: *>    1-by-1 and 2-by-2 diagonal blocks.
   66: *>
   67: *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
   68: *>    returns with INFO = i. Otherwise, the factored form of A is used
   69: *>    to estimate the condition number of the matrix A.  If the
   70: *>    reciprocal of the condition number is less than machine precision,
   71: *>    INFO = N+1 is returned as a warning, but the routine still goes on
   72: *>    to solve for X and compute error bounds as described below.
   73: *>
   74: *> 3. The system of equations is solved for X using the factored form
   75: *>    of A.
   76: *>
   77: *> 4. Iterative refinement is applied to improve the computed solution
   78: *>    matrix and calculate error bounds and backward error estimates
   79: *>    for it.
   80: *> \endverbatim
   81: *
   82: *  Arguments:
   83: *  ==========
   84: *
   85: *> \param[in] FACT
   86: *> \verbatim
   87: *>          FACT is CHARACTER*1
   88: *>          Specifies whether or not the factored form of A has been
   89: *>          supplied on entry.
   90: *>          = 'F':  On entry, AF and IPIV contain the factored form
   91: *>                  of A.  A, AF and IPIV will not be modified.
   92: *>          = 'N':  The matrix A will be copied to AF and factored.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] UPLO
   96: *> \verbatim
   97: *>          UPLO is CHARACTER*1
   98: *>          = 'U':  Upper triangle of A is stored;
   99: *>          = 'L':  Lower triangle of A is stored.
  100: *> \endverbatim
  101: *>
  102: *> \param[in] N
  103: *> \verbatim
  104: *>          N is INTEGER
  105: *>          The number of linear equations, i.e., the order of the
  106: *>          matrix A.  N >= 0.
  107: *> \endverbatim
  108: *>
  109: *> \param[in] NRHS
  110: *> \verbatim
  111: *>          NRHS is INTEGER
  112: *>          The number of right hand sides, i.e., the number of columns
  113: *>          of the matrices B and X.  NRHS >= 0.
  114: *> \endverbatim
  115: *>
  116: *> \param[in] A
  117: *> \verbatim
  118: *>          A is COMPLEX*16 array, dimension (LDA,N)
  119: *>          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
  120: *>          upper triangular part of A contains the upper triangular part
  121: *>          of the matrix A, and the strictly lower triangular part of A
  122: *>          is not referenced.  If UPLO = 'L', the leading N-by-N lower
  123: *>          triangular part of A contains the lower triangular part of
  124: *>          the matrix A, and the strictly upper triangular part of A is
  125: *>          not referenced.
  126: *> \endverbatim
  127: *>
  128: *> \param[in] LDA
  129: *> \verbatim
  130: *>          LDA is INTEGER
  131: *>          The leading dimension of the array A.  LDA >= max(1,N).
  132: *> \endverbatim
  133: *>
  134: *> \param[in,out] AF
  135: *> \verbatim
  136: *>          AF is COMPLEX*16 array, dimension (LDAF,N)
  137: *>          If FACT = 'F', then AF is an input argument and on entry
  138: *>          contains the block diagonal matrix D and the multipliers used
  139: *>          to obtain the factor U or L from the factorization
  140: *>          A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF.
  141: *>
  142: *>          If FACT = 'N', then AF is an output argument and on exit
  143: *>          returns the block diagonal matrix D and the multipliers used
  144: *>          to obtain the factor U or L from the factorization
  145: *>          A = U*D*U**T or A = L*D*L**T.
  146: *> \endverbatim
  147: *>
  148: *> \param[in] LDAF
  149: *> \verbatim
  150: *>          LDAF is INTEGER
  151: *>          The leading dimension of the array AF.  LDAF >= max(1,N).
  152: *> \endverbatim
  153: *>
  154: *> \param[in,out] IPIV
  155: *> \verbatim
  156: *>          IPIV is INTEGER array, dimension (N)
  157: *>          If FACT = 'F', then IPIV is an input argument and on entry
  158: *>          contains details of the interchanges and the block structure
  159: *>          of D, as determined by ZSYTRF.
  160: *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  161: *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
  162: *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
  163: *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  164: *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
  165: *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
  166: *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  167: *>
  168: *>          If FACT = 'N', then IPIV is an output argument and on exit
  169: *>          contains details of the interchanges and the block structure
  170: *>          of D, as determined by ZSYTRF.
  171: *> \endverbatim
  172: *>
  173: *> \param[in] B
  174: *> \verbatim
  175: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  176: *>          The N-by-NRHS right hand side matrix B.
  177: *> \endverbatim
  178: *>
  179: *> \param[in] LDB
  180: *> \verbatim
  181: *>          LDB is INTEGER
  182: *>          The leading dimension of the array B.  LDB >= max(1,N).
  183: *> \endverbatim
  184: *>
  185: *> \param[out] X
  186: *> \verbatim
  187: *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
  188: *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
  189: *> \endverbatim
  190: *>
  191: *> \param[in] LDX
  192: *> \verbatim
  193: *>          LDX is INTEGER
  194: *>          The leading dimension of the array X.  LDX >= max(1,N).
  195: *> \endverbatim
  196: *>
  197: *> \param[out] RCOND
  198: *> \verbatim
  199: *>          RCOND is DOUBLE PRECISION
  200: *>          The estimate of the reciprocal condition number of the matrix
  201: *>          A.  If RCOND is less than the machine precision (in
  202: *>          particular, if RCOND = 0), the matrix is singular to working
  203: *>          precision.  This condition is indicated by a return code of
  204: *>          INFO > 0.
  205: *> \endverbatim
  206: *>
  207: *> \param[out] FERR
  208: *> \verbatim
  209: *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
  210: *>          The estimated forward error bound for each solution vector
  211: *>          X(j) (the j-th column of the solution matrix X).
  212: *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
  213: *>          is an estimated upper bound for the magnitude of the largest
  214: *>          element in (X(j) - XTRUE) divided by the magnitude of the
  215: *>          largest element in X(j).  The estimate is as reliable as
  216: *>          the estimate for RCOND, and is almost always a slight
  217: *>          overestimate of the true error.
  218: *> \endverbatim
  219: *>
  220: *> \param[out] BERR
  221: *> \verbatim
  222: *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
  223: *>          The componentwise relative backward error of each solution
  224: *>          vector X(j) (i.e., the smallest relative change in
  225: *>          any element of A or B that makes X(j) an exact solution).
  226: *> \endverbatim
  227: *>
  228: *> \param[out] WORK
  229: *> \verbatim
  230: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
  231: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  232: *> \endverbatim
  233: *>
  234: *> \param[in] LWORK
  235: *> \verbatim
  236: *>          LWORK is INTEGER
  237: *>          The length of WORK.  LWORK >= max(1,2*N), and for best
  238: *>          performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
  239: *>          NB is the optimal blocksize for ZSYTRF.
  240: *>
  241: *>          If LWORK = -1, then a workspace query is assumed; the routine
  242: *>          only calculates the optimal size of the WORK array, returns
  243: *>          this value as the first entry of the WORK array, and no error
  244: *>          message related to LWORK is issued by XERBLA.
  245: *> \endverbatim
  246: *>
  247: *> \param[out] RWORK
  248: *> \verbatim
  249: *>          RWORK is DOUBLE PRECISION array, dimension (N)
  250: *> \endverbatim
  251: *>
  252: *> \param[out] INFO
  253: *> \verbatim
  254: *>          INFO is INTEGER
  255: *>          = 0: successful exit
  256: *>          < 0: if INFO = -i, the i-th argument had an illegal value
  257: *>          > 0: if INFO = i, and i is
  258: *>                <= N:  D(i,i) is exactly zero.  The factorization
  259: *>                       has been completed but the factor D is exactly
  260: *>                       singular, so the solution and error bounds could
  261: *>                       not be computed. RCOND = 0 is returned.
  262: *>                = N+1: D is nonsingular, but RCOND is less than machine
  263: *>                       precision, meaning that the matrix is singular
  264: *>                       to working precision.  Nevertheless, the
  265: *>                       solution and error bounds are computed because
  266: *>                       there are a number of situations where the
  267: *>                       computed solution can be more accurate than the
  268: *>                       value of RCOND would suggest.
  269: *> \endverbatim
  270: *
  271: *  Authors:
  272: *  ========
  273: *
  274: *> \author Univ. of Tennessee
  275: *> \author Univ. of California Berkeley
  276: *> \author Univ. of Colorado Denver
  277: *> \author NAG Ltd.
  278: *
  279: *> \ingroup complex16SYsolve
  280: *
  281: *  =====================================================================
  282:       SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
  283:      $                   LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
  284:      $                   RWORK, INFO )
  285: *
  286: *  -- LAPACK driver routine --
  287: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  288: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  289: *
  290: *     .. Scalar Arguments ..
  291:       CHARACTER          FACT, UPLO
  292:       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
  293:       DOUBLE PRECISION   RCOND
  294: *     ..
  295: *     .. Array Arguments ..
  296:       INTEGER            IPIV( * )
  297:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
  298:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  299:      $                   WORK( * ), X( LDX, * )
  300: *     ..
  301: *
  302: *  =====================================================================
  303: *
  304: *     .. Parameters ..
  305:       DOUBLE PRECISION   ZERO
  306:       PARAMETER          ( ZERO = 0.0D+0 )
  307: *     ..
  308: *     .. Local Scalars ..
  309:       LOGICAL            LQUERY, NOFACT
  310:       INTEGER            LWKOPT, NB
  311:       DOUBLE PRECISION   ANORM
  312: *     ..
  313: *     .. External Functions ..
  314:       LOGICAL            LSAME
  315:       INTEGER            ILAENV
  316:       DOUBLE PRECISION   DLAMCH, ZLANSY
  317:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANSY
  318: *     ..
  319: *     .. External Subroutines ..
  320:       EXTERNAL           XERBLA, ZLACPY, ZSYCON, ZSYRFS, ZSYTRF, ZSYTRS
  321: *     ..
  322: *     .. Intrinsic Functions ..
  323:       INTRINSIC          MAX
  324: *     ..
  325: *     .. Executable Statements ..
  326: *
  327: *     Test the input parameters.
  328: *
  329:       INFO = 0
  330:       NOFACT = LSAME( FACT, 'N' )
  331:       LQUERY = ( LWORK.EQ.-1 )
  332:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
  333:          INFO = -1
  334:       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
  335:      $          THEN
  336:          INFO = -2
  337:       ELSE IF( N.LT.0 ) THEN
  338:          INFO = -3
  339:       ELSE IF( NRHS.LT.0 ) THEN
  340:          INFO = -4
  341:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  342:          INFO = -6
  343:       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
  344:          INFO = -8
  345:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  346:          INFO = -11
  347:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  348:          INFO = -13
  349:       ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
  350:          INFO = -18
  351:       END IF
  352: *
  353:       IF( INFO.EQ.0 ) THEN
  354:          LWKOPT = MAX( 1, 2*N )
  355:          IF( NOFACT ) THEN
  356:             NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 )
  357:             LWKOPT = MAX( LWKOPT, N*NB )
  358:          END IF
  359:          WORK( 1 ) = LWKOPT
  360:       END IF
  361: *
  362:       IF( INFO.NE.0 ) THEN
  363:          CALL XERBLA( 'ZSYSVX', -INFO )
  364:          RETURN
  365:       ELSE IF( LQUERY ) THEN
  366:          RETURN
  367:       END IF
  368: *
  369:       IF( NOFACT ) THEN
  370: *
  371: *        Compute the factorization A = U*D*U**T or A = L*D*L**T.
  372: *
  373:          CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
  374:          CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
  375: *
  376: *        Return if INFO is non-zero.
  377: *
  378:          IF( INFO.GT.0 )THEN
  379:             RCOND = ZERO
  380:             RETURN
  381:          END IF
  382:       END IF
  383: *
  384: *     Compute the norm of the matrix A.
  385: *
  386:       ANORM = ZLANSY( 'I', UPLO, N, A, LDA, RWORK )
  387: *
  388: *     Compute the reciprocal of the condition number of A.
  389: *
  390:       CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
  391: *
  392: *     Compute the solution vectors X.
  393: *
  394:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
  395:       CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
  396: *
  397: *     Use iterative refinement to improve the computed solutions and
  398: *     compute error bounds and backward error estimates for them.
  399: *
  400:       CALL ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
  401:      $             LDX, FERR, BERR, WORK, RWORK, INFO )
  402: *
  403: *     Set INFO = N+1 if the matrix is singular to working precision.
  404: *
  405:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
  406:      $   INFO = N + 1
  407: *
  408:       WORK( 1 ) = LWKOPT
  409: *
  410:       RETURN
  411: *
  412: *     End of ZSYSVX
  413: *
  414:       END