la_herfsx_extended(3)

📦
version 3.12.1 · 9 versions
, in arch rolling

okay so bestie, this is la_herfsx_extended(3), and the manual is Library Functions Manual — welcome to the group chat. 💅

💖 NAME 💖

la_herfsx_extended - la_herfsx_extended: step in herfsx. That’s the name, no cap.

✨ SYNOPSIS ✨

📚 Functions

cla_herfsx_extended

subroutine cla_herfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

CLA_HERFSX_EXTENDED improves the computed solution to a system of linear equations for Hermitian indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

cla_syrfsx_extended

subroutine cla_syrfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

CLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

dla_syrfsx_extended

subroutine dla_syrfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

DLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

sla_syrfsx_extended

subroutine sla_syrfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

SLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

zla_herfsx_extended

subroutine zla_herfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

ZLA_HERFSX_EXTENDED improves the computed solution to a system of linear equations for Hermitian indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

zla_syrfsx_extended

subroutine zla_syrfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)

ZLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

📝 Detailed Description 📝

lowkey the detailed description is just the function docs below — no extra lore, no side quests. 🫶

🛠️ Function Documentation 🛠️

cla_herfsx_extended

subroutine cla_herfsx_extended (integer prec_type, character uplo, integer
    n, integer nrhs, complex, dimension( lda, * ) a, integer lda, complex,
    dimension( ldaf, * ) af, integer ldaf, integer, dimension( * ) ipiv,
    logical colequ, real, dimension( * ) c, complex, dimension( ldb, * ) b,
    integer ldb, complex, dimension( ldy, * ) y, integer ldy, real,
    dimension( * ) berr_out, integer n_norms, real, dimension( nrhs, * )
    err_bnds_norm, real, dimension( nrhs, * ) err_bnds_comp, complex,
    dimension( * ) res, real, dimension( * ) ayb, complex, dimension( * )
    dy, complex, dimension( * ) y_tail, real rcond, integer ithresh, real
    rthresh, real dz_ub, logical ignore_cwise, integer info)

CLA_HERFSX_EXTENDED improves the computed solution to a system of linear equations for Hermitian indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

🎯 Purpose 🎯

CLA_HERFSX_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 CHERFSX 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 responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP.

🧾 Parameters 🧾

PREC_TYPE

PREC_TYPE is INTEGER. It 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
  • P = 'D': Double
  • P = 'I': Indigenous
  • P = 'X' or 'E': Extra

UPLO

UPLO is CHARACTER*1.

  • = 'U': Upper triangle of A is stored;
  • = 'L': Lower triangle of A is stored.

N

N is INTEGER. It’s the number of linear equations, i.e., the order of the matrix A. N >= 0.

NRHS

NRHS is INTEGER. It’s the number of right-hand-sides, i.e., the number of columns of the matrix B.

A

A is COMPLEX array, dimension (LDA,N). On entry, the N-by-N matrix A.

LDA

LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,N).

AF

AF is COMPLEX array, dimension (LDAF,N). The block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by CHETRF.

LDAF

LDAF is INTEGER. The leading dimension of the array AF. LDAF >= max(1,N).

IPIV

IPIV is INTEGER array, dimension (N). Details of the interchanges and the block structure of D as determined by CHETRF.

COLEQU

COLEQU is 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

C is 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

B is COMPLEX array, dimension (LDB,NRHS). The right-hand-side matrix B.

LDB

LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).

Y

Y is COMPLEX array, dimension (LDY,NRHS). On entry, the solution matrix X, as computed by CHETRS. On exit, the improved solution matrix Y.

LDY

LDY is INTEGER. The leading dimension of the array Y. LDY >= max(1,N).

BERR_OUT

BERR_OUT is 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

N_NORMS is INTEGER. It 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

🎀 ERR_BNDS_NORM 🎀

ERR_BNDS_NORM is a REAL array, dimension (NRHS, N_ERR_BNDS). For each right-hand side, this array is carrying info about various error bounds and condition numbers corresponding to the normwise relative error, which is defined like this:

        max_j (abs(XTRUE(j,i) - X(j,i)))
       ------------------------------
             max_j abs(X(j,i))

The array is indexed by the type of error info 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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 💖

ERR_BNDS_COMP is a REAL array, dimension (NRHS, N_ERR_BNDS). For each right-hand side, this array is holding info about various error bounds and condition numbers corresponding to the componentwise relative error, defined like this:

               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 info 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 < 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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 🌟

