LAPACK Archives

[Lapack] Proposed interfaces for extra-precise refinement.

Appended are the proposed interfaces for extra-precise
refinement in single precision.  Other precisions and
matrix types are similar.  The interfaces are, alas, 
complicated.  Comments appreciated.

  SGESVXX : Expert driver that factors, solves, and refines.

  SGERFSX : Refinement routine itself.

Still in flux:
  1) How to query for the workspace size.  The complex 
     version has integer, real, and complex workspace;
     adding three LWORK-like parameters feels a bit
     strange.

  2) The actual workspace size.  If we remove the
     cautious scaling from the condition number 
     estimator, we need either N or 2*N less
     workspace.

Jason

-------------- next part --------------
      SUBROUTINE SGESVXX(FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
     $     EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS,
     $     ERR_BNDS, NPARAMS, PARAMS, WORK, IWORK, INFO)
*
*     -- LAPACK driver routine (version 3.X) --
*     Copyright (c) 2005-2006, Regents of the University of California,
*     Univ. of Tennessee, NAG Ltd., Courant Institute, Argonne National
*     Lab, and Rice University
*     See COPYING file for license.
*
*     .. Scalar Arguments ..
      CHARACTER          EQUED, FACT, TRANS
      INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
     $                   N_ERR_BNDS
      REAL               RCOND, RPVGRW
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IWORK( * )
      REAL               A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
     $                   R( * ), C( * ), X( LDX, * ), PARAMS(*),
     $                   ERR_BNDS(NRHS, 2, N_ERR_BNDS), WORK ( * )

