CXML
duskyr
Unsymmetric sparse iterative refinement using skyline storage
scheme
FORMAT
DUSKYR
(n, au, auf, iaudiag, nau, al, alf, ialdiag, nal, b, ldb, x, ldx, ferr,
berr, nbx, iparam, rparam, iwrk, rwrk, ierror)
Arguments
n integer*4
On entry, the order of the matrix A.
On exit, n is unchanged.
au real*8
On entry, an array containing information on the matrix
A. If istore = 1 or 2, then au contains the upper
triangular part, including the diagonal, of the matrix
A, stored in the profile-in or diagonal-out mode,
respectively. Array AU is of length at least nau, where
nau is the envelope size of the upper triangular part
of A, including the diagonal. If istore = 3, then au
contains the matrix A, stored in the structurally
symmetric, profile-in storage mode. In this case, array
AU is of length at least nau, where nau is the envelope
size of the matrix A.
On exit, au is unchanged.
auf real*8
On entry, if istore = 1 or 2, auf contains the
factors U and D of the L*D*U factorization of the
matrix A. Array AUF is of length at least nau, where
nau is the envelope size of the upper triangular part
of A, including the diagonal. If istore = 3, then auf
contains the L*D*U factorization of the matrix A. In
this case, array AUF is of length at least nau, where
nau is the envelope size of the matrix A. The L*D*U
factorization has been obtained by a prior call to the
routine DUSKYF. auf must remain unchanged between
calls to the routines DUSKYF and DUSKYR.
On exit, auf is unchanged.
iaudiag integer*4
On entry, an array containing the pointers to the
locations of the diagonal elements in the arrays AU and
AUF. iaudiag is of length at least n for the profile-
in and the structurally symmetric profile-in storage
modes. iaudiag is of length at least (n+1) for the
diagonal-out storage mode.
On exit, iaudiag is unchanged.
nau integer*4
On entry, the number of elements stored in array AU.
If istore = 1 or 2, then nau is the envelope size of
the upper triangular part of the matrix A. If istore =
3, then nau is the envelope size of the matrix A. For
the profile-in and the structurally symmetric profile-
in storage modes, nau = IAUDIAG(n). For the
diagonal-out storage mode, nau = IAUDIAG(n+1) - 1.
On exit, nau is unchanged.
al real*8
On entry, an array containing information on the matrix
A. If istore = 1 or 2, then al contains the lower
triangular part, including the diagonal, of the matrix
A, stored in the profile-in or diagonal-out mode,
respectively. Storage is allocated for the diagonal
elements, though the elements themselves are not
stored. Array AL is of length at least nal, where nal
is the envelope size of the lower triangular part of A,
including the diagonal. If istore = 3, then al is a
dummy argument.
On exit, al is unchanged.
alf real*8
On entry, if istore = 1 or 2, alf contains the factor
L of the L*D*U factorization of the matrix A. Array ALF
is of length at least nal, where nal is the envelope
size of the lower triangular part of A, including the
diagonal. If istore = 3, then alf is a dummy
argument. The L*D*U factorization is obtained from a
prior call to the routine DUSKYF. alf must remain
unchanged between calls to the routines DUSKYF and
DUSKYR.
On exit, alf is unchanged.
ialdiag integer*4
On entry, an array containing the pointers to the
locations of the diagonal elements in the arrays AL and
ALF. ialdiag is of length at least n for the profile-
in storage mode. ialdiag is of length at least (n+1)
for the diagonal-out storage mode. If istore = 3,
then ialdiag is a dummy argument.
On exit, ialdiag is unchanged.
nal integer*4
On entry, the number of elements stored in array AL.
If istore = 1 or 2, then nal is the envelope size of
the lower triangular part of the matrix A. For the
profile-in storage mode, nal = IALDIAG(n). For the
diagonal-out storage mode, nal = IALDIAG(n+1) - 1. If
istore = 3, then nal is a dummy argument.
On exit, nal is unchanged.
b real*8
On entry, a two dimensional array B of order ldb by at
least nbx, containing the nbx right sides.
On exit, b is unchanged.
ldb integer*4
On entry, the leading dimension of array B.
On exit, ldb is unchanged.
x real*8
On entry, a two dimensional array X of order ldx by at
least nbx, containing the nbx solution vectors obtained
after a call to the routine DUSKYS.
On exit, x contains the improved solutions obtained
after iterative refinement.
ldx integer*4
On entry, the leading dimension of array X. ldx >=n.
On exit, ldx is unchanged.
ferr real*8
On entry, an array FERR of length at least nbx, whose
elements are unspecified variables.
On exit, ferr contains the estimated error bounds for
each of the nbz solution vectors.
berr real*8
On entry, an array BERR of length at least nbx, whose
elements are unspecified variables.
On exit, berr contains the component-wise relative
backward error for each of the nbz solution vectors.
nbx integer*4
On entry, the number of right sides.
On exit, nbx is unchanged.
iparam integer*4
An array of length at least 100, containing the integer
parameters for the iterative refinement and error
bounds calculation.
iparam(1): niparam
On entry, defines the length of the array IPARAM.
niparam >= 100.
On exit, iparam(1) is unchanged.
iparam(2): nrparam
On entry, defines the length of the array RPARAM. As
the real parameters array is not used at present,
nrparam may be unspecified.
On exit, iparam(2) is unchanged.
iparam(3): niwrk
On entry, defines the size of the integer work array,
IWRK. niwrk >=5n.
On exit, iparam(3) is unchanged.
iparam(4): nrwrk
On entry, defines the size of the real work array,
RWRK. nrwrk >=3n.
On exit, iparam(4) is unchanged.
iparam(5): iounit
On entry, defines the I/O unit number for printing
error messages and information from the routine DUSKYR.
The I/O unit must be opened in the calling subprogram.
If iounit <= 0, no output is generated.
On exit, iparam(5) is unchanged.
iparam(6): iolevel
On entry, defines the message level that determines the
amount of information printed out to iounit, when
iounit > 0.
iolevel = 0 : fatal error messages only
iolevel = 1 : error messages and minimal information
iolevel = 2 : error messages and detailed information
On exit, iparam(6) is unchanged.
iparam(7): idefault
On entry, defines if the default values should be used
in arrays IPARAM and RPARAM. If idefault = 0, then the
following default values are assigned:
IPARAM(1) = niparam = 100
IPARAM(6) = iolevel = 0
IPARAM(8) = istore = 1
IPARAM(9) = itrans = 0
IPARAM(10) = itmax = 5
If idefault = 1, then you must assign values to the
above variables before the call to the DUSKYR routine.
On exit, iparam(7) is unchanged.
iparam(8): istore
On entry, defines the type of storage scheme used for
the skyline matrix. If istore = 1, the unsymmetric
matrix A is stored using the profile-in storage mode;
if istore = 2, the unsymmetric matrix A is stored
using the diagonal-out storage mode; if istore = 3,
the unsymmetric matrix A is stored using the
structurally symmetric profile-in storage mode.
Default: istore = 1.
On exit, iparam(8) is unchanged.
iparam(9): itrans
On entry, defines the form of matrix used in the
iterative refinement. If itrans = 0, the system
refined is A X = B; if itrans = 1, the system refined
is trans(A)*X = B. Default: itrans = 0.
On exit, iparam(9) is unchanged.
iparam(10): itmax
On entry, defines the maximum number of iterations for
the iterative refinement process. Default: itmax = 5.
On exit, iparam(10) is unchanged.
rparam real*8
An array of length at least 100, containing the real
parameters for the iterative refinement and error
bounds calculation.
On exit, rparam is unchanged. rparam is not used by
the routine DUSKYR at present, but is reserved for
future use. It can be a dummy variable.
iwrk integer*4
On entry, an array of length at least 5n used for
integer workspace. The first 4n elements of the array
IWRK, generated by the routine DUSKYF, should be passed
unchanged to the routine DUSKYR.
On exit, the first 4n elements of iwrk are unchanged.
rwrk real*8
On entry, an array of length at least 3n used for real
workspace.
On exit, the first 3n elements of rwrk are
overwritten.
ierror integer*4
On entry, an unspecified variable.
On exit, ierror contains the error flag. A value of
zero indicates a normal exit from the routine DUSKYR.
Description
DUSKYR obtains an improved solution to the system
AX = B
or
trans(A) X = B
via iterative refinement. This is done by calculating the matrix of
residuals R using the matrix of solutions X_hat obtained from DUSKYS, and
obtaining a new matrix of solutions X(new) as follows:
For itrans = 0:
R = B - A * X_hat
delta_X = inverse(A) * R
and
X(new) = X_hat + delta_X
For itrans = 1:
R = B - trans(A) * X_hat
delta_X = inv_transp(A) * R
and
X(new) = X_hat + delta_X
The process of iterative refinement therefore requires both the original
matrix A as well as the L*D*U factorization obtained via the routine
DUSKYF. Since this routine overwrites the matrix A by the factorization, a
copy of the matrix must be made before the call to DUSKYF. Further, both
the right sides B and the solution vectors X_hat are required during
iterative refinement. Since the solution process in the routine DUSKYS
overwrites the right sides with the solution vectors, a copy of the right
sides must be made before the call to the routine DUSKYS.
In addition to the iterative refinement of the solution vectors, the
routine DUSKYR also provides the component-wise relative backward error,
berr and the estimated forward error bound, ferr, for each solution vector
[Arioli, Demmel, Duff 1989, Anderson et. al. 1992]. berr is the smallest
relative change in any entry of A or B that makes x_hat an exact solution.
ferr bounds the magnitude of the largest entry in x_hat - x(true) divided
by the magnitude of the largest entry in x_hat.
The process of iterative refinement is continued as long as all of the
following conditions are satisfied [Arioli, Demmel, Duff 1989]:
• The number of iterations of the iterative refinement process is less
than IPARAM(10) = itmax.
• berr is reduced by at least a factor of 2 during the previous iteration.
• berr is larger than the machine precision.
The routine DUSKYR is called after a call to the routine DSSKYF to obtain
the L*D*U factorization and a call to the routine DUSKYS to obtain the
solution x_hat. The first 4n elements of the integer workspace array IWRK,
generated by DUSKYF, contain information for use by DUSKYR and therefore
must remain unchanged between the calls to the routines DUSKYF and DUSKYR.
The storage scheme used in the routines DUSKYF, DUSKYS, and DUSKYR must be
identical. The value of itrans must be the same in the calls to the
routines DUSKYS and DUSKYR.
CXML Home Page Index of CXML Routines