RES is COMPLEX array, dimension (N). Workspace to hold the intermediate residual.

🌈 AYB 🌈

AYB is REAL array, dimension (N). Workspace.

🦋 DY 🦋

DY is COMPLEX array, dimension (N). Workspace to hold the intermediate solution.

🎀 Y_TAIL 🎀

Y_TAIL is COMPLEX array, dimension (N). Workspace to hold the trailing bits of the intermediate solution.

💅 RCOND 💅

RCOND is 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 🔥

ITHRESH is 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 🌸

RTHRESH is 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 🧚

DZ_UB is REAL. Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solution Y is stable, which we define 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 🙈

IGNORE_CWISE is LOGICAL. If .TRUE. then ignore componentwise convergence. Default value is .FALSE..

💀 INFO 💀

INFO is INTEGER

  • = 0: Successful exit.
  • < 0: if INFO = -i, the ith argument to CLA_HERFSX_EXTENDED had an illegal value

Author

  • Univ. of Tennessee
  • Univ. of California Berkeley
  • Univ. of Colorado Denver
  • NAG Ltd.

📍 Definition at line 388 of file cla_herfsx_extended.f

🦋 cla_syrfsx_extended 🦋

subroutine cla_syrfsx_extended (integer prec_type, character uplo, integer
    n, integer nrhs, complex, dimension( lda, * ) a, integer lda, complex,
    dimension( ldaf, * ) af, integer ldaf, integer, dimension( * ) ipiv,
    logical colequ, real, dimension( * ) c, complex, dimension( ldb, * ) b,
    integer ldb, complex, dimension( ldy, * ) y, integer ldy, real,
    dimension( * ) berr_out, integer n_norms, real, dimension( nrhs, * )
    err_bnds_norm, real, dimension( nrhs, * ) err_bnds_comp, complex,
    dimension( * ) res, real, dimension( * ) ayb, complex, dimension( * )
    dy, complex, dimension( * ) y_tail, real rcond, integer ithresh, real
    rthresh, real dz_ub, logical ignore_cwise, integer info)

CLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

🎯 Purpose 🎯

CLA_SYRFSX_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 CSYRFSX 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 responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP.