*
*     ..
*
*     Purpose
*     =======
*
*     SGESVXX uses the LU factorization to compute the solution to a
*     real system of linear equations  A * X = B,  where A is an
*     N-by-N matrix and X and B are N-by-NRHS matrix.
*
*     If iterative refinement is performed, both normwise and maximum
*     componentwise error bound is returned, as well as a condition
*     estimate.  SGESVXX will return a solution with O(eps) error (where
*     eps is the working machine precision) unless the matrix is very
*     ill-conditioned.  Relevant condition numbers also are calculated
*     and returned.
*
*     Description
*     ===========
*
*     The following steps are performed:
*
*     1. If FACT = 'E', real scaling factors are computed to equilibrate
*     the system:
*
*       TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
*       TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
*       TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
*
*     Whether or not the system will be equilibrated depends on the
*     scaling of the matrix A, but if equilibration is used, A is
*     overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
*     or diag(C)*B (if TRANS = 'T' or 'C').
*
*     2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
*     matrix A (after equilibration if FACT = 'E') as
*
*       A = P * L * U,
*
*     where P is a permutation matrix, L is a unit lower triangular
*     matrix, and U is upper triangular.
*
*     3. If some U(i,i)=0, so that U is exactly singular, then the routine
*     returns with INFO = i. Otherwise, the factored form of A is used
*     to estimate the condition number of the matrix A.  If the
*     reciprocal of the condition number is less than machine precision,
*     INFO = N+1 is returned as a warning, but the routine still goes on
*     to solve for X and compute error bounds as described below.
*
*     4. The system of equations is solved for X using the factored form
*     of A.
*
*     5. PARAMS(ITREF_PARAM) specifies whether we are performing
*     iterative refinement.  Refinement calculates the residual to
*     at least twice the working precision.
*
*     6. If equilibration was used, the matrix X is premultiplied by
*     diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
*     that it solves the original system before equilibration.
*
*     Arguments
*     =========
*
*     Some parameters are bundled in the PARAMS array.  These settings
*     affect other behavior, so we describe them first.  If the defaults
*     are acceptable, users can pass NPARAMS = 0 and the code will never
*     access the PARAMS argument.
*
*     NPARAMS (input) INTEGER
*     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
*     PARAMS array is never referenced.
*
*     PARAMS  (input / output) REAL array, dimension NPARAMS
*     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
*     that entry will be filled with default value used for that
*     parameter.  Only positions up to NPARAMS are accessed; defaults
*     are used for unavailable parameters.  PARAMETER declarations for
*     the index parameters below are specified in sgerfsx_params.fh.
*
*       PARAMS(ITREF_PARAM = 1) : Whether to perform iterative 
*            refinement or not.
*         Default: 1.0
*            = 0.0 : No refinement is perfomed, and no error bounds are
*                    computed.  Error arguments are not touched.
*            = 1.0 : Use the double-precision refinement algorithm, 
*                    possibly with doubled-single computations if the 
*                    compilation environment does not support DOUBLE 
*                    PRECISION.
*              (other values are reserved for future use)
*
*       PARAMS(RCONDTHRESH_PARAM = 2) : Reciprocal condition number 
*            threshold where the error estimates are no longer 
*            considered trustworthy.
*            Change with extreme caution.
*         Default: SQRT(N) * SLAMCH('Epsilon')
*
*       PARAMS(ITHRESH_PARAM = 3) : Total number of residual computations
*            allowed for refinement.
*         Default: 10 for double-precision, 5 for single-precision
*         Aggressive: Set to 100 for our "aggresive" settings.
*
*       PARAMS(COMPONENTWISE_PARAM = 4) : Flag determining if the code
*            will attempt to find a solution with small componentwise
*            relative error in the double-precision algorithm.  Non-zero
*            is true, 0.0 is false.
*         Default: 1.0 (attempt componentwise convergenge)
*
*       PARAMS(RTHRESH_PARAM = 5) : Ratio of consecutive "step sizes"
*            required to continue relative normwise or componentwise
*            refinement.  If the norm of the previous step (dx or dz) is
*            SPREV and the norm of the current step is SCURR, then SCURR
*            / SPREV must be .LE. PARAMS(RTHRESH_PARAM) to continue.
*         Default: 0.5
*         Aggressive: Set to 0.9 for our "aggresive" settings.
*
*       PARAMS(DZTHRESH_PARAM = 6) : Threshold where the solution is
*            considered stable enough for computing componentwise
*            measurements.
*         Default: 0.25
*
*     FACT    (input) CHARACTER*1
*     Specifies whether or not the factored form of the matrix A is
*     supplied on entry, and if not, whether the matrix A should be
*     equilibrated before it is factored.
*       = 'F':  On entry, AF and IPIV contain the factored form of A.
*               If EQUED is not 'N', the matrix A has been
*               equilibrated with scaling factors given by R and C.
*               A, AF, and IPIV are not modified.
*       = 'N':  The matrix A will be copied to AF and factored.
*       = 'E':  The matrix A will be equilibrated if necessary, then
*               copied to AF and factored.
*
*     TRANS   (input) CHARACTER*1
*     Specifies the form of the system of equations:
*       = 'N':  A * X = B     (No transpose)
*       = 'T':  A**T * X = B  (Transpose)
*       = 'C':  A**H * X = B  (Transpose)
*
*     N       (input) INTEGER
*     The number of linear equations, i.e., the order of the
*     matrix A.  N >= 0.
*
*     NRHS    (input) INTEGER
*     The number of right hand sides, i.e., the number of columns
*     of the matrices B and X.  NRHS >= 0.
*
*     A       (input/output) REAL array, dimension (LDA,N)
*     On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
*     not 'N', then A must have been equilibrated by the scaling
*     factors in R and/or C.  A is not modified if FACT = 'F' or
*     'N', or if FACT = 'E' and EQUED = 'N' on exit.
*
*     On exit, if EQUED .ne. 'N', A is scaled as follows:
*     EQUED = 'R':  A := diag(R) * A
*     EQUED = 'C':  A := A * diag(C)
*     EQUED = 'B':  A := diag(R) * A * diag(C).
*
*     LDA     (input) INTEGER
*     The leading dimension of the array A.  LDA >= max(1,N).
*
*     AF      (input or output) REAL array, dimension (LDAF,N)
*     If FACT = 'F', then AF is an input argument and on entry
*     contains the factors L and U from the factorization
*     A = P*L*U as computed by SGETRF.  If EQUED .ne. 'N', then
*     AF is the factored form of the equilibrated matrix A.
*
*     If FACT = 'N', then AF is an output argument and on exit
*     returns the factors L and U from the factorization A = P*L*U
*     of the original matrix A.
*
*     If FACT = 'E', then AF is an output argument and on exit
*     returns the factors L and U from the factorization A = P*L*U
*     of the equilibrated matrix A (see the description of A for
*     the form of the equilibrated matrix).
*
*     LDAF    (input) INTEGER
*     The leading dimension of the array AF.  LDAF >= max(1,N).
*
*     IPIV    (input or output) INTEGER array, dimension (N)
*     If FACT = 'F', then IPIV is an input argument and on entry
*     contains the pivot indices from the factorization A = P*L*U
*     as computed by SGETRF; row i of the matrix was interchanged
*     with row IPIV(i).
*
*     If FACT = 'N', then IPIV is an output argument and on exit
*     contains the pivot indices from the factorization A = P*L*U
*     of the original matrix A.
*
*     If FACT = 'E', then IPIV is an output argument and on exit
*     contains the pivot indices from the factorization A = P*L*U
*     of the equilibrated matrix A.
*
*     EQUED   (input or output) CHARACTER*1
*     Specifies the form of equilibration that was done.
*       = 'N':  No equilibration (always true if FACT = 'N').
*       = 'R':  Row equilibration, i.e., A has been premultiplied by
*               diag(R).
*       = 'C':  Column equilibration, i.e., A has been postmultiplied
*               by diag(C).
*       = 'B':  Both row and column equilibration, i.e., A has been
*               replaced by diag(R) * A * diag(C).
*     EQUED is an input argument if FACT = 'F'; otherwise, it is an
*     output argument.
*
*     R       (input or output) REAL array, dimension (N)
*     The row scale factors for A.  If EQUED = 'R' or 'B', A is
*     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
*     is not accessed.  R is an input argument if FACT = 'F';
*     otherwise, R is an output argument.  If FACT = 'F' and
*     EQUED = 'R' or 'B', each element of R must be positive.
*     If R is output, each element of R is a power of the radix.
*     If R is input, and more than working precision is used,
*     then each element of R should be a power of the radix
*     for the solution and error bound to be reliable.
*
*     C       (input or output) REAL array, dimension (N)
*     The column scale factors for A.  If EQUED = 'C' or 'B', A is
*     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
*     is not accessed.  C is an input argument if FACT = 'F';
*     otherwise, C is an output argument.  If FACT = 'F' and
*     EQUED = 'C' or 'B', each element of C must be positive.
*     If C is output, each element of C is a power of the radix.
*     If C is input, and more than working precision is used,
*     then each element of C should be a power of the radix
*     for the solution and error bound to be reliable.
*
*     B       (input/output) REAL array, dimension (LDB,NRHS)
*     On entry, the N-by-NRHS right hand side matrix B.
*     On exit,
*     if EQUED = 'N', B is not modified;
*     if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
*        diag(R)*B;
*     if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
*        overwritten by diag(C)*B.
*
*     LDB     (input) INTEGER
*     The leading dimension of the array B.  LDB >= max(1,N).
*
*     X       (output) REAL array, dimension (LDX,NRHS)
*     If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
*     to the original system of equations.  Note that A and B are
*     modified on exit if EQUED .ne. 'N', and the solution to the
*     equilibrated system is inv(diag(C))*X if TRANS = 'N' and
*     EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
*     and EQUED = 'R' or 'B'.
*
*     LDX     (input) INTEGER
*     The leading dimension of the array X.  LDX >= max(1,N).
*
*     RCOND   (output) REAL
*     Reciprocal scaled condition number.  This is an estimate of the
*     reciprocal condition number of the matrix A after equilibration
*     (if done).  If this is less than the machine precision (in
*     particular, if it is zero), the matrix is singular to working
*     precision.  This condition is indicated by a return code of
*     INFO > 0.  Note that the error may still be small even if this
*     number is very small and the matrix appears ill-conditioned.
*
*     RPVGRW  (output) REAL
*     Reciprocal pivot growth.  On exit, this contains the reciprocal
*     pivot growth factor norm(A)/norm(U). The "max absolute element"
*     norm is used.  If this is much less than 1, then the stability of
*     the LU factorization of the (equilibrated) matrix A could be poor.
*     This also means that the solution X, estimated condition numbers,
*     and error bounds could be unreliable. If factorization fails with
*     0<INFO<=N, then this contains the reciprocal pivot growth factor
*     for the leading INFO columns of A.  In SGESVX, this quantity is
*     returned in WORK(1).
*
*     BERR    (output) REAL array, dimension (NRHS)
*     Componentwise relative backward error.  Otherwise, this is the
*     componentwise relative backward error of each solution vector X(j)
*     (i.e., the smallest relative change in any element of A or B that
*     makes X(j) an exact solution).
*
*     N_ERR_BNDS (input) INTEGER
*     Number of error-bound entries to return for each right hand side
*     and each type (normwise or componentwise).  See ERR_BNDS below.
*
*     ERR_BNDS  (output) REAL array, dimension (NRHS, 2, N_ERR_BNDS)
*     The array containing information about various error
*     bounds and condition numbers.  The array is indexed by the right
*     -hand side j, the error "norm", and the type of error information 
*     as described below.  There currently are up to three pieces of 
*     information returned for each right-hand side and "norm".  If 
*     componentwise accuracy is disabled in PARAMS, the (:,2,:) entries 
*     are not accessed.  If N_ERR_BNDS .LT. 3, then at most the first
*     (:,:,N_ERR_BNDS) entries are accessed.
*
*     The first index in ERR_BNDS(j,:,:) corresponds to the jth right-
*     hand side.
*
*     The second index in ERR_BNDS(:,nrm,:) corresponds to the "norm".
*     If XTRUE(j) is the true solution for X(j) and we let 0/0 = 0, then
*     the norms are
*
*        = 1 Normwise relative:
*             max_j (abs(XTRUE(j) - X(j)))
*            ------------------------------
*                  max_j abs(X(j))
*
*        = 2 Componentwise relative
*                    abs(XTRUE(j) - X(j))
*             max_j ----------------------
*                         abs(X(j))
*
*     The final index in ERR_BNDS(:,:,err) is one of:
*        = 1 "Guaranteed" error bound: The estimated forward error,
*             almost certainly within a factor of 10 of the true error
*             so long as the next entry is greater than the threshold
*             PARAMS(RCONDTHRESH_PARAM).  If the reciprocal condition
*             number is less than that threshold, the estimate is
*             considered untrustworthy and is set to at least 1.0.
*
*        = 2  Reciprocal condition number: Estimated reciprocal condition
*             number scaled for the corresponding norm.  Compared with
*             PARAMS(RCONDTHRESH_PARAM) to determine if the error
*             estimate is "guaranteed".  (That parameter defaults to
*             sqrt(n) * slamch('Epsilon').)  These reciprocal condition
*             numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
*             appropriately scaled matrix Z.
*              Normwise: Let Z = S*A, where S scales each row by a power
*                        of the radix so all row sums of Z are
*                        approximately 1.
*              Componentwise: Let Z = S*(A*diag(x)), where x is the
*                        solution for the current right-hand side and S
*                        scales each row of A*diag(x) by a power of the
*                        radix so all row sums of Z are approximately 1.
*
*        = 3  Raw error estimate: Error estimated as detailed in LAPACK
*             Working Note 165.  This estimate will never be greater
*             than the first, "guaranteed" entry.  General systems
*             should not use or trust this error estimate.  Only use
*             this if there is some known structure to your problem that
*             is not considered in our condition number.
*     
*     See Lapack Working Note 165 for further details and extra
*     cautions.
*
*     WORK    (workspace) REAL array, dimension (5*N)
*
*     IWORK   (workspace) INTEGER array, dimension (N)
*
*     INFO    (output) INTEGER
*       = 0:  successful exit
*       < 0:  if INFO = -i, the i-th argument had an illegal value
*       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
*         has been completed, but the factor U is exactly singular, so
*         the solution and error bounds could not be computed. RCOND = 0
*         is returned.
*       = N+1: U is nonsingular, but RCOND is less than machine precision,
*         meaning that the matrix is singular to working precision.
*         Nevertheless, the solution and error bounds are computed because
*         there are a number of situations where the computed solution can
*         be more accurate than the value of RCOND would suggest.  However,
*         if INFO=N+1+J for some right-hand side J, then INFO is never set
*         to N+1.
*       = N+1+J: The iterative refinement for the J-th right-hand side did
*         not produce a small error bound.  Other right-hand sides K with 
*         K > J may have not converged as well, but only the first 
*         right-hand side to not converge is reported.
*
*     =====================================================================
      END
