Updated documentation for EPIR routines. Changed ERRS_{N,C} variable names. Other cosmetic changes.

1 files changed, 256 insertions, 8 deletions

diff --git a/SRC/cla_gbrfsx_extended.f b/SRC/cla_gbrfsx_extended.f
index 1dce776c..552bf47e 100644
--- a/SRC/cla_gbrfsx_extended.f
+++ b/SRC/cla_gbrfsx_extended.f
@@ -1,10 +1,10 @@
       SUBROUTINE CLA_GBRFSX_EXTENDED ( PREC_TYPE, TRANS_TYPE, N, KL, KU,
      $                                NRHS, AB, LDAB, AFB, LDAFB, IPIV,
      $                                COLEQU, C, B, LDB, Y, LDY,
-     $                                BERR_OUT, N_NORMS, ERRS_N, ERRS_C,
-     $                                RES, AYB, DY, Y_TAIL, RCOND,
-     $                                ITHRESH, RTHRESH, DZ_UB,
-     $                                IGNORE_CWISE, INFO )
+     $                                BERR_OUT, N_NORMS, ERR_BNDS_NORM,
+     $                                ERR_BNDS_COMP, RES, AYB, DY,
+     $                                Y_TAIL, RCOND, ITHRESH, RTHRESH,
+     $                                DZ_UB, IGNORE_CWISE, INFO )
 *
 *     -- LAPACK routine (version 3.2)                                 --
 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
@@ -27,17 +27,263 @@
       COMPLEX            AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
      $                   Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
       REAL               C( * ), AYB(*), RCOND, BERR_OUT( * ),
-     $                   ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
+     $                   ERR_BNDS_NORM( NRHS, * ),
+     $                   ERR_BNDS_COMP( NRHS, * )
 *     ..
 *
 *  Purpose
 *  =======
 *
-*  CLA_GBRFSX_EXTENDED computes ... .
+*  CLA_GBRFSX_EXTENDED improves the computed solution to a system of
+*  linear equations by performing extra-precise iterative refinement
+*  and provides error bounds and backward error estimates for the solution.
+*  This subroutine is called by CGBRFSX to perform iterative refinement.
+*  In addition to normwise error bound, the code provides maximum
+*  componentwise error bound if possible. See comments for ERR_BNDS_NORM
+*  and ERR_BNDS_COMP for details of the error bounds. Note that this
+*  subroutine is only resonsible for setting the second fields of
+*  ERR_BNDS_NORM and ERR_BNDS_COMP.
 *
 *  Arguments
 *  =========
 *