🌸 Parameters 🌸

  • PREC_TYPE — 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' or 'E': Extra
  • UPLO — CHARACTER*1
    • = 'U': Upper triangle of A is stored;
    • = 'L': Lower triangle of A is stored.
  • N — INTEGER. The number of linear equations, i.e., the order of the matrix A. N >= 0.
  • NRHS — INTEGER. The number of right-hand-sides, i.e., the number of columns of the matrix B.
  • A — COMPLEX array, dimension (LDA,N). On entry, the N-by-N matrix A.
  • LDA
  • LDA — LDA is INTEGER. It's the leading dimension of the array A, and the vibe check it must pass is LDA >= max(1,N). Non-negotiable bestie.
  • AF — AF is COMPLEX array, dimension (LDAF,N). This one's carrying the block diagonal matrix D plus the multipliers used to obtain the factor U or L, exactly as computed by CSYTRF. She's the pre-processed version, basically.
  • LDAF — LDAF is INTEGER. Leading dimension of the array AF, and yes, same rule: LDAF >= max(1,N).
  • IPIV — IPIV is INTEGER array, dimension (N). It holds the details of the interchanges and the block structure of D, as determined by CSYTRF. Receipts for the pivoting, iykyk.
  • COLEQU — COLEQU is LOGICAL. If .TRUE. then column equilibration was done to A before calling this routine. This is genuinely needed to compute the solution and error bounds correctly, so don't be out here lying about it 💀
  • C — C is REAL array, dimension (N). The column scale factors for A. If COLEQU = .FALSE., C is not accessed at all (she's not even invited). 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 — and that is a whole mess you do not want 😭
  • B — B is COMPLEX array, dimension (LDB,NRHS). It's the right-hand-side matrix B. The moment. The whole reason we're here.
  • LDB — LDB is INTEGER. Leading dimension of the array B, and it's giving LDB >= max(1,N), as always.
  • Y — Y is COMPLEX array, dimension (LDY,NRHS). On entry, this is the solution matrix X, as computed by CSYTRS. On exit, it's the improved solution matrix Y. Glow-up arc, truly ✨
  • LDY — LDY is INTEGER. Leading dimension of the array Y. You already know the drill: LDY >= max(1,N).
  • BERR_OUT — BERR_OUT is REAL array, dimension (NRHS). On exit, BERR_OUT(j) contains the componentwise relative backward error for right-hand-side j, from this 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, so shoutout to her 🫶

  • N_NORMS — N_NORMS is INTEGER, and it 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 — ERR_BNDS_NORM is 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 boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 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 . 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.
Okay so that one sentence in err = 3 literally just stops at "is ." in the source — I'm not gonna invent the missing word for you 🙈 the manpage said what it said.

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 — ERR_BNDS_COMP is 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 < 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 boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 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 . 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.
Same situation as above, bestie — "is ." is verbatim what's in the source. Not me filling in blanks 💅

This subroutine is only responsible for setting the second field above. See Lapack Working Note 165 for further details and extra cautions.

  • RES — RES is COMPLEX array, dimension (N). Workspace to hold the intermediate residual. Just vibing in the back, doing her job.
  • AYB — AYB is REAL array, dimension (N). Workspace. That's it. That's the description 😌
  • DY — DY is COMPLEX array, dimension (N). Workspace to hold the intermediate solution.
  • Y_TAIL — Y_TAIL is COMPLEX array, dimension (N). Workspace to hold the trailing bits of the intermediate solution. The little leftovers, the crumbs, the detail 🌸
  • RCOND — RCOND is 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 — so like, don't panic-panic, but do read the whole line 👑
  • ITHRESH — ITHRESH is 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. Big warning, no cap ⚠️
  • RTHRESH — RTHRESH is 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 — DZ_UB is REAL. Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solution Y is stable, which we define 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 — IGNORE_CWISE is LOGICAL. If .TRUE. then ignore componentwise convergence. Default value is .FALSE.. Sometimes ignoring your problems is the option, and here it's literally a flag 🦋
  • INFO — INFO is INTEGER.
    • = 0: Successful exit. Slay.
    • < 0: if INFO = -i, the ith argument to CLA_SYRFSX_EXTENDED had an illegal value.

🌈 Author

Univ. of Tennessee

📍 Definition at line 388 of file cla_syrfsx_extended.f.

🦋 dla_syrfsx_extended 🦋

subroutine dla_syrfsx_extended (integer prec_type, character uplo, integer
    n, integer nrhs, double precision, dimension( lda, * ) a, integer lda,
    double precision, dimension( ldaf, * ) af, integer ldaf, integer,
    dimension( * ) ipiv, logical colequ, double precision, dimension( * ) c,
    double precision, dimension( ldb, * ) b, integer ldb, double precision,
    dimension( ldy, * ) y, integer ldy, double precision, dimension( * )
    berr_out, integer n_norms, double precision, dimension( nrhs, * )
    err_bnds_norm, double precision, dimension( nrhs, * ) err_bnds_comp,
    double precision, dimension( * ) res, double precision, dimension( * )
    ayb, double precision, dimension( * ) dy, double precision, dimension( *
    ) y_tail, double precision rcond, integer ithresh, double precision
    rthresh, double precision dz_ub, logical ignore_cwise, integer info)

DLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. okay so it's the glow-up routine for your solution, basically 💅

🌟 Purpose 🌟

!> !> !> DLA_SYRFSX_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 DSYRFSX 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 responsible for setting the second fields of !> ERR_BNDS_NORM and ERR_BNDS_COMP. !>

Translation bestie: it's the muscle behind DSYRFSX. It refines, it bounds the error, and it only fills in field two of those two ERR_BNDS arrays — the second fields, not the whole thing. Don't @ me, that's literally what it says. 🫶

📋 Parameters 📋

PREC_TYPE — INTEGER

PREC_TYPE is 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' or 'E': Extra

UPLO — CHARACTER1*

UPLO is CHARACTER*1

  • = 'U': Upper triangle of A is stored;
  • = 'L': Lower triangle of A is stored.

N — INTEGER

N is INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.

NRHS — INTEGER

NRHS is INTEGER The number of right-hand-sides, i.e., the number of columns of the matrix B.

A — DOUBLE PRECISION array, dimension (LDA,N)

A is DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N matrix A.

LDA — INTEGER

LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).

AF — DOUBLE PRECISION array, dimension (LDAF,N)

AF is DOUBLE PRECISION array, dimension (LDAF,N) The block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by DSYTRF.

LDAF — INTEGER

LDAF is INTEGER The leading dimension of the array AF. LDAF >= max(1,N).

IPIV — INTEGER array, dimension (N)

IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D as determined by DSYTRF.

COLEQU — LOGICAL

