Developer Reference

Contents

?gbsvxx

Uses extra precise iterative refinement to compute the solution to the system of linear equations with a banded coefficient matrix A and multiple right-hand sides

Syntax

lapack_int LAPACKE_sgbsvxx
(
int
matrix_layout
,
char
fact
,
char
trans
,
lapack_int
n
,
lapack_int
kl
,
lapack_int
ku
,
lapack_int
nrhs
,
float*
ab
,
lapack_int
ldab
,
float*
afb
,
lapack_int
ldafb
,
lapack_int*
ipiv
,
char*
equed
,
float*
r
,
float*
c
,
float*
b
,
lapack_int
ldb
,
float*
x
,
lapack_int
ldx
,
float*
rcond
,
float*
rpvgrw
,
float*
berr
,
lapack_int
n_err_bnds
,
float*
err_bnds_norm
,
float*
err_bnds_comp
,
lapack_int
nparams
,
const float*
params
);
lapack_int LAPACKE_dgbsvxx
(
int
matrix_layout
,
char
fact
,
char
trans
,
lapack_int
n
,
lapack_int
kl
,
lapack_int
ku
,
lapack_int
nrhs
,
double*
ab
,
lapack_int
ldab
,
double*
afb
,
lapack_int
ldafb
,
lapack_int*
ipiv
,
char*
equed
,
double*
r
,
double*
c
,
double*
b
,
lapack_int
ldb
,
double*
x
,
lapack_int
ldx
,
double*
rcond
,
double*
rpvgrw
,
double*
berr
,
lapack_int
n_err_bnds
,
double*
err_bnds_norm
,
double*
err_bnds_comp
,
lapack_int
nparams
,
const double*
params
);
lapack_int LAPACKE_cgbsvxx
(
int
matrix_layout
,
char
fact
,
char
trans
,
lapack_int
n
,
lapack_int
kl
,
lapack_int
ku
,
lapack_int
nrhs
,
lapack_complex_float*
ab
,
lapack_int
ldab
,
lapack_complex_float*
afb
,
lapack_int
ldafb
,
lapack_int*
ipiv
,
char*
equed
,
float*
r
,
float*
c
,
lapack_complex_float*
b
,
lapack_int
ldb
,
lapack_complex_float*
x
,
lapack_int
ldx
,
float*
rcond
,
float*
rpvgrw
,
float*
berr
,
lapack_int
n_err_bnds
,
float*
err_bnds_norm
,
float*
err_bnds_comp
,
lapack_int
nparams
,
const float*
params
);
lapack_int LAPACKE_zgbsvxx
(
int
matrix_layout
,
char
fact
,
char
trans
,
lapack_int
n
,
lapack_int
kl
,
lapack_int
ku
,
lapack_int
nrhs
,
lapack_complex_double*
ab
,
lapack_int
ldab
,
lapack_complex_double*
afb
,
lapack_int
ldafb
,
lapack_int*
ipiv
,
char*
equed
,
double*
r
,
double*
c
,
lapack_complex_double*
b
,
lapack_int
ldb
,
lapack_complex_double*
x
,
lapack_int
ldx
,
double*
rcond
,
double*
rpvgrw
,
double*
berr
,
lapack_int
n_err_bnds
,
double*
err_bnds_norm
,
double*
err_bnds_comp
,
lapack_int
nparams
,
const double*
params
);
Include Files
  • mkl.h