+*     PREC_TYPE      (input) INTEGER
+*     Specifies the intermediate precision to be used in refinement.
+*     The value is defined by ILAPREC(P) where P is a CHARACTER and
+*     P    = 'S':  Single
+*          = 'D':  Double
+*          = 'I':  Indigenous
+*          = 'X', 'E':  Extra
+*
+*     TRANS_TYPE     (input) INTEGER
+*     Specifies the transposition operation on A.
+*     The value is defined by ILATRANS(T) where T is a CHARACTER and
+*     T    = 'N':  No transpose
+*          = 'T':  Transpose
+*          = 'C':  Conjugate transpose
+*
+*     N              (input) INTEGER
+*     The number of linear equations, i.e., the order of the
+*     matrix A.  N >= 0.
+*
+*     KL             (input) INTEGER
+*     The number of subdiagonals within the band of A.  KL >= 0.
+*
+*     KU             (input) INTEGER
+*     The number of superdiagonals within the band of A.  KU >= 0
+*
+*     NRHS           (input) INTEGER
+*     The number of right-hand-sides, i.e., the number of columns of the
+*     matrix B.
+*
+*     AB             (input) COMPLEX array, dimension (LDA,N)
+*     On entry, the N-by-N matrix A.
+*
+*     LDAB           (input) INTEGER
+*     The leading dimension of the array A.  LDA >= max(1,N).
+*
+*     AFB            (input) COMPLEX array, dimension (LDAF,N)
+*     The factors L and U from the factorization
+*     A = P*L*U as computed by CGBTRF.
+*
+*     LDAFB          (input) INTEGER
+*     The leading dimension of the array AF.  LDAF >= max(1,N).
+*
+*     IPIV           (input) INTEGER array, dimension (N)
+*     The pivot indices from the factorization A = P*L*U
+*     as computed by CGBTRF; row i of the matrix was interchanged
+*     with row IPIV(i).
+*
+*     COLEQU         (input) LOGICAL
+*     If .TRUE. then column equilibration was done to A before calling
+*     this routine. This is needed to compute the solution and error
+*     bounds correctly.
+*
+*     C              (input) REAL array, dimension (N)
+*     The column scale factors for A. If COLEQU = .FALSE., C
+*     is not accessed. If C is input, each element of C should be a power
+*     of the radix to ensure a reliable solution and error estimates.
+*     Scaling by powers of the radix does not cause rounding errors unless
+*     the result underflows or overflows. Rounding errors during scaling
+*     lead to refining with a matrix that is not equivalent to the
+*     input matrix, producing error estimates that may not be
+*     reliable.
+*
+*     B              (input) COMPLEX array, dimension (LDB,NRHS)
+*     The right-hand-side matrix B.
+*
+*     LDB            (input) INTEGER
+*     The leading dimension of the array B.  LDB >= max(1,N).
+*
+*     Y              (input/output) COMPLEX array, dimension (LDY,NRHS)
+*     On entry, the solution matrix X, as computed by CGBTRS.
+*     On exit, the improved solution matrix Y.
+*
+*     LDY            (input) INTEGER
+*     The leading dimension of the array Y.  LDY >= max(1,N).
+*
+*     BERR_OUT       (output) REAL array, dimension (NRHS)
+*     On exit, BERR_OUT(j) contains the componentwise relative backward
+*     error for right-hand-side j from the formula
+*         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
+*     where abs(Z) is the componentwise absolute value of the matrix
+*     or vector Z. This is computed by CLA_LIN_BERR.
+*
+*     N_NORMS        (input) INTEGER
+*     Determines which error bounds to return (see ERR_BNDS_NORM
+*     and ERR_BNDS_COMP).
+*     If N_NORMS >= 1 return normwise error bounds.
+*     If N_NORMS >= 2 return componentwise error bounds.
+*
+*     ERR_BNDS_NORM  (input/output) REAL array, dimension
+*                    (NRHS, N_ERR_BNDS)
+*     For each right-hand side, this array contains information about
+*     various error bounds and condition numbers corresponding to the
+*     normwise relative error, which is defined as follows:
+*
+*     Normwise relative error in the ith solution vector:
+*             max_j (abs(XTRUE(j,i) - X(j,i)))
+*            ------------------------------
+*                  max_j abs(X(j,i))
+*
+*     The array is indexed by the type of error information as described
+*     below. There currently are up to three pieces of information
+*     returned.
+*
+*     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
+*     right-hand side.
+*
+*     The second index in ERR_BNDS_NORM(:,err) contains the following
+*     three fields:
+*     err = 1 "Trust/don't trust" boolean. Trust the answer if the
+*              reciprocal condition number is less than the threshold
+*              sqrt(n) * slamch('Epsilon').
+*
+*     err = 2 "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
+*              sqrt(n) * slamch('Epsilon'). This error bound should only
+*              be trusted if the previous boolean is true.
+*
+*     err = 3  Reciprocal condition number: Estimated normwise
+*              reciprocal condition number.  Compared with the threshold
+*              sqrt(n) * slamch('Epsilon') to determine if the error
+*              estimate is "guaranteed". These reciprocal condition
+*              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
+*              appropriately scaled matrix Z.
+*              Let Z = S*A, where S scales each row by a power of the
+*              radix so all absolute row sums of Z are approximately 1.
+*
+*     This subroutine is only responsible for setting the second field
+*     above.
+*     See Lapack Working Note 165 for further details and extra
+*     cautions.
+*
+*     ERR_BNDS_COMP  (input/output) REAL array, dimension
+*                    (NRHS, N_ERR_BNDS)
+*     For each right-hand side, this array contains information about
+*     various error bounds and condition numbers corresponding to the
+*     componentwise relative error, which is defined as follows:
+*
+*     Componentwise relative error in the ith solution vector:
+*                    abs(XTRUE(j,i) - X(j,i))
+*             max_j ----------------------
+*                         abs(X(j,i))
+*
+*     The array is indexed by the right-hand side i (on which the
+*     componentwise relative error depends), and the type of error
+*     information as described below. There currently are up to three
+*     pieces of information returned for each right-hand side. If
+*     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
+*     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
+*     the first (:,N_ERR_BNDS) entries are returned.
+*
+*     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
+*     right-hand side.
+*
+*     The second index in ERR_BNDS_COMP(:,err) contains the following
+*     three fields:
+*     err = 1 "Trust/don't trust" boolean. Trust the answer if the
+*              reciprocal condition number is less than the threshold
+*              sqrt(n) * slamch('Epsilon').
+*
+*     err = 2 "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
+*              sqrt(n) * slamch('Epsilon'). This error bound should only
+*              be trusted if the previous boolean is true.
+*
+*     err = 3  Reciprocal condition number: Estimated componentwise
+*              reciprocal condition number.  Compared with the threshold
+*              sqrt(n) * slamch('Epsilon') to determine if the error
+*              estimate is "guaranteed". These reciprocal condition
+*              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
+*              appropriately scaled matrix Z.
+*              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 absolute row
+*              sums of Z are approximately 1.
+*
+*     This subroutine is only responsible for setting the second field
+*     above.
+*     See Lapack Working Note 165 for further details and extra
+*     cautions.
+*
+*     RES            (input) COMPLEX array, dimension (N)
+*     Workspace to hold the intermediate residual.
+*
+*     AYB            (input) REAL array, dimension (N)
+*     Workspace.
+*
+*     DY             (input) COMPLEX array, dimension (N)
+*     Workspace to hold the intermediate solution.
+*
+*     Y_TAIL         (input) COMPLEX array, dimension (N)
+*     Workspace to hold the trailing bits of the intermediate solution.
+*
+*     RCOND          (input) REAL
+*     Reciprocal scaled condition number.  This is an estimate of the
+*     reciprocal Skeel 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.  Note that the error may still be small even
+*     if this number is very small and the matrix appears ill-
+*     conditioned.
+*
+*     ITHRESH        (input) INTEGER
+*     The maximum number of residual computations allowed for
+*     refinement. The default is 10. For 'aggressive' set to 100 to
+*     permit convergence using approximate factorizations or
+*     factorizations other than LU. If the factorization uses a
+*     technique other than Gaussian elimination, the guarantees in
+*     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
+*
+*     RTHRESH        (input) REAL
+*     Determines when to stop refinement if the error estimate stops
+*     decreasing. Refinement will stop when the next solution no longer
+*     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
+*     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
+*     default value is 0.5. For 'aggressive' set to 0.9 to permit
+*     convergence on extremely ill-conditioned matrices. See LAWN 165
+*     for more details.
+*
+*     DZ_UB          (input) REAL
+*     Determines when to start considering componentwise convergence.
+*     Componentwise convergence is only considered after each component
+*     of the solution Y is stable, which we definte as the relative
+*     change in each component being less than DZ_UB. The default value
+*     is 0.25, requiring the first bit to be stable. See LAWN 165 for
+*     more details.
+*
+*     IGNORE_CWISE   (input) LOGICAL
+*     If .TRUE. then ignore componentwise convergence. Default value
+*     is .FALSE..
+*
+*     INFO           (output) INTEGER
+*       = 0:  Successful exit.
+*       < 0:  if INFO = -i, the ith argument to CGBTRS had an illegal
+*             value
+*
 *  =====================================================================
 *
 *     .. Local Scalars ..
@@ -282,10 +528,12 @@
 *     Compute error bounds.
 *
          IF ( N_NORMS .GE. 1 ) THEN
-            ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX)
+            ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) =
+     $           FINAL_DX_X / (1 - DXRATMAX)
          END IF
          IF ( N_NORMS .GE. 2 ) THEN
-            ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX)
+            ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) =
+     $           FINAL_DZ_Z / (1 - DZRATMAX)
          END IF
 *
 *     Compute componentwise relative backward error from formula