COLEQU is 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. (so like, if you equilibrated and didn't tell it, the vibes are off and the bounds are wrong — set it honestly 💀)

C — DOUBLE PRECISION array, dimension (N)

C is DOUBLE PRECISION 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.

🚨 real talk: powers of the radix ONLY, bestie. Rounding errors during scaling mean you're refining against a matrix that isn't even your matrix anymore, and then those error estimates are giving delusion, not accuracy.

B — DOUBLE PRECISION array, dimension (LDB,NRHS)

B is DOUBLE PRECISION array, dimension (LDB,NRHS) The right-hand-side matrix B.

LDB — INTEGER

LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).

Y — DOUBLE PRECISION array, dimension (LDY,NRHS)

Y is DOUBLE PRECISION array, dimension (LDY,NRHS) On entry, the solution matrix X, as computed by DSYTRS. On exit, the improved solution matrix Y. ✨ the glow-up is real: goes in as X, comes out as Y ✨

LDY — INTEGER

LDY is INTEGER The leading dimension of the array Y. LDY >= max(1,N).

BERR_OUT — DOUBLE PRECISION array, dimension (NRHS)

BERR_OUT is DOUBLE PRECISION 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 DLA_LIN_BERR. iykyk 🧮

N_NORMS — INTEGER

N_NORMS is 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.

So like a little menu: 1 gets you normwise, 2 gets you both. Choose your fighter 💅

ERR_BNDS_NORM — DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)

ERR_BNDS_NORM is DOUBLE PRECISION 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:

        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 *i*th right-hand side.

The second index in ERR_BNDS_NORM(:,err) contains the following three fields:

  • err = 1 boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 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 . 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.
Heads up: the source doc genuinely trails off at "the error estimate is ." — not my typo, bestie, that's just how it's written 😭 I'm not about to invent words for it. Everything else in that field is exactly as documented.

This subroutine is only responsible for setting the second field above. See Lapack Working Note 165 for further details and extra cautions.

LWN 165 is the assigned reading, no skipping 👑

ERR_BNDS_COMP — and this is exactly where this chunk of the docs cuts off 🦋 the componentwise twin of ERR_BNDS_NORM — its description picks up right where this part ends, so hold that thought 💜

  • ERR_BNDS_COMP — DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) 💖 for every right-hand side, this one carries all the tea about error bounds and condition numbers tied to the componentwise relative error, which is defined like this:
                   abs(XTRUE(j,i) - X(j,i))
            max_j ----------------------
                         abs(X(j,i))

