Return to zsysvx.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
    2:      $                   LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
    3:      $                   RWORK, INFO )
    4: *
    5: *  -- LAPACK driver 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: *     .. Scalar Arguments ..
   11:       CHARACTER          FACT, UPLO
   12:       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
   13:       DOUBLE PRECISION   RCOND
   14: *     ..
   15: *     .. Array Arguments ..
   16:       INTEGER            IPIV( * )
   17:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
   18:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
   19:      $                   WORK( * ), X( LDX, * )
   20: *     ..
   21: *
   22: *  Purpose
   23: *  =======
   24: *
   25: *  ZSYSVX uses the diagonal pivoting factorization to compute the
   26: *  solution to a complex system of linear equations A * X = B,
   27: *  where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
   28: *  matrices.
   29: *
   30: *  Error bounds on the solution and a condition estimate are also
   31: *  provided.
   32: *
   33: *  Description
   34: *  ===========
   35: *
   36: *  The following steps are performed:
   37: *
   38: *  1. If FACT = 'N', the diagonal pivoting method is used to factor A.
   39: *     The form of the factorization is
   40: *        A = U * D * U**T,  if UPLO = 'U', or
   41: *        A = L * D * L**T,  if UPLO = 'L',
   42: *     where U (or L) is a product of permutation and unit upper (lower)
   43: *     triangular matrices, and D is symmetric and block diagonal with
   44: *     1-by-1 and 2-by-2 diagonal blocks.
   45: *
   46: *  2. If some D(i,i)=0, so that D is exactly singular, then the routine
   47: *     returns with INFO = i. Otherwise, the factored form of A is used
   48: *     to estimate the condition number of the matrix A.  If the
   49: *     reciprocal of the condition number is less than machine precision,
   50: *     INFO = N+1 is returned as a warning, but the routine still goes on
   51: *     to solve for X and compute error bounds as described below.
   52: *
   53: *  3. The system of equations is solved for X using the factored form
   54: *     of A.
   55: *
   56: *  4. Iterative refinement is applied to improve the computed solution
   57: *     matrix and calculate error bounds and backward error estimates
   58: *     for it.
   59: *
   60: *  Arguments
   61: *  =========
   62: *
   63: *  FACT    (input) CHARACTER*1
   64: *          Specifies whether or not the factored form of A has been
   65: *          supplied on entry.
   66: *          = 'F':  On entry, AF and IPIV contain the factored form
   67: *                  of A.  A, AF and IPIV will not be modified.
   68: *          = 'N':  The matrix A will be copied to AF and factored.
   69: *
   70: *  UPLO    (input) CHARACTER*1
   71: *          = 'U':  Upper triangle of A is stored;
   72: *          = 'L':  Lower triangle of A is stored.
   73: *
   74: *  N       (input) INTEGER
   75: *          The number of linear equations, i.e., the order of the
   76: *          matrix A.  N >= 0.
   77: *
   78: *  NRHS    (input) INTEGER
   79: *          The number of right hand sides, i.e., the number of columns
   80: *          of the matrices B and X.  NRHS >= 0.
   81: *
   82: *  A       (input) COMPLEX*16 array, dimension (LDA,N)
   83: *          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
   84: *          upper triangular part of A contains the upper triangular part
   85: *          of the matrix A, and the strictly lower triangular part of A
   86: *          is not referenced.  If UPLO = 'L', the leading N-by-N lower
   87: *          triangular part of A contains the lower triangular part of
   88: *          the matrix A, and the strictly upper triangular part of A is
   89: *          not referenced.
   90: *
   91: *  LDA     (input) INTEGER
   92: *          The leading dimension of the array A.  LDA >= max(1,N).
   93: *
   94: *  AF      (input or output) COMPLEX*16 array, dimension (LDAF,N)
   95: *          If FACT = 'F', then AF is an input argument and on entry
   96: *          contains the block diagonal matrix D and the multipliers used
   97: *          to obtain the factor U or L from the factorization
   98: *          A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF.
   99: *
  100: *          If FACT = 'N', then AF is an output argument and on exit
  101: *          returns the block diagonal matrix D and the multipliers used
  102: *          to obtain the factor U or L from the factorization
  103: *          A = U*D*U**T or A = L*D*L**T.
  104: *
  105: *  LDAF    (input) INTEGER
  106: *          The leading dimension of the array AF.  LDAF >= max(1,N).
  107: *
  108: *  IPIV    (input or output) INTEGER array, dimension (N)
  109: *          If FACT = 'F', then IPIV is an input argument and on entry
  110: *          contains details of the interchanges and the block structure
  111: *          of D, as determined by ZSYTRF.
  112: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  113: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
  114: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
  115: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  116: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
  117: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
  118: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  119: *
  120: *          If FACT = 'N', then IPIV is an output argument and on exit
  121: *          contains details of the interchanges and the block structure
  122: *          of D, as determined by ZSYTRF.
  123: *
  124: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
  125: *          The N-by-NRHS right hand side matrix B.
  126: *
  127: *  LDB     (input) INTEGER
  128: *          The leading dimension of the array B.  LDB >= max(1,N).
  129: *
  130: *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
  131: *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
  132: *
  133: *  LDX     (input) INTEGER
  134: *          The leading dimension of the array X.  LDX >= max(1,N).
  135: *
  136: *  RCOND   (output) DOUBLE PRECISION
  137: *          The estimate of the reciprocal condition number of the matrix
  138: *          A.  If RCOND is less than the machine precision (in
  139: *          particular, if RCOND = 0), the matrix is singular to working
  140: *          precision.  This condition is indicated by a return code of
  141: *          INFO > 0.
  142: *
  143: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  144: *          The estimated forward error bound for each solution vector
  145: *          X(j) (the j-th column of the solution matrix X).
  146: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
  147: *          is an estimated upper bound for the magnitude of the largest
  148: *          element in (X(j) - XTRUE) divided by the magnitude of the
  149: *          largest element in X(j).  The estimate is as reliable as
  150: *          the estimate for RCOND, and is almost always a slight
  151: *          overestimate of the true error.
  152: *
  153: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  154: *          The componentwise relative backward error of each solution
  155: *          vector X(j) (i.e., the smallest relative change in
  156: *          any element of A or B that makes X(j) an exact solution).
  157: *
  158: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
  159: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  160: *
  161: *  LWORK   (input) INTEGER
  162: *          The length of WORK.  LWORK >= max(1,2*N), and for best
  163: *          performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
  164: *          NB is the optimal blocksize for ZSYTRF.
  165: *
  166: *          If LWORK = -1, then a workspace query is assumed; the routine
  167: *          only calculates the optimal size of the WORK array, returns
  168: *          this value as the first entry of the WORK array, and no error
  169: *          message related to LWORK is issued by XERBLA.
  170: *
  171: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
  172: *
  173: *  INFO    (output) INTEGER
  174: *          = 0: successful exit
  175: *          < 0: if INFO = -i, the i-th argument had an illegal value
  176: *          > 0: if INFO = i, and i is
  177: *                <= N:  D(i,i) is exactly zero.  The factorization
  178: *                       has been completed but the factor D is exactly
  179: *                       singular, so the solution and error bounds could
  180: *                       not be computed. RCOND = 0 is returned.
  181: *                = N+1: D is nonsingular, but RCOND is less than machine
  182: *                       precision, meaning that the matrix is singular
  183: *                       to working precision.  Nevertheless, the
  184: *                       solution and error bounds are computed because
  185: *                       there are a number of situations where the
  186: *                       computed solution can be more accurate than the
  187: *                       value of RCOND would suggest.
  188: *
  189: *  =====================================================================
  190: *
  191: *     .. Parameters ..
  192:       DOUBLE PRECISION   ZERO
  193:       PARAMETER          ( ZERO = 0.0D+0 )
  194: *     ..
  195: *     .. Local Scalars ..
  196:       LOGICAL            LQUERY, NOFACT
  197:       INTEGER            LWKOPT, NB
  198:       DOUBLE PRECISION   ANORM
  199: *     ..
  200: *     .. External Functions ..
  201:       LOGICAL            LSAME
  202:       INTEGER            ILAENV
  203:       DOUBLE PRECISION   DLAMCH, ZLANSY
  204:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANSY
  205: *     ..
  206: *     .. External Subroutines ..
  207:       EXTERNAL           XERBLA, ZLACPY, ZSYCON, ZSYRFS, ZSYTRF, ZSYTRS
  208: *     ..
  209: *     .. Intrinsic Functions ..
  210:       INTRINSIC          MAX
  211: *     ..
  212: *     .. Executable Statements ..
  213: *
  214: *     Test the input parameters.
  215: *
  216:       INFO = 0
  217:       NOFACT = LSAME( FACT, 'N' )
  218:       LQUERY = ( LWORK.EQ.-1 )
  219:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
  220:          INFO = -1
  221:       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
  222:      $          THEN
  223:          INFO = -2
  224:       ELSE IF( N.LT.0 ) THEN
  225:          INFO = -3
  226:       ELSE IF( NRHS.LT.0 ) THEN
  227:          INFO = -4
  228:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  229:          INFO = -6
  230:       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
  231:          INFO = -8
  232:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  233:          INFO = -11
  234:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  235:          INFO = -13
  236:       ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
  237:          INFO = -18
  238:       END IF
  239: *
  240:       IF( INFO.EQ.0 ) THEN
  241:          LWKOPT = MAX( 1, 2*N )
  242:          IF( NOFACT ) THEN
  243:             NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 )
  244:             LWKOPT = MAX( LWKOPT, N*NB )
  245:          END IF
  246:          WORK( 1 ) = LWKOPT
  247:       END IF
  248: *
  249:       IF( INFO.NE.0 ) THEN
  250:          CALL XERBLA( 'ZSYSVX', -INFO )
  251:          RETURN
  252:       ELSE IF( LQUERY ) THEN
  253:          RETURN
  254:       END IF
  255: *
  256:       IF( NOFACT ) THEN
  257: *
  258: *        Compute the factorization A = U*D*U' or A = L*D*L'.
  259: *
  260:          CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
  261:          CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
  262: *
  263: *        Return if INFO is non-zero.
  264: *
  265:          IF( INFO.GT.0 )THEN
  266:             RCOND = ZERO
  267:             RETURN
  268:          END IF
  269:       END IF
  270: *
  271: *     Compute the norm of the matrix A.
  272: *
  273:       ANORM = ZLANSY( 'I', UPLO, N, A, LDA, RWORK )
  274: *
  275: *     Compute the reciprocal of the condition number of A.
  276: *
  277:       CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
  278: *
  279: *     Compute the solution vectors X.
  280: *
  281:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
  282:       CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
  283: *
  284: *     Use iterative refinement to improve the computed solutions and
  285: *     compute error bounds and backward error estimates for them.
  286: *
  287:       CALL ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
  288:      $             LDX, FERR, BERR, WORK, RWORK, INFO )
  289: *
  290: *     Set INFO = N+1 if the matrix is singular to working precision.
  291: *
  292:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
  293:      $   INFO = N + 1
  294: *
  295:       WORK( 1 ) = LWKOPT
  296: *
  297:       RETURN
  298: *
  299: *     End of ZSYSVX
  300: *
  301:       END