Description
The routine uses the LU factorization to compute the solution to a real or complex system of linear equations
A*X
=
B
,
A
T
*
X
=
B
, or
A
H
*
X
=
B
, where
A
is an
n
-by-
n
banded matrix, the columns of the matrix
B
are individual right-hand sides, and the columns of
X
are the corresponding solutions.
Both normwise and maximum componentwise error bounds are also provided on request. The routine returns a solution with a small guaranteed error (
O(eps)
, where
eps
is the working machine precision) unless the matrix is very ill-conditioned, in which case a warning is returned. Relevant condition numbers are also calculated and returned.
The routine accepts user-provided factorizations and equilibration factors; see definitions of the
fact
and
equed
options. Solving with refinement and using a factorization from a previous call of the routine also produces a solution with
O(eps)
errors or warnings but that may not be true for general user-provided factorizations and equilibration factors if they differ from what the routine would itself produce.
The routine
?gbsvxx
performs the following steps:
  1. If
    fact
    =
    'E'
    , scaling factors
    r
    and
    c
    are computed to equilibrate the system:
    trans
    =
    'N'
    :
    diag
    (
    r
    )*
    A
    *
    diag
    (
    c
    )*inv(
    diag
    (
    c
    ))*
    X
    =
    diag
    (
    r
    )*
    B
    trans
    =
    'T'
    :
    (
    diag
    (
    r
    )*
    A
    *
    diag
    (
    c
    ))
    T
    *inv(
    diag
    (
    r
    ))*
    X
    =
    diag
    (
    c
    )*
    B
    trans
    =
    'C'
    :
    (
    diag
    (
    r
    )*
    A
    *
    diag
    (
    c
    ))
    H
    *inv(
    diag
    (
    r
    ))*
    X
    =
    diag
    (
    c
    )*
    B
    Whether or not the system will be equilibrated depends on the scaling of the matrix
    A
    , but if equilibration is used,
    A
    is overwritten by
    diag
    (
    r
    )*
    A
    *
    diag
    (
    c
    )
    and
    B
    by
    diag
    (
    r
    )*
    B
    (if
    trans
    =
    'N'
    )
    or
    diag
    (
    c
    )*
    B
    (if
    trans
    =
    'T'
    or
    'C'
    ).
  2. If
    fact
    =
    'N'
    or
    'E'
    , the
    LU
    decomposition is used to factor the matrix
    A
    (after equilibration if
    fact
    =
    'E'
    ) as
    A
    =
    P*L*U
    , where
    P
    is a permutation matrix,
    L
    is a unit lower triangular matrix, and
    U
    is upper triangular.
  3. If some
    U
    i
    ,
    i
    = 0, so that
    U
    is exactly singular, then the routine returns with
    info
    =
    i
    . Otherwise, the factored form of
    A
    is used to estimate the condition number of the matrix
    A
    (see the
    rcond
    parameter). If the reciprocal of the condition number is less than machine precision, the routine still goes on to solve for
    X
    and compute error bounds.
  4. The system of equations is solved for
    X
    using the factored form of
    A
    .
  5. By default, unless
    params[0]
    is set to zero, the routine applies iterative refinement to improve the computed solution matrix and calculate error bounds. Refinement calculates the residual to at least twice the working precision.
  6. If equilibration was used, the matrix
    X
    is premultiplied by
    diag
    (
    c
    )
    (if
    trans
    =
    'N'
    ) or
    diag
    (
    r
    )
    (if
    trans
    =
    'T'
    or
    'C'
    ) so that it solves the original system before equilibration.