The array is indexed by the right-hand side i (because that's what the componentwise relative error depends on, obvi), and then by the type of error info, as described below. Right now there are up to three pieces of info returned per right-hand side. If componentwise accuracy was not requested (PARAMS(3) = 0.0), then ERR_BNDS_COMP is not accessed at all — it's giving "you didn't ask, so I'm not telling" energy. And if N_ERR_BNDS < 3, then at most the first (:,N_ERR_BNDS) entries get returned. Lowkey thrifty but valid.

The first index in ERR_BNDS_COMP(i,:) corresponds to the ith right-hand side 🫶

The second index in ERR_BNDS_COMP(:,err) contains these three fields:

  • err = 1 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon'). Basically a vibe check on whether you can believe the number ✨
  • err = 2 — 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'). And no cap, this error bound should only be trusted if the previous boolean is true. Package deal, bestie 💅
  • err = 3 — reciprocal condition number: the estimated componentwise reciprocal condition number. Compared with the threshold sqrt(n) * slamch('Epsilon') to determine if the error estimate is . 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. Math girlie behaviour honestly 🌸

This subroutine is only responsible for setting the second field above. That's it, that's the whole job. See Lapack Working Note 165 for further details and extra cautions 👑

  • RES — DOUBLE PRECISION array, dimension (N). Workspace to hold the intermediate residual. Scrappy little scratchpad 🎀
  • AYB — DOUBLE PRECISION array, dimension (N). Workspace. And here's the slay: this can be the same workspace you passed for Y_TAIL — sharing is caring, we love a sustainable queen 🦋
  • DY — DOUBLE PRECISION array, dimension (N). Workspace to hold the intermediate solution. Mid-refinement era, not the final form 💜
  • Y_TAIL — DOUBLE PRECISION array, dimension (N). Workspace to hold the trailing bits of the intermediate solution. The leftovers, the crumbs, and honestly they matter 😭
  • RCOND — DOUBLE PRECISION. 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 — that's the dramatic part, don't skip it. Note that the error may still be small even if this number is very small and the matrix appears ill-conditioned 🫶 (yes, it can look rough and still eat).
  • ITHRESH — INTEGER. The maximum number of residual computations allowed for refinement. The default is 10. For 'aggressive' set it to 100 to permit convergence using approximate factorizations or factorizations other than LU. Caveat time: 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 ⚠️ like, proceed but know what you're doing.
  • RTHRESH — DOUBLE PRECISION. 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 it to 0.9 to permit convergence on extremely ill-conditioned matrices. See LAWN 165 for more details 🔥
  • DZ_UB — DOUBLE PRECISION. Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solution Y is stable, which we define 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 — LOGICAL. If .TRUE. then ignore componentwise convergence. Default value is .FALSE.. Setting it true is the "I don't want to know" option 🙈
  • INFO — INTEGER:
    • = 0: Successful exit. We love a clean finish 💖
    • < 0: if INFO = -i, the ith argument to DLA_SYRFSX_EXTENDED had an illegal value. Not iconic, fix your args bestie 💀

🌸 Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 389 of file dla_syrfsx_extended.f.

🌈 subroutine sla_syrfsx_extended 🌈

subroutine sla_syrfsx_extended (integer prec_type, character uplo, integer
n, integer nrhs, real, dimension( lda, * ) a, integer lda, real,
dimension( ldaf, * ) af, integer ldaf, integer, dimension( * ) ipiv,
logical colequ, real, dimension( * ) c, real, dimension( ldb, * ) b,
integer ldb, real, dimension( ldy, * ) y, integer ldy, real, dimension(
* ) berr_out, integer n_norms, real, dimension( nrhs, * ) err_bnds_norm,
real, dimension( nrhs, * ) err_bnds_comp, real, dimension( * ) res,
real, dimension( * ) ayb, real, dimension( * ) dy, real, dimension( * )
y_tail, real rcond, integer ithresh, real rthresh, real dz_ub, logical
ignore_cwise, integer info)

SLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement, and provides error bounds and backward error estimates for the solution. She does it all 💅

✨ Purpose ✨

SLA_SYRFSX_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 SSYRFSX to perform iterative refinement. In addition to the normwise error bound, the code provides a maximum componentwise error bound if possible. See the comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the error bounds. Note that this subroutine is only responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP — just the second ones, don't ask her for more 🎀

💜 Parameters 💜

  • PREC_TYPE — 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' or 'E': Extra

Pick your fighter ✨

  • UPLO — CHARACTER*1:
    • = 'U': Upper triangle of A is stored;
    • = 'L': Lower triangle of A is stored.

Half the matrix is all you get, and that's fine 🌸

  • N — INTEGER. The number of linear equations, i.e., the order of the matrix A. N >= 0.
  • NRHS — INTEGER. The number of right-hand-sides, i.e., the number of columns of the matrix B. Multiple problems, one queen 👑
  • A — REAL array, dimension (LDA,N). On entry, the N-by-N matrix A. The main character herself 💖
  • LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N).
  • AF — REAL array, dimension (LDAF,N). The block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by SSYTRF.
  • LDAF — INTEGER. The leading dimension of the array AF. LDAF >= max(1,N).
  • IPIV — INTEGER array, dimension (N). Details of the interchanges and the block structure of D as determined by SSYTRF. The receipts of what got swapped where 🧚
  • COLEQU — 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 — skip it and everything's just wrong, no cap ⚠️
  • C — REAL array, dimension (N). The column scale factors for A. If COLEQU = .FALSE., C is not accessed (she's off duty). 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. So bestie: powers of the radix, or risk the consequences 😭
  • B — REAL array, dimension (LDB,NRHS). The right-hand-side matrix B.
  • LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
  • Y — REAL array, dimension (LDY,NRHS). On entry, the solution matrix X, as computed by SSYTRS. On exit, the improved solution matrix Y. Glow-up arc, literally 🌈
  • LDY — INTEGER. The leading dimension of the array Y. LDY >= max(1,N).
  • BERR_OUT —
  • BERR_OUT — 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 SLA_LIN_BERR. 💅

  • N_NORMS — 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 — 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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 — 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 < 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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 — REAL array, dimension (N). Workspace to hold the intermediate residual.
  • AYB — REAL array, dimension (N). Workspace. This can be the same workspace passed for Y_TAIL.
  • DY — REAL array, dimension (N). Workspace to hold the intermediate solution.
  • Y_TAIL — REAL array, dimension (N). Workspace to hold the trailing bits of the intermediate solution.
  • RCOND — 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 — 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 — 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 — REAL. Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solution Y is stable, which we define 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 — LOGICAL. If .TRUE. then ignore componentwise convergence. Default value is .FALSE..
  • INFO — INTEGER.
    • = 0: Successful exit.
    • < 0: if INFO = -i, the ith argument to SLA_SYRFSX_EXTENDED had an illegal value

👑 Author 👑

  • Univ. of Tennessee
  • Univ. of California Berkeley
  • Univ. of Colorado Denver
  • NAG Ltd.

Definition at line 389 of file sla_syrfsx_extended.f.

okay bestie, next up is the complex Hermitian version, and she is iconic: zla_herfsx_extended. ✨

subroutine zla_herfsx_extended (integer prec_type, character uplo, integer
       n, integer nrhs, complex*16, dimension( lda, * ) a, integer lda,
       complex*16, dimension( ldaf, * ) af, integer ldaf, integer, dimension( *
       ) ipiv, logical colequ, double precision, dimension( * ) c, complex*16,
       dimension( ldb, * ) b, integer ldb, complex*16, dimension( ldy, * ) y,
       integer ldy, double precision, dimension( * ) berr_out, integer n_norms,
       double precision, dimension( nrhs, * ) err_bnds_norm, double precision,
       dimension( nrhs, * ) err_bnds_comp, complex*16, dimension( * ) res,
       double precision, dimension( * ) ayb, complex*16, dimension( * ) dy,
       complex*16, dimension( * ) y_tail, double precision rcond, integer
       ithresh, double precision rthresh, double precision dz_ub, logical
       ignore_cwise, integer info)

ZLA_HERFSX_EXTENDED improves the computed solution to a system of linear equations for Hermitian indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. 💜

Purpose

ZLA_HERFSX_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 ZHERFSX 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 responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP. 🎀

Parameters

  • PREC_TYPE —
  • PREC_TYPE (INTEGER) — this one tells the routine which intermediate precision to use while it's doing the refinement 💅 it's defined by ILAPREC(P) where P is a CHARACTER, and the values are:
    • = 'S' → Single
    • = 'D' → Double
    • = 'I' → Indigenous
    • = 'X' or 'E' → Extra
  • UPLO (CHARACTER*1) — which triangle of A we actually stored, bestie:
    • = 'U' → the Upper triangle of A is stored
    • = 'L' → the Lower triangle of A is stored ✨
  • N (INTEGER) — how many linear equations we're solving, aka the order of the matrix A. Gotta be N >= 0.
  • NRHS (INTEGER) — the number of right-hand-sides, meaning the number of columns of the matrix B. So many b's, so little time 🫶
  • A (COMPLEX*16 array, dimension (LDA,N)) — on entry, this is the N-by-N matrix A. The main character, honestly 👑
  • LDA (INTEGER) — the leading dimension of the array A. Rule: LDA >= max(1,N).
  • AF (COMPLEX*16 array, dimension (LDAF,N)) — the block diagonal matrix D plus the multipliers used to obtain the factor U or L, all as computed by ZHETRF. It's the factorization, already served 🍓
  • LDAF (INTEGER) — the leading dimension of the array AF. Rule: LDAF >= max(1,N).
  • IPIV (INTEGER array, dimension (N)) — the deets on the interchanges and the block structure of D, as determined by ZHETRF. The receipts of who got swapped with whom 💅
  • COLEQU (LOGICAL) — if .TRUE., then column equilibration was already done to A before this routine got called. This matters because it's needed to compute the solution and the error bounds correctly — so don't lie to it, k? 🙈
  • C (DOUBLE PRECISION array, dimension (N)) — the column scale factors for A. If COLEQU = .FALSE., then C is not accessed at all, it's just vibing. If C is input, each element of C should be a power of the radix so you get a reliable solution and reliable error estimates. Scaling by powers of the radix doesn't cause rounding errors unless the result underflows or overflows — and rounding errors during scaling lead to refining with a matrix that is not equivalent to the input matrix, which produces error estimates that may not be reliable. So like… be careful, fr fr 😭
  • B (COMPLEX*16 array, dimension (LDB,NRHS)) — the right-hand-side matrix B. The whole reason we're here 💖
  • LDB (INTEGER) — the leading dimension of the array B. Rule: LDB >= max(1,N).
  • Y (COMPLEX*16 array, dimension (LDY,NRHS)) — on entry, this holds the solution matrix X as computed by ZHETRS; on exit, it's the improved solution matrix Y. Glow-up arc ✨🔥
  • LDY (INTEGER) — the leading dimension of the array Y. Rule: LDY >= max(1,N).
  • BERR_OUT (DOUBLE PRECISION array, dimension (NRHS)) — on exit, BERR_OUT(j) contains the componentwise relative backward error for right-hand-side j, straight from this 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 ZLA_LIN_BERR. Math slay 🧮👑

  • N_NORMS (INTEGER) — decides which error bounds come back to you (peep 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 (DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)) — for each right-hand side, this array is packed with info about various error bounds and condition numbers for the normwise relative error, defined like this:
    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, 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 these three fields:

  • err = 1 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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.

Heads up: 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 (DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)) — for each right-hand side, this array holds info about various error bounds and condition numbers for the componentwise relative error, defined as:
    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 (since the componentwise relative error depends on it) and by the type of error information 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 < 3, then at most the first (:,N_ERR_BNDS) entries are returned — so you might get a partial drop 🙈

