doubleGBcomputational(3)
π doubleGBcomputational(3) β LAPACK, but make it make sense π
okay bestie, welcome to the LAPACK corner of the universe where the matrices are banded and the vibes are immaculate. This is the manual page for the doubleGBcomputational squad β yes, that's the whole manual name, it's giving very "my full government name" energy and honestly? we respect it. π
π NAME π
doubleGBcomputational β it's the double precision era, and it's all about GB matrices (that's general band matrices, iykyk) π¦
β¨ SYNOPSIS β¨
π Functions π
dgbbrdβsubroutine dgbbrd (VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q, LDQ, PT, LDPT, C, LDC, WORK, INFO)β DGBBRDdgbconβsubroutine dgbcon (NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, WORK, IWORK, INFO)β DGBCONdgbequβsubroutine dgbequ (M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)β DGBEQUdgbequbβsubroutine dgbequb (M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)β DGBEQUBdgbrfsβsubroutine dgbrfs (TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)β DGBRFSdgbrfsxβsubroutine dgbrfsx (TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)β DGBRFSXdgbtf2βsubroutine dgbtf2 (M, N, KL, KU, AB, LDAB, IPIV, INFO)β DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.dgbtrfβsubroutine dgbtrf (M, N, KL, KU, AB, LDAB, IPIV, INFO)β DGBTRFdgbtrsβsubroutine dgbtrs (TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)β DGBTRSdggbakβsubroutine dggbak (JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)β DGGBAKdggbalβsubroutine dggbal (JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)β DGGBALdla_gbamvβsubroutine dla_gbamv (TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, INCX, BETA, Y, INCY)β DLA_GBAMV performs a matrix-vector operation to calculate error bounds.dla_gbrcondβdouble precision function dla_gbrcond (TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, IPIV, CMODE, C, INFO, WORK, IWORK)β DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.dla_gbrfsx_extendedβsubroutine dla_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, ERR_BNDS_NORM, ERR_BNDS_COMP, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)β DLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.dla_gbrpvgrwβdouble precision function dla_gbrpvgrw (N, KL, KU, NCOLS, AB, LDAB, AFB, LDAFB)β DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.dorgbrβsubroutine dorgbr (VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)β DORGBR
π«Ά Detailed Description π«Ά
This is the group of double computational functions for GB matrices β¨
π Function Documentation π
dgbbrd β DGBBRD π¦
subroutine dgbbrd (character VECT, integer M, integer N, integer NCC,
integer KL, integer KU, double precision, dimension( ldab, * ) AB,
integer LDAB, double precision, dimension( * ) D, double precision,
dimension( * ) E, double precision, dimension( ldq, * ) Q, integer LDQ,
double precision, dimension( ldpt, * ) PT, integer LDPT, double
precision, dimension( ldc, * ) C, integer LDC, double precision,
dimension( * ) WORK, integer INFO)
π Purpose:
DGBBRD reduces a real general m-by-n band matrix A to upper bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
The routine computes B, and optionally forms Q or P**T, or computes Q**T*C for a given matrix C.
π Parameters π
- VECT β
VECT is CHARACTER*1. Specifies whether or not the matrices Q and P**T are to be formed.= 'N': do not form Q or P**T;= 'Q': form Q only;= 'P': form P**T only;= 'B': form both.
- M β
M is INTEGER. The number of rows of the matrix A.M >= 0. - N β
N is INTEGER. The number of columns of the matrix A.N >= 0. - NCC β
NCC is INTEGER. The number of columns of the matrix C.NCC >= 0. - KL β
KL is INTEGER. The number of subdiagonals of the matrix A.KL >= 0. - KU β
KU is INTEGER. The number of superdiagonals of the matrix A.KU >= 0. - AB β
AB is DOUBLE PRECISION array, dimension (LDAB,N). On entry, the m-by-n band matrix A, stored in rows 1 toKL+KU+1. The j-th column of A is stored in the j-th column of the array AB as follows:AB(ku+1+i-j,j) = A(i,j)formax(1,j-ku)<=i<=min(m,j+kl). On exit, A is overwritten by values generated during the reduction. - LDAB β
LDAB is INTEGER. The leading dimension of the array A.LDAB >= KL+KU+1. - D β
D is DOUBLE PRECISION array, dimension (min(M,N)). The diagonal elements of the bidiagonal matrix B. - E β
E is DOUBLE PRECISION array, dimension (min(M,N)-1). The superdiagonal elements of the bidiagonal matrix B. - Q β
Q is DOUBLE PRECISION array, dimension (LDQ,M). IfVECT = 'Q'or'B', the m-by-m orthogonal matrix Q. IfVECT = 'N'or'P', the array Q is not referenced. - LDQ β
LDQ is INTEGER. The leading dimension of the array Q.LDQ >= max(1,M)ifVECT = 'Q'or'B';LDQ >= 1otherwise. - PT β
PT is DOUBLE PRECISION array, dimension (LDPT,N). IfVECT = 'P'or'B', the n-by-n orthogonal matrix P'. IfVECT = 'N'or'Q', the array PT is not referenced. - LDPT β
LDPT is INTEGER. The leading dimension of the array PT.LDPT >= max(1,N)ifVECT = 'P'or'B';LDPT >= 1otherwise. - C β
C is DOUBLE PRECISION array, dimension (LDC,NCC). On entry, an m-by-ncc matrix C. On exit, C is overwritten byQ**T*C. C is not referenced ifNCC = 0. - LDC β
LDC is INTEGER. The leading dimension of the array C.LDC >= max(1,M)ifNCC > 0;LDC >= 1ifNCC = 0. - WORK β
WORK is DOUBLE PRECISION array, dimension (2*max(M,N)) - INFO β
INFO is INTEGER= 0: successful exit.< 0: ifINFO = -i, the i-th argument had an illegal value.
πΈ Author πΈ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
dgbcon β DGBCON π
subroutine dgbcon (character NORM, integer N, integer KL, integer KU, double
precision, dimension( ldab, * ) AB, integer LDAB, integer, dimension( *
) IPIV, double precision ANORM, double precision RCOND, double
precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer
INFO)
π Purpose:
DGBCON estimates the reciprocal of the condition number of a real general band matrix A, in either the 1-norm or the infinity-norm, using the LU factorization computed by DGBTRF.
An estimate is obtained for norm(inv(A)), and the reciprocal of the condition number is computed as
RCOND = 1 / ( norm(A) * norm(inv(A)) ).
π Parameters π
- NORM β
NORM is CHARACTER*1. Specifies whether the 1-norm condition number or the infinity-norm condition number is required:= '1'or'O': 1-norm;= 'I': Infinity-norm.
- N β
N is INTEGER. The order of the matrix A.N >= 0. - KL β
KL is INTEGER. The number of subdiagonals within the band of A.KL >= 0. - KU β
KU is INTEGER. The number of superdiagonals within the band of A.KU >= 0. - AB β
AB is DOUBLE PRECISION array, dimension (LDAB,N). Details of the LU factorization of the band matrix A, as computed by DGBTRF. U is stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows 1 toKL+KU+1, and the multipliers used during the factorization are stored in rowsKL+KU+2to2*KL+KU+1. - LDAB β
LDAB is INTEGER. The leading dimension of the array AB.LDAB >= 2*KL+KU+1. - IPIV β
IPIV is INTEGER array, dimension (N). The pivot indices; for1 <= i <= N, row i of the matrix was interchanged with rowIPIV(i). - ANORM β
ANORM is DOUBLE PRECISION. IfNORM = '1'or'O', the 1-norm of the original matrix A. IfNORM = 'I', the infinity-norm of the original matrix A. - RCOND β
RCOND is DOUBLE PRECISION. The reciprocal of the condition number of the matrix A, computed asRCOND = 1/(norm(A) * norm(inv(A))). - WORK β
WORK is DOUBLE PRECISION array, dimension (3*N) - IWORK β
IWORK is INTEGER array, dimension (N) - INFO β
INFO is INTEGER= 0: successful exit< 0: ifINFO = -i, the i-th argument had an illegal value
πΈ Author πΈ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
dgbequ β DGBEQU π
subroutine dgbequ (integer M, integer N, integer KL, integer KU, double
precision, dimension( ldab, * ) AB, integer LDAB, double precision,
dimension( * ) R, double precision, dimension( * ) C, double precision
ROWCND, double precision COLCND, double precision AMAX, integer INFO)
π Purpose:
DGBEQU computes row and column scalings intended to equilibrate an M-by-N band matrix A and reduce its condition number. R returns the row scale factors and C the column scale factors, chosen to try to make the largest element in each row and column of the matrix B with elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
R(i) and C(j) are restricted to be between SMLNUM = smallest safe number and BIGNUM = largest safe number. Use of these scaling factors is not guaranteed to reduce the condition number of A but works well in practice.
π Parameters π
- M β
M is INTEGER. The number of rows of the matrix A.M >= 0. - N β
N is INTEGER. The number of columns of the matrix A.N >= 0. - KL β
KL is INTEGER. The number of subdiagonals within the band of A.KL >= 0. - KU β
KU is INTEGER. The number of subdiagonals within the band of A.KL >= 0. KUβKU is INTEGER, the number of superdiagonals living inside the band ofA. Constraint:KU >= 0. πABβAB is DOUBLE PRECISION array, dimension (LDAB,N). This is the band matrixA, chillin' in rows1toKL+KU+1. The j-th column ofAlives in the j-th column of the arrayAB, like this:
AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
LDABβLDAB is INTEGER, the leading dimension of the arrayAB. Constraint:LDAB >= KL+KU+1. πRβR is DOUBLE PRECISION array, dimension (M). IfINFO = 0, orINFO > M, thenRis holding the row scale factors forA. β¨CβC is DOUBLE PRECISION array, dimension (N). IfINFO = 0, thenCcontains the column scale factors forA.ROWCNDβROWCND is DOUBLE PRECISION. IfINFO = 0orINFO > M,ROWCNDcontains the ratio of the smallestR(i)to the largestR(i). And bestie, listen: ifROWCND >= 0.1andAMAXis neither too large nor too small, it is not worth scaling byR. periodt.COLCNDβCOLCND is DOUBLE PRECISION. IfINFO = 0,COLCNDcontains the ratio of the smallestC(i)to the largestC(i). IfCOLCND >= 0.1, it is not worth scaling byC. πAMAXβAMAX is DOUBLE PRECISION. This is the absolute value of the largest matrix element. IfAMAXis very close to overflow or very close to underflow, the matrix should be scaled. π¨INFOβINFO is INTEGER:= 0: successful exit β¨< 0: ifINFO = -i, the i-th argument had an illegal value π> 0: ifINFO = i, andiis:<= M: the i-th row ofAis exactly zero> M: the (i-M)-th column ofAis exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
subroutine dgbequb β DGBEQUB π
subroutine dgbequb (integer M, integer N, integer KL, integer KU, double
precision, dimension( ldab, * ) AB, integer LDAB, double precision,
dimension( * ) R, double precision, dimension( * ) C, double precision
ROWCND, double precision COLCND, double precision AMAX, integer INFO)
Purpose
DGBEQUB computes row and column scalings intended to equilibrate an M-by-N matrix A and reduce its condition number. R returns the row scale factors and C the column scale factors, chosen to try to make the largest element in each row and column of the matrix B with elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most the radix. π₯
R(i) and C(j) are restricted to be a power of the radix between SMLNUM = smallest safe number and BIGNUM = largest safe number. Use of these scaling factors is not guaranteed to reduce the condition number of A but works well in practice. no cap.
This routine differs from DGEEQU by restricting the scaling factors to a power of the radix. Barring over- and underflow, scaling by these factors introduces no additional rounding errors. However, the scaled entries' magnitudes are no longer approximately 1 but lie between sqrt(radix) and 1/sqrt(radix). π€―
Parameters
MβM is INTEGER, the number of rows of the matrixA. Constraint:M >= 0.NβN is INTEGER, the number of columns of the matrixA. Constraint:N >= 0.KLβKL is INTEGER, the number of subdiagonals within the band ofA. Constraint:KL >= 0.KUβKU is INTEGER, the number of superdiagonals within the band ofA. Constraint:KU >= 0. πABβAB is DOUBLE PRECISION array, dimension (LDAB,N). On entry, the matrixAin band storage, in rows1toKL+KU+1. The j-th column ofAis stored in the j-th column of the arrayABas follows:
AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
LDABβLDAB is INTEGER, the leading dimension of the arrayA. Constraint:LDAB >= max(1,M). πRβR is DOUBLE PRECISION array, dimension (M). IfINFO = 0orINFO > M,Rcontains the row scale factors forA. β¨CβC is DOUBLE PRECISION array, dimension (N). IfINFO = 0,Ccontains the column scale factors forA.ROWCNDβROWCND is DOUBLE PRECISION. IfINFO = 0orINFO > M,ROWCNDcontains the ratio of the smallestR(i)to the largestR(i). IfROWCND >= 0.1andAMAXis neither too large nor too small, it is not worth scaling byR.COLCNDβCOLCND is DOUBLE PRECISION. IfINFO = 0,COLCNDcontains the ratio of the smallestC(i)to the largestC(i). IfCOLCND >= 0.1, it is not worth scaling byC. πAMAXβAMAX is DOUBLE PRECISION. Absolute value of largest matrix element. IfAMAXis very close to overflow or very close to underflow, the matrix should be scaled. π¨INFOβINFO is INTEGER:= 0: successful exit β¨< 0: ifINFO = -i, the i-th argument had an illegal value π> 0: ifINFO = i, andiis:<= M: the i-th row ofAis exactly zero> M: the (i-M)-th column ofAis exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
subroutine dgbrfs β DGBRFS π¦
subroutine dgbrfs (character TRANS, integer N, integer KL, integer KU,
integer NRHS, double precision, dimension( ldab, * ) AB, integer LDAB,
double precision, dimension( ldafb, * ) AFB, integer LDAFB, integer,
dimension( * ) IPIV, double precision, dimension( ldb, * ) B, integer
LDB, double precision, dimension( ldx, * ) X, integer LDX, double
precision, dimension( * ) FERR, double precision, dimension( * ) BERR,
double precision, dimension( * ) WORK, integer, dimension( * ) IWORK,
integer INFO)
Purpose
DGBRFS improves the computed solution to a system of linear equations when the coefficient matrix is banded, and provides error bounds and backward error estimates for the solution. it's giving accuracy glow-up. π
Parameters
TRANSβTRANS is 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)
NβN is INTEGER, the order of the matrixA. Constraint:N >= 0.KLβKL is INTEGER, the number of subdiagonals within the band ofA. Constraint:KL >= 0.KUβKU is INTEGER, the number of superdiagonals within the band ofA. Constraint:KU >= 0.NRHSβNRHS is INTEGER, the number of right hand sides, i.e., the number of columns of the matricesBandX. Constraint:NRHS >= 0.ABβAB is DOUBLE PRECISION array, dimension (LDAB,N). The original band matrixA, stored in rows1toKL+KU+1. The j-th column ofAis stored in the j-th column of the arrayABas follows:
AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
LDABβLDAB is INTEGER, the leading dimension of the arrayAB. Constraint:LDAB >= KL+KU+1. πAFBβAFB is DOUBLE PRECISION array, dimension (LDAFB,N). Details of the LU factorization of the band matrixA, as computed byDGBTRF.Uis stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows1toKL+KU+1, and the multipliers used during the factorization are stored in rowsKL+KU+2to2*KL+KU+1. β¨LDAFBβLDAFB is INTEGER, the leading dimension of the arrayAFB. Constraint:LDAFB >= 2*KL*KU+1.IPIVβIPIV is INTEGER array, dimension (N). The pivot indices fromDGBTRF; for1<=i<=N, rowiof the matrix was interchanged with rowIPIV(i).BβB is DOUBLE PRECISION array, dimension (LDB,NRHS). The right hand side matrixB.LDBβLDB is INTEGER, the leading dimension of the arrayB. Constraint:LDB >= max(1,N).XβX is DOUBLE PRECISION array, dimension (LDX,NRHS). On entry, the solution matrixX, as computed byDGBTRS. On exit, the improved solution matrixX. πLDXβLDX is INTEGER, the leading dimension of the arrayX. Constraint:LDX >= max(1,N).FERRβFERR is DOUBLE PRECISION array, dimension (NRHS). The estimated forward error bound for each solution vectorX(j)(the j-th column of the solution matrixX). IfXTRUEis the true solution corresponding toX(j),FERR(j)is an estimated upper bound for the magnitude of the largest element in(X(j) - XTRUE)divided by the magnitude of the largest element inX(j). The estimate is as reliable as the estimate forRCOND, and is almost always a slight overestimate of the true error. πBERRβBERR is DOUBLE PRECISION array, dimension (NRHS). The componentwise relative backward error of each solution vectorX(j)(i.e., the smallest relative change in any element ofAorBthat makesX(j)an exact solution).WORKβWORK is DOUBLE PRECISION array, dimension (3*N)IWORKβIWORK is INTEGER array, dimension (N)INFOβINFO is INTEGER:= 0: successful exit β¨< 0: ifINFO = -i, the i-th argument had an illegal value π
Internal Parameters:
ITMAX is the maximum number of steps of iterative refinement. π₯
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π subroutine dgbrfsx
okay bestie so this one is a subroutine, meaning she doesn't hand you back a value β she just does the work and slides the results into the arguments you gave her. Here's the full squad, character-for-character, exactly how the compiler needs to see them:
subroutine dgbrfsx (character TRANS, character EQUED, integer N, integer KL,
integer KU, integer NRHS, double precision, dimension( ldab, * ) AB,
integer LDAB, double precision, dimension( ldafb, * ) AFB, integer
LDAFB, integer, dimension( * ) IPIV, double precision, dimension( * ) R,
double precision, dimension( * ) C, double precision, dimension( ldb, *
) B, integer LDB, double precision, dimension( ldx , * ) X, integer LDX,
double precision RCOND, double precision, dimension( * ) BERR, integer
N_ERR_BNDS, double precision, dimension( nrhs, * ) ERR_BNDS_NORM, double
precision, dimension( nrhs, * ) ERR_BNDS_COMP, integer NPARAMS, double
precision, dimension( * ) PARAMS, double precision, dimension( * ) WORK,
integer, dimension( * ) IWORK, integer INFO)
β¨ Purpose β¨
DGBRFSX is the glow-up routine. She takes the solution you already computed to a system of linear equations, makes it better, and then hands you receipts: error bounds and backward error estimates for the solution. And she's not one-note about it β on top of the normwise error bound, she'll also give you the maximum componentwise error bound whenever that's possible. (Peep the comments for ERR_BNDS_NORM and ERR_BNDS_COMP if you want the full tea on how those bounds work.)
Now, plot twist: your original system of linear equations might have already been equilibrated before you ever called her β that's the whole EQUED, R and C situation described down below. If that's what happened, don't stress: the solution and error bounds she returns are for the original, unequilibrated system. She un-does the glow-up math for you. π
And there's a lil' bundle of optional settings living inside the PARAMS array. Those knobs decide how the refinement gets performed β but honestly? The defaults are usually giving exactly what you need. If the defaults are fine for you, you can just pass NPARAMS = 0, which stops the source code from even touching the PARAMS argument. Lazy girl math but make it correct. π
π Parameters π
TRANSβTRANS is CHARACTER*1. Tells us which form of the system of equations we're solving:= 'N':A * X = B(no transpose)= 'T':A**T * X = B(transpose)= 'C':A**H * X = B(conjugate transpose = transpose)
EQUEDβEQUED is CHARACTER*1. Which flavor of equilibration got applied toAbefore this routine was called. She genuinely needs this to compute the solution and the error bounds correctly, so don't lie to her:= 'N': no equilibration= 'R': row equilibration, i.e.,Ahas been premultiplied bydiag(R).= 'C': column equilibration, i.e.,Ahas been postmultiplied bydiag(C).= 'B': both row and column equilibration, i.e.,Ahas been replaced bydiag(R) * A * diag(C). The right hand sideBhas been changed accordingly.
NβN is INTEGER. The order of the matrixA. Constraint:N >= 0.KLβKL is INTEGER. How many subdiagonals live inside the band ofA. Constraint:KL >= 0.KUβKU is INTEGER. How many superdiagonals live inside the band ofA. Constraint:KU >= 0.NRHSβNRHS is INTEGER. The number of right hand sides, aka the number of columns of the matricesBandX. Constraint:NRHS >= 0.ABβAB is DOUBLE PRECISION array, dimension (LDAB,N). The original band matrixA, stored in rows1toKL+KU+1. The j-th column ofAis stored in the j-th column of the arrayABlike so:
AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
LDABβLDAB is INTEGER. The leading dimension of the arrayAB. Constraint:LDAB >= KL+KU+1.AFBβAFB is DOUBLE PRECISION array, dimension (LDAFB,N). The details of the LU factorization of the band matrixA, exactly as computed byDGBTRF.Uis stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows1toKL+KU+1, and the multipliers that got used during the factorization are parked in rowsKL+KU+2to2*KL+KU+1.LDAFBβLDAFB is INTEGER. The leading dimension of the arrayAFB. Constraint:LDAFB >= 2*KL*KU+1.IPIVβIPIV is INTEGER array, dimension (N). The pivot indices fromDGETRF; for1<=i<=N, rowiof the matrix was interchanged with rowIPIV(i).RβR is DOUBLE PRECISION array, dimension (N). The row scale factors forA. IfEQUED = 'R'or'B',Ais multiplied on the left bydiag(R); ifEQUED = 'N'or'C',Ris not accessed at all.Ris an input argument ifFACT = 'F'; otherwise,Ris an output argument. IfFACT = 'F'andEQUED = 'R'or'B', every element ofRmust be positive. IfRis output, each element ofRis a power of the radix. IfRis input, each element ofRshould 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.CβC is DOUBLE PRECISION array, dimension (N). The column scale factors forA. IfEQUED = 'C'or'B',Ais multiplied on the right bydiag(C); ifEQUED = 'N'or'R',Cis not accessed at all.Cis an input argument ifFACT = 'F'; otherwise,Cis an output argument. IfFACT = 'F'andEQUED = 'C'or'B', every element ofCmust be positive. IfCis output, each element ofCis a power of the radix. IfCis input, each element ofCshould 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.
lowkey important: those "should be a power of the radix" lines onRandCare not decoration, bestie. Ignore them and your error estimates may stop being trustworthy. π
BβB is DOUBLE PRECISION array, dimension (LDB,NRHS). The right hand side matrixB.LDBβLDB is INTEGER. The leading dimension of the arrayB. Constraint:LDB >= max(1,N).XβX is DOUBLE PRECISION array, dimension (LDX,NRHS). On entry, the solution matrixX, as computed byDGETRS. On exit, the improved solution matrixX. β¨ glow-up complete.LDXβLDX is INTEGER. The leading dimension of the arrayX. Constraint:LDX >= max(1,N).RCONDβRCOND is DOUBLE PRECISION. Reciprocal scaled condition number. This is an estimate of the reciprocal Skeel condition number of the matrixAafter equilibration (if any happened). If this comes back less than the machine precision β and yes, in particular if it's zero β then the matrix is singular to working precision. But heads up: even if this number is tiny and the matrix is looking all ill-conditioned, the error can still be small. Don't panic-read this one. π«ΆBERRβBERR is DOUBLE PRECISION array, dimension (NRHS). Componentwise relative backward error. That's the componentwise relative backward error of each solution vectorX(j)β i.e., the smallest relative change in any element ofAorBthat makesX(j)an exact solution.N_ERR_BNDSβN_ERR_BNDS is INTEGER. How many error bounds to return for each right hand side and each type (normwise or componentwise). SeeERR_BNDS_NORMandERR_BNDS_COMPright below. πERR_BNDS_NORMβERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS). For each right-hand side, this array is stuffed with info about the various error bounds and condition numbers corresponding to the normwise relative error, which is 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 the following three fields:
err = 1β the 'Trust/don't trust' boolean. Trust the answer if the reciprocal condition number is less than the thresholdsqrt(n) * dlamch('Epsilon').err = 2β the 'Guaranteed' error bound: the estimated forward error, almost certainly within a factor of10of the true error, so long as the next entry is greater than the thresholdsqrt(n) * dlamch('Epsilon'). This error bound should only be trusted if the previous boolean is true.err = 3β Reciprocal condition number: the estimated normwise reciprocal condition number. This gets compared with the thresholdsqrt(n) * dlamch('Epsilon')to determine whether the error estimate is 'guaranteed'. These reciprocal condition numbers are1 / (norm(Z^{-1},inf) * norm(Z,inf))for some appropriately scaled matrixZ. LetZ = S*A, whereSscales each row by a power of the radix so all absolute row sums ofZare approximately1.
See Lapack Working Note 165 for further details and extra cautions.
ERR_BNDS_COMPβERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS). For each right-hand side, this array carries info about the various error bounds and condition numbers corresponding to the componentwise relative error, which is defined like this:
Componentwise relative error in the ith solution vector:
abs(XTRUE(j,i) - X(j,i))
max_j ----------------------
abs(X(j,i))
okay so <pink>ERR_BNDS_COMP</pink> is indexed by the right-hand side i β the side where the componentwise relative error actually depends on things β and then by the type of error info, which we're about to break all the way down. There are currently up to three pieces of info returned for each right-hand side. If componentwise accuracy was never requested (that's PARAMS(3) = 0.0), then <hotpink>ERR_BNDS_COMP</hotpink> is not accessed at all. And if N_ERR_BNDS < 3, then at most the first (:,N_ERR_BNDS) entries are returned. That's the whole deal, no cap.
The first index in ERR_BNDS_COMP(i,:) corresponds to the ith right-hand side. Simple, cute, iconic.
The second index in ERR_BNDS_COMP(:,err) contains the following three fields:
err = 1β the<purple>'Trust/don't trust'</purple>boolean. Trust the answer if the reciprocal condition number is less than the thresholdsqrt(n) * dlamch('Epsilon').err = 2β the<purple>'Guaranteed'</purple>error bound: the estimated forward error, almost certainly within a factor of10of the true error so long as the next entry is greater than the thresholdsqrt(n) * dlamch('Epsilon'). And lowkey β this error bound should only be trusted if the previous boolean is true. Otherwise bestie it is giving nothing.err = 3β reciprocal condition number: the estimated componentwise reciprocal condition number. You compare it with the thresholdsqrt(n) * dlamch('Epsilon')to determine if the error estimate is'guaranteed'. These reciprocal condition numbers are1 / (norm(Z^{-1},inf) * norm(Z,inf))for some appropriately scaled matrixZ. LetZ = S*(A*diag(x)), wherexis the solution for the current right-hand side andSscales each row ofA*diag(x)by a power of the radix so all absolute row sums ofZare approximately1.
π Want the full lore plus extra cautions? See Lapack Working Note 165. That's the receipts, fr fr.
π« NPARAMS
NPARAMS is INTEGER. It specifies the number of parameters set in PARAMS. If <= 0, the PARAMS array is never referenced and default values are used. Periodt.
π« PARAMS
PARAMS is DOUBLE PRECISION array, dimension (NPARAMS). It specifies algorithm parameters. If an entry is < 0.0, then that entry will be filled with the default value used for that parameter. Only positions up to NPARAMS are accessed; defaults are used for higher-numbered parameters.
PARAMS(LA_LINRX_ITREF_I = 1)β whether to perform iterative refinement or not.- Default:
1.0D+0 = 0.0: no refinement is performed, and no error bounds are computed.= 1.0: use the double-precision refinement algorithm, possibly with doubled-single computations if the compilation environment does not supportDOUBLE PRECISION.- (other values are reserved for future use)
- Default:
PARAMS(LA_LINRX_ITHRESH_I = 2)β maximum number of residual computations allowed for refinement.- Default:
10 - Aggressive: set to
100to permit convergence using approximate factorizations or factorizations other than LU. π¨ If the factorization uses a technique other than Gaussian elimination, the guarantees inerr_bnds_normanderr_bnds_compmay no longer be trustworthy.
- Default:
PARAMS(LA_LINRX_CWISE_I = 3)β 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.0is false.- Default:
1.0(attempt componentwise convergence)
- Default:
π« WORK
WORK is DOUBLE PRECISION array, dimension (4*N)
π« IWORK
IWORK is INTEGER array, dimension (N)
π« INFO
INFO is INTEGER
= 0: successful exit. The solution to every right-hand side is guaranteed. π< 0: ifINFO = -i, thei-th argument had an illegal value. Rude, but okay.> 0and<= N:U(INFO,INFO)is exactly zero. The factorization has been completed, but the factorUis exactly singular, so the solution and error bounds could not be computed.RCOND = 0is returned.= N+J: the solution corresponding to theJth right-hand side is not guaranteed. The solutions corresponding to other right-hand sidesKwithK > Jmay not be guaranteed as well, but only the first such right-hand side is reported. If a small componentwise error is not requested (PARAMS(3) = 0.0) then theJth right-hand side is the first with a normwise error bound that is not guaranteed (the smallestJsuch thatERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) theJth right-hand side is the first with either a normwise or componentwise error bound that is not guaranteed (the smallestJsuch that eitherERR_BNDS_NORM(J,1) = 0.0orERR_BNDS_COMP(J,1) = 0.0). See the definition ofERR_BNDS_NORM(:,1)andERR_BNDS_COMP(:,1). To get information about all of the right-hand sides, checkERR_BNDS_NORMorERR_BNDS_COMP.
π Author π
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
π¦ subroutine dgbtf2 π¦
subroutine dgbtf2 (integer M, integer N, integer KL, integer KU, double
precision, dimension( ldab, * ) AB, integer LDAB, integer, dimension( *
) IPIV, integer INFO)
DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm. Yes, the unblocked one, we're going vintage.
β¨ Purpose: β¨
DGBTF2 computes an LU factorization of a real m-by-n band matrix A using partial pivoting with row interchanges.
This is the unblocked version of the algorithm, calling Level 2 BLAS.
β¨ Parameters β¨
MβMisINTEGER. The number of rows of the matrixA.M >= 0.NβNisINTEGER. The number of columns of the matrixA.N >= 0.KLβKLisINTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβKUisINTEGER. The number of superdiagonals within the band ofA.KU >= 0.ABβABisDOUBLE PRECISION array, dimension (LDAB,N). On entry, the matrixAin band storage, in rowsKL+1to2*KL+KU+1; rows1toKLof the array need not be set. Thej-th column ofAis stored in thej-th column of the arrayABas follows:
AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
On exit, details of the factorization: U is stored as an upper triangular band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and the multipliers used during the factorization are stored in rows KL+KU+2 to 2*KL+KU+1. See below for further details.
LDABβLDABisINTEGER. The leading dimension of the arrayAB.LDAB >= 2*KL+KU+1.IPIVβIPIVisINTEGER array, dimension (min(M,N)). The pivot indices; for1 <= i <= min(M,N), rowiof the matrix was interchanged with rowIPIV(i).INFOβINFOisINTEGER= 0: successful exit< 0: ifINFO = -i, thei-th argument had an illegal value> 0: ifINFO = +i,U(i,i)is exactly zero. The factorization has been completed, but the factorUis exactly singular, and division by zero will occur if it is used to solve a system of equations. π
β¨ Author β¨
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
β¨ Further Details: β¨
The band storage scheme is illustrated by the following example, when M = N = 6, KL = 2, KU = 1:
On entry: On exit:
* * * + + + * * * u14 u25 u36
* * + + + + * * u13 u24 u35 u46
* a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
a31 a42 a53 a64 * * m31 m42 m53 m64 * *
Array elements marked * are not used by the routine; elements marked + need not be set on entry, but are required by the routine to store elements of U, because of fill-in resulting from the row interchanges.
π¦ subroutine dgbtrf π¦
subroutine dgbtrf (integer M, integer N, integer KL, integer KU, double
precision, dimension( ldab, * ) AB, integer LDAB, integer, dimension( *
) IPIV, integer INFO)
DGBTRF β that's literally the whole one-line description, bestie. No notes.
β¨ Purpose: β¨
DGBTRF computes an LU factorization of a real m-by-n band matrix A using partial pivoting with row interchanges.
This is the blocked version of the algorithm, calling Level 3 BLAS. Glow up from its unblocked sibling. β¨
β¨ Parameters β¨
MβMisINTEGER. The number of rows of the matrixA.M >= 0.NβNisINTEGER. The number of columns of the matrixA.N >= 0.KLβKLisINTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβKUisINTEGER. The number of superdiagonals within the band ofA.KU >= 0.ABβABisDOUBLE PRECISION array, dimension (LDAB,N). On entry, the matrixAin band storage, in rowsKL+1to2*KL+KU+1; rows1toKLof the array need not be set. Thej-th column ofAis stored in thej-th column of the arrayABas follows:
AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
On exit, details of the factorization: U is stored as an upper triangular band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and the multipliers used during the factorization are stored in rows KL+KU+2 to 2*KL+KU+1. See below for further details.
LDABLDABβLDABis INTEGER. It's the leading dimension of the arrayAB, and bestie,LDAB >= 2*KL+KU+1. That's the rule, no cap. πIPIVβIPIVis INTEGER array, dimension(min(M,N)). These are the pivot indices: for1 <= i <= min(M,N), row i of the matrix was interchanged with rowIPIV(i). She's been shuffled around, andIPIVremembers every single swap. π₯INFOβINFOis INTEGER, and it's the vibe check on how everything went:= 0: successful exit. Slay. β¨< 0: ifINFO = -i, the i-th argument had an illegal value. Not it. π> 0: ifINFO = +i,U(i,i)is exactly zero. The factorization has been completed, but the factorUis exactly singular, and division by zero will occur if it is used to solve a system of equations. π«
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π¦ Further Details:
The band storage scheme is illustrated by the following example, when M = N = 6, KL = 2, KU = 1:
On entry: On exit:
* * * + + + * * * u14 u25 u36
* * + + + + * * u13 u24 u35 u46
* a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
a31 a42 a53 a64 * * m31 m42 m53 m64 * *
Array elements marked * are not used by the routine; elements marked + need not be set on entry, but are required by the routine to store elements of U because of fill-in resulting from the row interchanges. Okay so the * girlies are just chillin', untouched, and the + slots are empty at first but the routine is gonna need them later for the fill-in glow-up. πΈ
β¨ subroutine dgbtrs β¨
subroutine dgbtrs (character TRANS, integer N, integer KL, integer KU,
integer NRHS, double precision, dimension( ldab, * ) AB, integer LDAB,
integer, dimension( * ) IPIV, double precision, dimension( ldb, * ) B,
integer LDB, integer INFO)
DGBTRS
π Purpose:
DGBTRS solves a system of linear equations
A * X = B or A**T * X = B
with a general band matrix A using the LU factorization computed by DGBTRF. She's the second half of the duo β DGBTRF does the factoring, and DGBTRS actually solves the thing. Dynamic duo behavior. π«Ά
π Parameters
TRANSβTRANSisCHARACTER*1. Specifies the form of the system of equations.= 'N':A * X = B(No transpose)= 'T':A**T* X = B(Transpose)= 'C':A**T* X = B(Conjugate transpose = Transpose)
NβNis INTEGER. The order of the matrixA.N >= 0.KLβKLis INTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβKUis INTEGER. The number of superdiagonals within the band ofA.KU >= 0.NRHSβNRHSis INTEGER. The number of right hand sides, i.e., the number of columns of the matrixB.NRHS >= 0.ABβABis DOUBLE PRECISION array, dimension(LDAB,N). Details of the LU factorization of the band matrixA, as computed byDGBTRF.Uis stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows1toKL+KU+1, and the multipliers used during the factorization are stored in rowsKL+KU+2to2*KL+KU+1. It's all stacked in there, iconic organization honestly. πLDABβLDABis INTEGER. The leading dimension of the arrayAB.LDAB >= 2*KL+KU+1.IPIVβIPIVis INTEGER array, dimension(N). The pivot indices; for1 <= i <= N, row i of the matrix was interchanged with rowIPIV(i).BβBis DOUBLE PRECISION array, dimension(LDB,NRHS). On entry, the right hand side matrixB. On exit, the solution matrixX. She goes in as the problem and comes out as the answer, main character energy fr. π₯LDBβLDBis INTEGER. The leading dimension of the arrayB.LDB >= max(1,N).INFOβINFOis INTEGER.= 0: successful exit β¨< 0: ifINFO = -i, the i-th argument had an illegal value π
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π subroutine dggbak π
subroutine dggbak (character JOB, character SIDE, integer N, integer ILO,
integer IHI, double precision, dimension( * ) LSCALE, double precision,
dimension( * ) RSCALE, integer M, double precision, dimension( ldv, * )
V, integer LDV, integer INFO)
DGGBAK
π Purpose:
DGGBAK forms the right or left eigenvectors of a real generalized eigenvalue problem A*x = lambda*B*x, by backward transformation on the computed eigenvectors of the balanced pair of matrices output by DGGBAL. Basically: DGGBAL did the balancing, and now DGGBAK undoes it so your eigenvectors actually mean something in the original problem. Undo button, but make it math. π§
π Parameters
JOBβJOBisCHARACTER*1. Specifies the type of backward transformation required:= 'N': do nothing, return immediately;= 'P': do backward transformation for permutation only;= 'S': do backward transformation for scaling only;= 'B': do backward transformations for both permutation and scaling.
JOB must be the same as the argument JOB supplied to DGGBAL. Like, they have to match, okay? No mismatched energy allowed here. π
SIDEβSIDEisCHARACTER*1.= 'R':Vcontains right eigenvectors;= 'L':Vcontains left eigenvectors.
NβNis INTEGER. The number of rows of the matrixV.N >= 0.ILOβILOis INTEGERIHIβIHIis INTEGER
The integers ILO and IHI determined by DGGBAL. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
LSCALEβLSCALEis DOUBLE PRECISION array, dimension(N). Details of the permutations and/or scaling factors applied to the left side ofAandB, as returned byDGGBAL.RSCALEβRSCALEis DOUBLE PRECISION array, dimension(N). Details of the permutations and/or scaling factors applied to the right side ofAandB, as returned byDGGBAL.MβMis INTEGER. The number of columns of the matrixV.M >= 0.VβVis DOUBLE PRECISION array, dimension(LDV,M). On entry, the matrix of right or left eigenvectors to be transformed, as returned byDTGEVC. On exit,Vis overwritten by the transformed eigenvectors. Glow-up complete. β¨LDVβLDVis INTEGER. The leading dimension of the matrixV.LDV >= max(1,N).INFOβINFOis INTEGER.= 0: successful exit.< 0: ifINFO = -i, the i-th argument had an illegal value.
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π¦ Further Details:
See R.C. Ward, Balancing the generalized eigenvalue problem, SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
π subroutine dggbal π
subroutine dggbal (character JOB, integer N, double precision, dimension(
lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B,
integer LDB, integer ILO, integer IHI, double precision, dimension( * )
LSCALE, double precision, dimension( * ) RSCALE, double precision,
dimension( * ) WORK, integer INFO)
DGGBAL
π Purpose:
DGGBAL balances a pair of general real matrices (A,B). This involves, first, permuting A and B by similarity transformations to isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N elements on the diagonal; and second, applying a diagonal similarity transformation to rows and columns ILO to IHI to make the rows and columns as close in norm as possible. Both steps are optional. She's out here tidying up your matrices so the eigenvalue problem behaves. So thoughtful. π
Balancing may reduce the 1-norm of the matrices, and improve the accuracy of the computed eigenvalues and/or eigenvectors in the generalized eigenvalue problem A*x = lambda*B*x. So like, it's optional but it's also lowkey the reason your answers are accurate. Just saying. π
π Parameters
JOBβJOBisCHARACTER*1. Specifies the operations to be performed onAandB:= 'N': none: simply setILO = 1,IHI = N,LSCALE(I) = 1.0andRSCALE(I) = 1.0fori = 1,...,N.= 'P': permute only;= 'S': scale only;= 'B': both permute and scale.
NβNis INTEGER. The order of the matricesAandB.N >= 0.AβAis DOUBLE PRECISION array, dimension(LDA,N). On entry, the input matrixA. On exit,Ais overwritten by the balanced matrix. IfJOB = 'N',Ais not referenced. She didn't even look at her, periodt. πLDAβLDAis INTEGER. The leading dimension of the arrayA.LDA >= max(1,N).BβBis DOUBLE PRECISION array, dimension(LDB,N). On entry, the input matrixB. On exit,Bis overwritten by the balanced matrix. IfJOB = 'N',Bis not referenced.LDBβLDBis INTEGER. The leading dimension of the arrayB.LDB >= max(1,N).ILOβILOis INTEGERIHIβIHIis INTEGER
ILO and IHI are set to integers such that on exit A(i,j) = 0 and B(i,j) = 0 if i > j and j = 1,...,ILO-1 or i = IHI+1,...,N. If JOB = 'N' or 'S', ILO = 1 and IHI = N.
LSCALEβLSCALEis DOUBLE PRECISION array, dimension(N). Details of the permutations and scaling factors applied to the left side ofAandB. IfP(j)is the index of the row interchanged with rowj, andD(j)is the scaling factor applied to rowj, then
LSCALE(j) = P(j) for J = 1,...,ILO-1
= D(j) for J = ILO,...,IHI
= P(j) for J = IHI+1,...,N.
The order in which the interchanges are made is N to IHI+1, then 1 to ILO-1. Order matters bestie, don't come for her sequencing. π
RSCALERSCALEβ it's aDOUBLE PRECISIONarray, dimension(N). She holds the receipts on the permutations and scaling factors that got applied to the right side ofAandB. So ifP(j)is the index of the column that got swapped with columnj, andD(j)is the scaling factor that got applied to columnj, then bestie, here's the tea:
LSCALE(j) = P(j) for J = 1,...,ILO-1
= D(j) for J = ILO,...,IHI
= P(j) for J = IHI+1,...,N.
And the order in which those interchanges go down? It's N to IHI+1, then 1 to ILO-1. That's the sequence, don't fight it. π
WORKβWORKis aDOUBLE PRECISIONarray, dimension(lwork). Now listen:lworkhas to be at leastmax(1,6*N)whenJOB = 'S'or'B', and at least1whenJOB = 'N'or'P'. Don't skimp on the workspace, she needs her room. β¨INFOβINFOisINTEGER, and she's the vibe check at the end:= 0: successful exit, we ate, no notes. π< 0: ifINFO = -i, then the i-th argument had an illegal value. Somebody fumbled the input, not us.
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π Further Details:
See R.C. WARD, Balancing the generalized eigenvalue problem, SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
π¦ subroutine dla_gbamv
subroutine dla_gbamv (integer TRANS, integer M, integer N, integer KL,
integer KU, double precision ALPHA, double precision, dimension( ldab, *
) AB, integer LDAB, double precision, dimension( * ) X, integer INCX,
double precision BETA, double precision, dimension( * ) Y, integer INCY)
DLA_GBAMV does a matrix-vector operation to calculate error bounds. That's literally her whole job, and she's so good at it. π«Ά
π― Purpose:
So DLA_GBAMV performs one of these matrix-vector operations:
y := alpha*abs(A)*abs(x) + beta*abs(y),
or y := alpha*abs(A)**T*abs(x) + beta*abs(y),
where alpha and beta are scalars, x and y are vectors, and A is an m by n matrix.
This function is primarily used in calculating error bounds. And here's the underflow protection lore, bestie: to protect against underflow during evaluation, components in the resulting vector get perturbed away from zero by (N+1) times the underflow threshold. And to prevent unnecessarily large errors for block-structure embedded in general matrices, the 'symbolically' zero components do not get perturbed. A zero entry is considered 'symbolic' if all the multiplications involved in computing that entry have at least one zero multiplicand. iykyk π
βοΈ Parameters
TRANSβTRANSisINTEGER. On entry,TRANSspecifies the operation to be performed, as follows:
BLAS_NO_TRANS y := alpha*abs(A)*abs(x) + beta*abs(y)
BLAS_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
BLAS_CONJ_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
Unchanged on exit. Periodt. β¨
MβMisINTEGER. On entry,Mspecifies the number of rows of the matrixA.Mmust be at least zero. Unchanged on exit.NβNisINTEGER. On entry,Nspecifies the number of columns of the matrixA.Nmust be at least zero. Unchanged on exit.KLβKLisINTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβKUisINTEGER. The number of superdiagonals within the band ofA.KU >= 0.ALPHAβALPHAisDOUBLE PRECISION. On entry,ALPHAspecifies the scalaralpha. Unchanged on exit.ABβABis aDOUBLE PRECISIONarray, dimension( LDAB, n ). Before entry, the leadingm by npart of the arrayABmust contain the matrix of coefficients. Unchanged on exit.LDABβLDABisINTEGER. On entry,LDAspecifies the first dimension ofABas declared in the calling (sub) program.LDABmust be at leastmax( 1, m ). Unchanged on exit.XβXis aDOUBLE PRECISIONarray, dimension( 1 + ( n - 1 )*abs( INCX ) )whenTRANS = 'N'or'n', and at least( 1 + ( m - 1 )*abs( INCX ) )otherwise. Before entry, the incremented arrayXmust contain the vectorx. Unchanged on exit.INCXβINCXisINTEGER. On entry,INCXspecifies the increment for the elements ofX.INCXmust not be zero. Unchanged on exit.BETAβBETAisDOUBLE PRECISION. On entry,BETAspecifies the scalarbeta. WhenBETAis supplied as zero thenYneed not be set on input. Unchanged on exit.YβYis aDOUBLE PRECISIONarray, dimension( 1 + ( m - 1 )*abs( INCY ) )whenTRANS = 'N'or'n', and at least( 1 + ( n - 1 )*abs( INCY ) )otherwise. Before entry withBETAnon-zero, the incremented arrayYmust contain the vectory. On exit,Yis overwritten by the updated vectory. πINCYβINCYisINTEGER. On entry,INCYspecifies the increment for the elements ofY.INCYmust not be zero. Unchanged on exit.
Level 2 Blas routine.
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
πΈ double precision function dla_gbrcond
double precision function dla_gbrcond (character TRANS, integer N, integer
KL, integer KU, double precision, dimension( ldab, * ) AB, integer LDAB,
double precision, dimension( ldafb, * ) AFB, integer LDAFB, integer,
dimension( * ) IPIV, integer CMODE, double precision, dimension( * ) C,
integer INFO, double precision, dimension( * ) WORK, integer, dimension(
* ) IWORK)
DLA_GBRCOND estimates the Skeel condition number for a general banded matrix. Iconic behaviour honestly. β¨
π― Purpose:
DLA_GBRCOND Estimates the Skeel condition number of op(A) * op2(C)
where op2 is determined by CMODE as follows
CMODE = 1 op2(C) = C
CMODE = 0 op2(C) = I
CMODE = -1 op2(C) = inv(C)
The Skeel condition number cond(A) = norminf( |inv(A)||A| )
is computed by computing scaling factors R such that
diag(R)*A*op2(C) is row equilibrated and computing the standard
infinity-norm condition number.
βοΈ Parameters
TRANSβTRANSisCHARACTER*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)
NβNisINTEGER. The number of linear equations, i.e., the order of the matrixA.N >= 0.KLβKLisINTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβKUisINTEGER. The number of superdiagonals within the band ofA.KU >= 0.ABβABis aDOUBLE PRECISIONarray, dimension(LDAB,N). On entry, the matrixAin band storage, in rows1toKL+KU+1. The j-th column ofAis stored in the j-th column of the arrayABas follows:
AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
LDABβLDABisINTEGER. The leading dimension of the arrayAB.LDAB >= KL+KU+1.AFBβAFBis aDOUBLE PRECISIONarray, dimension(LDAFB,N). Details of the LU factorization of the band matrixA, as computed byDGBTRF.Uis stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows1toKL+KU+1, and the multipliers used during the factorization are stored in rowsKL+KU+2to2*KL+KU+1.LDAFBβLDAFBisINTEGER. The leading dimension of the arrayAFB.LDAFB >= 2*KL+KU+1.IPIVβIPIVis anINTEGERarray, dimension(N). The pivot indices from the factorizationA = P*L*Uas computed byDGBTRF; rowiof the matrix was interchanged with rowIPIV(i).CMODEβCMODEisINTEGER. Determinesop2(C)in the formulaop(A) * op2(C)as follows:CMODE = 1 op2(C) = CCMODE = 0 op2(C) = ICMODE = -1 op2(C) = inv(C)
CβCis aDOUBLE PRECISIONarray, dimension(N). The vectorCin the formulaop(A) * op2(C).INFOβINFOisINTEGER.= 0: Successful exit.i > 0: The ith argument is invalid. π
WORKβWORKis aDOUBLE PRECISIONarray, dimension(5*N). Workspace, bestie.IWORKβIWORKis anINTEGERarray, dimension(N). More workspace. π
π Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
π§ subroutine dla_gbrfsx_extended
subroutine dla_gbrfsx_extended (integer PREC_TYPE, integer TRANS_TYPE,
integer N, integer KL, integer KU, integer NRHS, double precision,
dimension( ldab, * ) AB, integer LDAB, double precision, dimension(
ldafb, * ) AFB, integer LDAFB, 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_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded matrices by performing extra-precise iterative refinement, and she also provides error bounds and backward error estimates for the solution. Multitasking queen. π
π― Purpose:
DLA_GBRFSX_EXTENDED is the glow-up routine β¨ it improves the computed solution to a system of linear equations by doing extra-precise iterative refinement, and it hands you error bounds plus backward error estimates for the solution. DGBRFSX is the one that calls it to run that iterative refinement. On top of the normwise error bound, the code also gives you the maximum componentwise error bound when it possibly can β peek the comments for ERR_BNDS_NORM and ERR_BNDS_COMP for all the deets on those bounds. But real talk: this subroutine is only responsible for setting the second fields of ERR_BNDS_NORM and ERR_BNDS_COMP. That's its whole job, no cap π
PREC_TYPEβINTEGER. Specifies the intermediate precision to be used in refinement. The value is defined byILAPREC(P)wherePis aCHARACTER, and:P = 'S': SingleP = 'D': DoubleP = 'I': IndigenousP = 'X'orP = 'E': Extra
TRANS_TYPEβINTEGER. Specifies the transposition operation onA. The value is defined byILATRANS(T)whereTis aCHARACTER, and:T = 'N': No transposeT = 'T': TransposeT = 'C': Conjugate transpose
NβINTEGER. The number of linear equations, aka the order of the matrixA.N >= 0.KLβINTEGER. The number of subdiagonals within the band ofA.KL >= 0.KUβINTEGER. The number of superdiagonals within the band ofA.KU >= 0NRHSβINTEGER. The number of right-hand-sides, i.e., the number of columns of the matrixB.ABβDOUBLE PRECISIONarray, dimension(LDAB,N). On entry, theN-by-NmatrixAB.LDABβINTEGER. The leading dimension of the arrayAB.LDBA >= max(1,N).AFBβDOUBLE PRECISIONarray, dimension(LDAFB,N). The factorsLandUfrom the factorizationA = P*L*Uas computed byDGBTRF.LDAFBβINTEGER. The leading dimension of the arrayAF.LDAFB >= max(1,N).IPIVβINTEGERarray, dimension(N). The pivot indices from the factorizationA = P*L*Uas computed byDGBTRF; rowiof the matrix was interchanged with rowIPIV(i).COLEQUβLOGICAL. If.TRUE.then column equilibration was done toAbefore calling this routine. This is needed to compute the solution and error bounds correctly β no shortcuts bestie πCβDOUBLE PRECISIONarray, dimension(N). The column scale factors forA. IfCOLEQU = .FALSE.,Cis not accessed. IfCis input, each element ofCshould 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 like, keep it to powers of the radix, okay?BβDOUBLE PRECISIONarray, dimension(LDB,NRHS). The right-hand-side matrixB.LDBβINTEGER. The leading dimension of the arrayB.LDB >= max(1,N).YβDOUBLE PRECISIONarray, dimension(LDY,NRHS). On entry, the solution matrixX, as computed byDGBTRS. On exit, the improved solution matrixY. the glow-up happened πLDYβINTEGER. The leading dimension of the arrayY.LDY >= max(1,N).BERR_OUTβDOUBLE PRECISIONarray, dimension(NRHS). On exit,BERR_OUT(j)contains the componentwise relative backward error for right-hand-sidejfrom 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.
N_NORMSβINTEGER. Determines which error bounds to return (seeERR_BNDS_NORMandERR_BNDS_COMP).- If
N_NORMS >= 1return normwise error bounds. - If
N_NORMS >= 2return componentwise error bounds.
- If
ERR_BNDS_NORMβDOUBLE PRECISIONarray, 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β the 'Trust/don't trust' boolean. Trust the answer if the reciprocal condition number is less than the thresholdsqrt(n) * slamch('Epsilon').err = 2β the 'Guaranteed' error bound: the estimated forward error, almost certainly within a factor of10of the true error so long as the next entry is greater than the thresholdsqrt(n) * slamch('Epsilon'). This error bound should only be trusted if the previous boolean is true.err = 3β the reciprocal condition number: estimated normwise reciprocal condition number. Compared with the thresholdsqrt(n) * slamch('Epsilon')to determine if the error estimate is 'guaranteed'. These reciprocal condition numbers are1 / (norm(Z^{-1},inf) * norm(Z,inf))for some appropriately scaled matrixZ. LetZ = S*A, whereSscales each row by a power of the radix so all absolute row sums ofZare approximately1.
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 PRECISIONarray, 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β the 'Trust/don't trust' boolean. Trust the answer if the reciprocal condition number is less than the thresholdsqrt(n) * slamch('Epsilon').err = 2β the 'Guaranteed' error bound: the estimated forward error, almost certainly within a factor of10of the true error so long as the next entry is greater than the thresholdsqrt(n) * slamch('Epsilon'). This error bound should only be trusted if the previous boolean is true.err = 3β the reciprocal condition number: estimated componentwise reciprocal condition number. Compared with the thresholdsqrt(n) * slamch('Epsilon')to determine if the error estimate is 'guaranteed'. These reciprocal condition numbers are1 / (norm(Z^{-1},inf) * norm(Z,inf))for some appropriately scaled matrixZ. LetZ = S*(A*diag(x)), wherexis the solution for the current right-hand side andSscales each row ofA*diag(x)by a power of the radix so all absolute row sums ofZare approximately1.
This subroutine is only responsible for setting the second field above.
See Lapack Working Note 165 for further details and extra cautions. π
RESβDOUBLE PRECISIONarray, dimension(N). Workspace to hold the intermediate residual.AYBβDOUBLE PRECISIONarray, dimension(N). Workspace. This can be the same workspace passed forY_TAIL.DYβDOUBLE PRECISIONarray, dimension(N). Workspace to hold the intermediate solution.Y_TAILβDOUBLE PRECISIONarray, dimension(N). Workspace to hold the trailing bits of the intermediate solution.RCONDβDOUBLE PRECISION. Reciprocal scaled condition number. This is an estimate of the reciprocal Skeel condition number of the matrixAafter 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. π don't panic just bc it's tiny, okay?ITHRESHITHRESHβ ITHRESH is INTEGER π The max number of residual computations allowed for refinement, bestie. Default is10. For 'aggressive' set it to100to permit convergence using approximate factorizations or factorizations other than LU. But lowkey be careful: if the factorization uses a technique other than Gaussian elimination, the guarantees inERR_BNDS_NORMandERR_BNDS_COMPmay no longer be trustworthy. πRTHRESHβ RTHRESH is DOUBLE PRECISION β¨ Determines when to stop refinement if the error estimate stops decreasing. Refinement will stop when the next solution no longer satisfiesnorm(dx_{i+1}) < RTHRESH * norm(dx_i)wherenorm(Z)is the infinity norm ofZ.RTHRESHsatisfies0 < RTHRESH <= 1. Default is0.5. For 'aggressive' set to0.9to permit convergence on extremely ill-conditioned matrices. See LAWN 165 for more details. πDZ_UBβ DZ_UB is DOUBLE PRECISION π¦ Determines when to start considering componentwise convergence. Componentwise convergence is only considered after each component of the solutionYis stable, which we define as the relative change in each component being less thanDZ_UB. Default value is0.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.. Periodt.INFOβ INFO is INTEGER π₯= 0: Successful exit. Yay!< 0: ifINFO = -i, the ith argument toDGBTRShad an illegal value. Not iconic.
π Author π
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
π« double precision function dla_gbrpvgrw π«
double precision function dla_gbrpvgrw (integer N, integer KL, integer KU,
integer NCOLS, double precision, dimension( ldab, * ) AB, integer LDAB,
double precision, dimension( ldafb, * ) AFB, integer LDAFB)
DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix. Slay.
π― Purpose π―
DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U). The 'max absolute element' norm is used. If this is much less than 1, 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. Big yikes. π
β¨ Parameters β¨
Nβ N is INTEGER π The number of linear equations, i.e., the order of the matrixA.N >= 0.KLβ KL is INTEGER πΈ The number of subdiagonals within the band ofA.KL >= 0.KUβ KU is INTEGER π¦ The number of superdiagonals within the band ofA.KU >= 0.NCOLSβ NCOLS is INTEGER π The number of columns of the matrixA.NCOLS >= 0.ABβ AB is DOUBLE PRECISION array, dimension (LDAB,N) π₯ On entry, the matrixAin band storage, in rows1toKL+KU+1. The j-th column ofAis stored in the j-th column of the arrayABas follows:
AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
Keep it exact, bestie.
LDABβ LDAB is INTEGER π The leading dimension of the arrayAB.LDAB >= KL+KU+1.AFBβ AFB is DOUBLE PRECISION array, dimension (LDAFB,N) π Details of the LU factorization of the band matrixA, as computed byDGBTRF.Uis stored as an upper triangular band matrix withKL+KUsuperdiagonals in rows1toKL+KU+1, and the multipliers used during the factorization are stored in rowsKL+KU+2to2*KL+KU+1.LDAFBβ LDAFB is INTEGER β The leading dimension of the arrayAFB.LDAFB >= 2*KL+KU+1.
π Author π
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
π subroutine dorgbr π
subroutine dorgbr (character VECT, integer M, integer N, integer K, double
precision, dimension( lda, * ) A, integer LDA, double precision,
dimension( * ) TAU, double precision, dimension( * ) WORK, integer
LWORK, integer INFO)
DORGBR
π― Purpose π―
DORGBR generates one of the real orthogonal matrices Q or P**T determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T are defined as products of elementary reflectors H(i) or G(i) respectively. Main character energy.
If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M:
- if
m >= k,Q = H(1) H(2) . . . H(k)andDORGBRreturns the firstncolumns ofQ, wherem >= n >= k; - if
m < k,Q = H(1) H(2) . . . H(m-1)andDORGBRreturnsQas an M-by-M matrix.
If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T is of order N:
- if
k < n,P**T = G(k) . . . G(2) G(1)andDORGBRreturns the firstmrows ofP**T, wheren >= m >= k; - if
k >= n,P**T = G(n-1) . . . G(2) G(1)andDORGBRreturnsP**Tas an N-by-N matrix.
β¨ Parameters β¨
VECTβ VECT is CHARACTER1* π Specifies whether the matrixQor the matrixP**Tis required, as defined in the transformation applied byDGEBRD:= 'Q': generateQ;= 'P': generateP**T.
Mβ M is INTEGER πΈ The number of rows of the matrixQorP**Tto be returned.M >= 0.Nβ N is INTEGER π¦ The number of columns of the matrixQorP**Tto be returned.N >= 0. IfVECT = 'Q',M >= N >= min(M,K); ifVECT = 'P',N >= M >= min(N,K).Kβ K is INTEGER π IfVECT = 'Q', the number of columns in the original M-by-K matrix reduced byDGEBRD. IfVECT = 'P', the number of rows in the original K-by-N matrix reduced byDGEBRD.K >= 0.Aβ A is DOUBLE PRECISION array, dimension (LDA,N) π₯ On entry, the vectors which define the elementary reflectors, as returned byDGEBRD. On exit, the M-by-N matrixQorP**T.LDAβ LDA is INTEGER π The leading dimension of the arrayA.LDA >= max(1,M).TAUβ TAU is DOUBLE PRECISION array, dimension π(min(M,K))ifVECT = 'Q'(min(N,K))ifVECT = 'P'TAU(i)must contain the scalar factor of the elementary reflectorH(i)orG(i), which determinesQorP**T, as returned byDGEBRDin its array argumentTAUQorTAUP.
WORKβ WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) β On exit, ifINFO = 0,WORK(1)returns the optimalLWORK.LWORKβ LWORK is INTEGER π The dimension of the arrayWORK.LWORK >= max(1,min(M,N)). For optimum performanceLWORK >= min(M,N)*NB, whereNBis the optimal blocksize. IfLWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of theWORKarray, returns this value as the first entry of theWORKarray, and no error message related toLWORKis issued byXERBLA. No cap.INFOβ INFO is INTEGER π= 0: successful exit< 0: ifINFO = -i, the i-th argument had an illegal value
π Author π
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
π€ Author π€
Generated automatically by Doxygen for LAPACK from the source code.
And that's the tea, bestie! This whole iconic page is from LAPACK Version 3.11.0, dated Sun Nov 27 2022. Stay slay. π β¨π
π¦ where this page comes from
- debian bookworm: liblapack-doc 3.11.0-2 Β· BSD-3-Clause
- π girlified by deepseek-ai/DeepSeek-V4.1-Flash on : a language model's rewrite of the original, under its terms