Input Parameters
matrix_layout
Specifies whether two-dimensional array storage is row major (
LAPACK_ROW_MAJOR
) or column major (
LAPACK_COL_MAJOR
).
fact
Must be
'F'
,
'N'
, or
'E'
.
Specifies whether or not the factored form of the matrix
A
is supplied on entry, and if not, whether the matrix
A
should be equilibrated before it is factored.
If
fact
=
'F'
, on entry,
afb
and
ipiv
contain the factored form of
A
. If
equed
is not
'N'
, the matrix
A
has been equilibrated with scaling factors given by
r
and
c
. Parameters
ab
,
afb
, and
ipiv
are not modified.
If
fact
=
'N'
, the matrix
A
will be copied to
afb
and factored.
If
fact
=
'E'
, the matrix
A
will be equilibrated, if necessary, copied to
afb
and factored.
trans
Must be
'N'
,
'T'
, or
'C'
.
Specifies the form of the system of equations:
If
trans
=
'N'
, the system has the form
A
*
X
=
B
(No transpose).
If
trans
=
'T'
, the system has the form
A
T
*
X
=
B
(Transpose).
If
trans
=
'C'
, the system has the form
A
H
*
X
=
B
(Conjugate Transpose = Transpose for real flavors, Conjugate Transpose for complex flavors).
n
The number of linear equations; the order of the matrix
A
;
n
0.
kl
The number of subdiagonals within the band of
A
;
kl
0.
ku
The number of superdiagonals within the band of
A
;
ku
0.
nrhs
The number of right-hand sides; the number of columns of the matrices
B
and
X
;
nrhs
0.
ab
,
afb
,
b
Arrays:
ab
(max(
ldab
*
n
))
,
afb
(max(
ldafb
*
n
))
,
b
(max(1,
ldb
*
nrhs
) for column major layout and max(1,
ldb
*
n
) for row major layout)
.
The array
ab
contains the matrix
A
in band storage.
If
fact
=
'F'
and
equed
is not
'N'
, then
AB
must have been equilibrated by the scaling factors in
r
and/or
c
.
The array
afb
is an input argument if
fact
=
'F'
. It contains the factored form of the banded matrix
A
, that is, the factors
L
and
U
from the factorization
A
=
P*L*U
as computed by
?gbtrf
.
U
is stored as an upper triangular banded matrix with
kl
+
ku
superdiagonals.
L
is stored as lower triangular band matrix with
kl
subdiagonals. If
equed
is not
'N'
, then
afb
is the factored form of the equilibrated matrix
A
.
The array
b
contains the matrix
B
whose columns are the right-hand sides for the systems of equations.
ldab
The leading dimension of the array
ab
;
ldab
kl
+
ku
+1
.
ldafb
The leading dimension of the array
afb
;
ldafb
2*
kl
+
ku
+1
.
ipiv
Array, size at least
max(1,
n
)
. The array
ipiv
is an input argument if
fact
=
'F'
. It contains the pivot indices from the factorization
A
=
P*L*U
as computed by
?gbtrf
; row
i
of the matrix was interchanged with row
ipiv
[
i
-1]
.
equed
Must be
'N'
,
'R'
,
'C'
, or
'B'
.
equed
is an input argument if
fact
=
'F'
. It specifies the form of equilibration that was done:
If
equed
=
'N'
, no equilibration was done (always true if
fact
=
'N'
).
If
equed
=
'R'
, row equilibration was done, that is,
A
has been premultiplied by
diag
(
r
).
If
equed
=
'C'
, column equilibration was done, that is,
A
has been postmultiplied by
diag
(
c
).
If
equed
=
'B'
, both row and column equilibration was done, that is,
A
has been replaced by
diag
(
r
)*
A
*
diag
(
c
)
.
r
,
c
Arrays:
r
(size
n
),
c
(size
n
). The array
r
contains the row scale factors for
A
, and the array
c
contains the column scale factors for
A
. These arrays are input arguments if
fact
=
'F'
only; otherwise they are output arguments.
If
equed
=
'R'
or
'B'
,
A
is multiplied on the left by
diag
(
r
); if
equed
=
'N'
or
'C'
,
r
is not accessed.
If
fact
=
'F'
and
equed
=
'R'
or
'B'
, each element of
r
must be positive.
If
equed
=
'C'
or
'B'
,
A
is multiplied on the right by
diag
(
c
); if
equed
=
'N'
or
'R'
,
c
is not accessed.
If
fact
=
'F'
and
equed
=
'C'
or
'B'
, each element of
c
must be positive.
Each element of
r
or
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.
ldb
The leading dimension of the array
b
;
ldb
max(1,
n
) for column major layout and
ldb
nrhs
for row major layout
.
ldx
The leading dimension of the output array
x
;
ldx
max(1,
n
) for column major layout and
ldx
nrhs
for row major layout
.
n_err_bnds
Number of error bounds to return for each right hand side and each type (normwise or componentwise). See
err_bnds_norm
and
err_bnds_comp
descriptions in
Output Arguments
section below.
nparams
Specifies the number of parameters set in
params
. If
0, the
params
array is never referenced and default values are used.
params
Array, size
max(1,
nparams)
. Specifies algorithm parameters. If an entry is less than 0.0, that entry is filled with the default value used for that parameter. Only positions up to
nparams
are accessed; defaults are used for higher-numbered parameters. If defaults are acceptable, you can pass
nparams
= 0, which prevents the source code from accessing the
params
argument.
params
[0]
: Whether to perform iterative refinement or not. Default: 1.0 (for single precision flavors), 1.0D+0 (for double precision flavors).
=0.0
No refinement is performed and no error bounds are computed.
=1.0
Use the extra-precise refinement algorithm.
(Other values are reserved for future use.)
params
[1]
: Maximum number of residual computations allowed for refinement.
Default
10.0
Aggressive
Set to 100.0 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.
params
[2]
: Flag determining if the code will attempt to find a solution with a small componentwise relative error in the double-precision algorithm. Positive is true, 0.0 is false. Default: 1.0 (attempt componentwise convergence).
Output Parameters
x
Array, size
max(1,
ldx
*
nrhs
) for column major layout and max(1,
ldx
*
n
) for row major layout
.
If
info
= 0
, the array
x
contains the solution
n
-by-
nrhs
matrix
X
to the
original
system of equations. Note that
A
and
B
are modified on exit if
equed
'N'
, and the solution to the equilibrated system is:
inv(
diag
(
c
))*
X
, if
trans
=
'N'
and
equed
=
'C'
or
'B'
; or
inv(
diag
(
r
))*
X
, if
trans
=
'T'
or
'C'
and
equed
=
'R'
or
'B'
.
ab
Array
ab
is not modified on exit if
fact
=
'F'
or
'N'
, or if
fact
=
'E'
and
equed
=
'N'
.
If
equed
'N'
,
A
is scaled on exit as follows:
equed
=
'R'
:
A
=
diag
(
r
)*
A
equed
=
'C'
:
A
=
A
*
diag
(
c
)
equed
=
'B'
:
A
=
diag
(
r
)*
A
*
diag
(
c
)
.
afb
If
fact
=
'N'
or
'E'
, then
afb
is an output argument and on exit returns the factors
L
and
U
from the factorization
A
=
PLU
of the original matrix
A
(if
fact
=
'N'
) or of the equilibrated matrix
A
(if
fact
=
'E'
).
b
Overwritten by
diag
(
r
)*
B
if
trans
=
'N'
and
equed
=
'R'
or
'B'
;
overwritten by
trans
=
'T'
or
'C'
and
equed
=
'C'
or
'B'
;
not changed if
equed
=
'N'
.
r
,
c
These arrays are output arguments if
fact
'F'
. Each element of these arrays is a power of the radix. See the description of
r
,
c
in
Input Arguments
section.
rcond
Reciprocal scaled condition number. An estimate of the reciprocal Skeel condition number of the matrix
A
after equilibration (if done). If
rcond
is less than the machine precision, in particular, if
rcond
= 0, 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.
rpvgrw
Contains the reciprocal pivot growth factor:
If this is much less than 1, the stability of the
LU
factorization of the (equlibrated) matrix
A
could be poor. This also means that the solution
X
, estimated condition numbers, and error bounds could be unreliable. If factorization fails with
0 <
info
n
, this parameter contains the reciprocal pivot growth factor for the leading
info
columns of
A
. In
?gbsvx
, this quantity is returned in
rpivot
.
berr
Array, size at least
max(1,
nrhs
)
. Contains the componentwise relative backward error for each solution vector
x
j
, that is, the smallest relative change in any element of
A
or
B
that makes
x
j
an exact solution.
err_bnds_norm
Array of size
nrhs
*
n_err_bnds
. For each right-hand side, 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
i
-th solution vector
The array is indexed by the type of error information as described below. There are currently up to three pieces of information returned.
err
=1
"Trust/don't trust" boolean. Trust the answer if the reciprocal condition number is less than the threshold
sqrt
(
n
)*
slamch
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors.
err
=2
"Guaranteed" error bound. The estimated forward error, almost certainly within a factor of 10 of the true error so long as the next entry is greater than the threshold
sqrt
(
n
)*
slamch
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors. 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
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors to determine if the error estimate is "guaranteed". These reciprocal condition numbers for some appropriately scaled matrix
Z
are:
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.
The information for right-hand side
i
, where 1
i
nrhs
, and type of error
err
is stored in
err_bnds_norm
[(
err
-1)*
nrhs
+
i
- 1]
.
err_bnds_comp
Array of size
nrhs
*
n_err_bnds
. For each right-hand side, 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
i
-th solution vector:
The array is indexed by the type of error information as described below. There are currently up to three pieces of information returned for each right-hand side. If componentwise accuracy is not requested (
params[2]
= 0.0), then
err_bnds_comp
is not accessed.
err
=1
"Trust/don't trust" boolean. Trust the answer if the reciprocal condition number is less than the threshold
sqrt
(
n
)*
slamch
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors.
err
=2
"Guaranteed" error bpound. 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
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors. 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
(
ε
)
for single precision flavors and
sqrt
(
n
)*
dlamch
(
ε
)
for double precision flavors to determine if the error estimate is "guaranteed". These reciprocal condition numbers for some appropriately scaled matrix
Z
are:
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.
The information for right-hand side
i
, where 1
i
nrhs
, and type of error
err
is stored in
err_bnds_comp
[(
err
-1)*
nrhs
+
i
- 1]
.
ipiv
If
fact
=
'N'
or