The first index in ERR_BNDS_COMP(i,:) corresponds to the ith right-hand side.

The second index in ERR_BNDS_COMP(:,err) contains these three fields:

  • err = 1 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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.

Again: this subroutine is only responsible for setting the second field above, and you should see Lapack Working Note 165 for further details and extra cautions 🎀

  • RES (COMPLEX*16 array, dimension (N)) — workspace that holds the intermediate residual. Just a lil scratch pad 🧚
  • AYB (DOUBLE PRECISION array, dimension (N)) — workspace. That's it, that's the whole description 💀
  • DY (COMPLEX*16 array, dimension (N)) — workspace that holds the intermediate solution.
  • Y_TAIL (COMPLEX*16 array, dimension (N)) — workspace that holds the trailing bits of the intermediate solution. The crumbs, the leftover sparkle ✨
  • RCOND (DOUBLE PRECISION) — the reciprocal scaled condition number. It's an estimate of the reciprocal Skeel condition number of the matrix A after equilibration (if that got done). If this is less than the machine precision (in particular, if it is zero), then the matrix is singular to working precision. Big note though: the error may still be small even if this number is very small and the matrix appears ill-conditioned. Don't panic prematurely, bestie 😭🫶
  • ITHRESH (INTEGER) — the max number of residual computations allowed for refinement. Default is 10. For 'aggressive' set it to 100 to permit convergence using approximate factorizations, or factorizations other than LU. ⚠️ But if the factorization uses a technique other than Gaussian elimination, then the guarantees in ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy. Slippery slope, no cap.
  • RTHRESH
  • RTHRESH — DOUBLE PRECISION</p>. This is the ✨when-do-I-give-up✨ knob: it decides when to stop refining if the error estimate stops decreasing. Refinement stops the moment 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. Default is 0.5. For the 'aggressive' vibes, set it to 0.9 to permit convergence on extremely ill-conditioned matrices. See LAWN 165 for more details. 💅
  • DZ_UB — DOUBLE PRECISION</p>. Determines when to start considering componentwise convergence. And no, bestie, componentwise convergence is not invited to the party early — it's only considered after each component of the solution Y is stable, which we define as the relative change in each component being less than DZ_UB. Default is 0.25, requiring the first bit to be stable. See LAWN 165 for more details. 🎀
  • IGNORE_CWISE — LOGICAL</p>. If .TRUE. then ignore componentwise convergence. Default value is .FALSE.. (so by default we do care, bestie 😌)
  • INFO — INTEGER</p>.
    • = 0: Successful exit. slay 🎉
    • < 0: if INFO = -i, the ith argument to ZLA_HERFSX_EXTENDED had an illegal value 💀