-------------- next part --------------
      SUBROUTINE SGERFSX(TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
     $     R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS,
     $     ERR_BNDS, NPARAMS, PARAMS, WORK, IWORK, INFO )
*
*     -- LAPACK routine (version 3.X) --
*     Copyright (c) 2005-2006, Regents of the University of California,
*     Univ. of Tennessee, NAG Ltd., Courant Institute, Argonne National
*     Lab, and Rice University
*     See COPYING file for license.
*
*     .. Scalar Arguments ..
      CHARACTER          TRANS, EQUED
      INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
     $                   N_ERR_BNDS
      REAL               RCOND, RPVGRW
*
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IWORK( * )
      REAL               A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
     $                   R( * ), C( * ), X( LDX, * ), PARAMS(*),
     $                   ERR_BNDS(NRHS, 2, N_ERR_BNDS), WORK ( * )
*
*     Purpose
*     =======
*
*     SGERFSX improves the computed solution to a system of linear
*     equations and provides error bounds and backward error estimates for
*     the solution.  In addition to normwise error bound, the code provides
*     maximum componentwise error bound if possible.  See comments for
*     ERR_BNDS for details of the error bounds.
*
*     The original system of linear equations may have been equilibrated
*     before calling this routine, as described by arguments EQUED, R and C
*     below. In this case, the solution and error bounds returned are for
*     the original unequilibrated system.
*
*     Arguments
*     =========
*
*     Some parameters are bundled in the PARAMS array.  These settings
*     affect other behavior, so we describe them first.  If the defaults
*     are acceptable, users can pass NPARAMS = 0 and the code will never
*     access the PARAMS argument.
*
*     NPARAMS (input) INTEGER
*     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
*     PARAMS array is never referenced.
*
*     PARAMS (input / output) REAL array, dimension NPARAMS
*     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
*     that entry will be filled with default value used for that
*     parameter.  Only positions up to NPARAMS are accessed; defaults
*     are used for unavailable parameters.  PARAMETER declarations for
*     the index parameters below are specified in sgerfsx_params.fh.
*
*       PARAMS(ITREF_PARAM = 1) : Whether to perform iterative 
*            refinement or not.
*         Default: 1.0
*            = 0.0 : No refinement is perfomed, and no error bounds are
*                    computed.  Error arguments are not touched.
*            = 1.0 : Use the double-precision refinement algorithm, 
*                    possibly with doubled-single computations if the 
*                    compilation environment does not support DOUBLE 
*                    PRECISION.
*              (other values are reserved for future use)
*
*       PARAMS(RCONDTHRESH_PARAM = 2) : Reciprocal condition number 
*            threshold where the error estimates are no longer 
*            considered trustworthy.
*            Change with extreme caution.
*         Default: SQRT(N) * SLAMCH('Epsilon')
*
*       PARAMS(ITHRESH_PARAM = 3) : Total number of residual computations
*            allowed for refinement.
*         Default: 10 for double-precision, 5 for single-precision
*         Aggressive: Set to 100 for our "aggresive" settings.
*
*       PARAMS(COMPONENTWISE_PARAM = 4) : Flag determining if the code
*            will attempt to find a solution with small componentwise
*            relative error in the double-precision algorithm.  Positive
*            is true, 0.0 is false.
*         Default: 1.0 (attempt componentwise convergenge)
*
*       PARAMS(RTHRESH_PARAM = 5) : Ratio of consecutive "step sizes"
*            required to continue relative normwise or componentwise
*            refinement.  If the norm of the previous step (dx or dz) is
*            SPREV and the norm of the current step is SCURR, then SCURR
*            / SPREV must be .LE. PARAMS(RTHRESH_PARAM) to continue.
*         Default: 0.5
*         Aggressive: Set to 0.9 for our "aggresive" settings.
*
*       PARAMS(DZTHRESH_PARAM = 6) : Threshold where the solution is
*            considered stable enough for computing componentwise
*            measurements.
*         Default: 0.25
*
*     TRANS   (input) CHARACTER*1
*     Specifies the form of the system of equations:
*       = 'N':  A * X = B     (No transpose)
*       = 'T':  A**T * X = B  (Transpose)
*       = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
*
*     EQUED   (input) CHARACTER*1
*     Specifies the form of equilibration that was done to A
*     before calling this routine. This is needed to compute
*     the solution and error bounds correctly.
*       = 'N':  No equilibration
*       = 'R':  Row equilibration, i.e., A has been premultiplied by diag(R).
*       = 'C':  Column equilibration, i.e., A has been postmultiplied
*               by diag(C).
*       = 'B':  Both row and column equilibration, i.e., A has been
*               replaced by diag(R) * A * diag(C).
*               The right hand side B has been changed accordingly.
*
*     N       (input) INTEGER
*     The order of the matrix A.  N >= 0.
*
*     NRHS    (input) INTEGER
*     The number of right hand sides, i.e., the number of columns
*     of the matrices B and X.  NRHS >= 0.
*
*     A       (input) REAL array, dimension (LDA,N)
*     The original N-by-N matrix A.
*
*     LDA     (input) INTEGER
*     The leading dimension of the array A.  LDA >= max(1,N).
*
*     AF      (input) REAL array, dimension (LDAF,N)
*     The factors L and U from the factorization A = P*L*U
*     as computed by SGETRF.
*
*     LDAF    (input) INTEGER
*     The leading dimension of the array AF.  LDAF >= max(1,N).
*
*     IPIV    (input) INTEGER array, dimension (N)
*     The pivot indices from SGETRF; for 1<=i<=N, row i of the
*     matrix was interchanged with row IPIV(i).
*
*     R       (input) REAL array, dimension (N)
*     The row scale factors for A.  If EQUED = 'R' or 'B', A is
*     multiplied on the left by diag(R). If EQUED = 'N' or 'C',
*     R is not accessed.  If R is input, and more than working
*     precision is used, then each element of R should be a power
*     of the radix for the solution and error bound to be reliable.
*
*     C       (input) REAL array, dimension (N)
*     The column scale factors for A.  If EQUED = 'C' or 'B', A is
*     multiplied on the right by diag(C). If EQUED = 'N' or 'R',
*     C is not accessed.  If C is input, and more than working
*     precision is used, then each element of C should be a power
*     of the radix for the solution and error bound to be reliable.
*
*     B       (input) REAL array, dimension (LDB,NRHS)
*     The right hand side matrix B.
*
*     LDB     (input) INTEGER
*     The leading dimension of the array B.  LDB >= max(1,N).
*
*     X       (input/output) REAL array, dimension (LDX,NRHS)
*     On entry, the solution matrix X, as computed by SGETRS.
*     On exit, the improved solution matrix X.
*
*     LDX     (input) INTEGER
*     The leading dimension of the array X.  LDX >= max(1,N).
*
*     RCOND   (output) REAL
*     Reciprocal scaled condition number.  This is an estimate of the
*     reciprocal condition number of the matrix A after equilibration
*     (if done).  If this is less than the machine precision (in
*     particular, if it is zero), the matrix is singular to working
*     precision.  This condition is indicated by a return code of
*     INFO > 0.  Note that the error may still be small even if this
*     number is very small and the matrix appears ill-conditioned.
*
*     RPVGRW  (output) REAL
*     Reciprocal pivot growth.  On exit, this contains the reciprocal
*     pivot growth factor norm(A)/norm(U). The "max absolute element"
*     norm is used.  If this is much less than 1, then the stability of
*     the LU factorization of the (equilibrated) matrix A could be poor.
*     This also means that the solution X, estimated condition numbers,
*     and error bounds could be unreliable. If factorization fails with
*     0<INFO<=N, then this contains the reciprocal pivot growth factor
*     for the leading INFO columns of A.  In SGESVX, this quantity is
*     returned in WORK(1).
*
*     BERR    (output) REAL array, dimension (NRHS)
*     Componentwise relative backward error.  Otherwise, this is the
*     componentwise relative backward error of each solution vector X(j)
*     (i.e., the smallest relative change in any element of A or B that
*     makes X(j) an exact solution).
*
*     N_ERR_BNDS (input) INTEGER
*     Number of error-bound entries to return for each right hand side
*     and each type (normwise or componentwise).  See ERR_BNDS below.
*
*     ERR_BNDS  (output) REAL array, dimension (NRHS, 2, N_ERR_BNDS)
*     The array containing information about various error
*     bounds and condition numbers.  The array is indexed by the right
*     -hand side j, the error "norm", and the type of error information 
*     as described below.  There currently are up to three pieces of 
*     information returned for each right-hand side and "norm".  If 
*     componentwise accuracy is disabled in PARAMS, the (:,2,:) entries 
*     are not accessed.  If N_ERR_BNDS .LT. 3, then at most the first
*     (:,:,N_ERR_BNDS) entries are accessed.
*
*     The first index in ERR_BNDS(j,:,:) corresponds to the jth right-
*     hand side.
*
*     The second index in ERR_BNDS(:,nrm,:) corresponds to the "norm".
*     If XTRUE(j) is the true solution for X(j) and we let 0/0 = 0, then
*     the norms are
*
*        = 1 Normwise relative:
*             max_j (abs(XTRUE(j) - X(j)))
*            ------------------------------
*                  max_j abs(X(j))
*
*        = 2 Componentwise relative
*                    abs(XTRUE(j) - X(j))
*             max_j ----------------------
*                         abs(X(j))
*
*     The final index in ERR_BNDS(:,:,err) is one of:
*        = 1 "Guaranteed" error bound: The estimated forward error,
*             almost certainly within a factor of 10 of the true error
*             so long as the next entry is greater than the threshold
*             PARAMS(RCONDTHRESH_PARAM).  If the reciprocal condition
*             number is less than that threshold, the estimate is
*             considered untrustworthy and is set to at least 1.0.
*
*        = 2  Reciprocal condition number: Estimated reciprocal condition
*             number scaled for the corresponding norm.  Compared with
*             PARAMS(RCONDTHRESH_PARAM) to determine if the error
*             estimate is "guaranteed".  (That parameter defaults to
*             sqrt(n) * slamch('Epsilon').)  These reciprocal condition
*             numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
*             appropriately scaled matrix Z.
*              Normwise: Let Z = S*A, where S scales each row by a power
*                        of the radix so all row sums of Z are
*                        approximately 1.
*              Componentwise: Let Z = S*(A*diag(x)), where x is the
*                        solution for the current right-hand side and S
*                        scales each row of A*diag(x) by a power of the
*                        radix so all row sums of Z are approximately 1.
*
*        = 3  Raw error estimate: Error estimated as detailed in LAPACK
*             Working Note 165.  This estimate will never be greater
*             than the first, "guaranteed" entry.  General systems
*             should not use or trust this error estimate.  Only use
*             this if there is some known structure to your problem that
*             is not considered in our condition number.
*     
*     See Lapack Working Note 165 for further details and extra
*     cautions.
*
*     WORK    (workspace) REAL array, dimension (3*N)
*
*     IWORK   (workspace) INTEGER array, dimension (N)
*
*     INFO    (output) INTEGER
*       = 0:  successful exit
*       < 0:  if INFO = -i, the i-th argument had an illegal value
*       = N+1+J: The iterative refinement for the J-th right-hand side 
*         did not produce small normwise error bound.  Other right-hand 
*         sides K with K > J may have not converged as well, but only 
*         the first right-hand side to not converge is reported.
*
*     For further details see LAPACK Working Note 165.
*
*     =====================================================================
      END

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or