🌸 Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 388 of file zla_herfsx_extended.f.

subroutine zla_syrfsx_extended

subroutine zla_syrfsx_extended (integer prec_type, character uplo, integer
    n, integer nrhs, complex*16, dimension( lda, * ) a, integer lda,
    complex*16, dimension( ldaf, * ) af, integer ldaf, integer, dimension( *
    ) ipiv, logical colequ, double precision, dimension( * ) c, complex*16,
    dimension( ldb, * ) b, integer ldb, complex*16, dimension( ldy, * ) y,
    integer ldy, double precision, dimension( * ) berr_out, integer n_norms,
    double precision, dimension( nrhs, * ) err_bnds_norm, double precision,
    dimension( nrhs, * ) err_bnds_comp, complex*16, dimension( * ) res,
    double precision, dimension( * ) ayb, complex*16, dimension( * ) dy,
    complex*16, dimension( * ) y_tail, double precision rcond, integer
    ithresh, double precision rthresh, double precision dz_ub, logical
    ignore_cwise, integer info)

ZLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. 🦋

Purpose

okay so, for real: ZLA_SYRFSX_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 ZSYRFSX to perform iterative refinement. In addition to the normwise error bound, the code provides the maximum componentwise error bound if possible. See the comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the error bounds. Note that this subroutine is only responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP. 👑 (yes, only the second field, she's not doing the whole thing, she's got boundaries)

Parameters

  • PREC_TYPE — INTEGER</p>. 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' or 'E': Extra
  • UPLO — CHARACTER*1</p>.
    • = 'U': Upper triangle of A is stored;
    • = 'L': Lower triangle of A is stored.
  • N — INTEGER</p>. The number of linear equations, i.e., the order of the matrix A. N >= 0.
  • NRHS — INTEGER</p>. The number of right-hand-sides, i.e., the number of columns of the matrix B.
  • A — COMPLEX*16 array, dimension (LDA,N)</p>. On entry, the N-by-N matrix A. 🫶
  • LDA — INTEGER</p>. The leading dimension of the array A. LDA >= max(1,N).
  • AF — COMPLEX*16 array, dimension (LDAF,N)</p>. The block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by ZSYTRF. literally her whole backstory 💅
  • LDAF — INTEGER</p>. The leading dimension of the array AF. LDAF >= max(1,N).
  • IPIV — INTEGER array, dimension (N)</p>. Details of the interchanges and the block structure of D as determined by ZSYTRF.
  • COLEQU — LOGICAL</p>. If .TRUE. then column equilibration was done to A before calling this routine. This is needed to compute the solution and error bounds correctly. (don't lie to her about this one, she will notice 🫣)
  • C — DOUBLE PRECISION array, dimension (N)</p>. 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. ⚠️ (read that twice, bestie, that's the cautionary tale of the whole section)
  • B — COMPLEX*16 array, dimension (LDB,NRHS)</p>. The right-hand-side matrix B.
  • LDB — INTEGER</p>. The leading dimension of the array B. LDB >= max(1,N).
  • Y — COMPLEX*16 array, dimension (LDY,NRHS)</p>. On entry, the solution matrix X, as computed by ZSYTRS. On exit, the improved solution matrix Y. ✨ glow-up ✨
  • LDY — INTEGER</p>. The leading dimension of the array Y. LDY >= max(1,N).
  • BERR_OUT — DOUBLE PRECISION array, dimension (NRHS)</p>. 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 ZLA_LIN_BERR.

  • N_NORMS — INTEGER</p>. 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 — DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)</p>. 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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. 💀 (yeah the sentence just trails off in the source, we don't gatekeep, we report it verbatim)

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

Okay bestie, the parameter section is still going — these are the final function-doc entries, and they're giving workspace + tuning-knob realness 💅

  • ERR_BNDS_COMP — is a DOUBLE PRECISION 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:
                   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 < 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 — boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n) * slamch('Epsilon').
  • err = 2 — 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 . 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 — is COMPLEX*16 array, dimension (N). Workspace to hold the intermediate residual.
  • AYB — is DOUBLE PRECISION array, dimension (N). Workspace.
  • DY — is COMPLEX*16 array, dimension (N). Workspace to hold the intermediate solution.
  • Y_TAIL — is COMPLEX*16 array, dimension (N). Workspace to hold the trailing bits of the intermediate solution.
  • RCOND — is DOUBLE PRECISION. 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 — is 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 — is DOUBLE PRECISION. 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 — is DOUBLE PRECISION. Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solution Y is stable, which we define 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 — is LOGICAL. If .TRUE. then ignore componentwise convergence. Default value is .FALSE..
  • INFO — is INTEGER.
    • = 0: Successful exit.
    • < 0: if INFO = -i, the ith argument to ZLA_HERFSX_EXTENDED had an illegal value.

Author 💖

  • Univ. of Tennessee
  • Univ. of California Berkeley
  • Univ. of Colorado Denver
  • NAG Ltd.

Definition at line 388 of file zla_syrfsx_extended.f.

Author 🧾

Generated automatically by Doxygen for LAPACK from the source code.

That's the whole page, bestie 🫶 — official footer sign-off: LAPACK Version 3.12.0, manual page la_herfsx_extended(3). That's the source/version + date-stamp energy at the very bottom, and we are done slaying this manpage 💜✨

📦 where this page comes from