CXML
duskyx
Unsymmetric sparse expert driver using skyline storage scheme
FORMAT
DUSKYX
(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 RPARAM(9) = ifactor = 0, auf is an
unspecified array of length at least nau. If ifactor
= 1, then auf contains information on the L*D*U
factorization of the matrix A. 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.
On exit, if ifactor = 0, then auf contains
information on the L*D*U factorization of the matrix A.
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. If ifactor = 1, then 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 (if ifactor = 0). iaudiag is of length at least n
for the profile-in and the structurally symmetric
profile-in storage modes. 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 IPARAM(9) = ifactor = 0, alf is an
unspecified array of length at least nal. If ifactor =
1, then alf contains information on the L*D*U
factorization of the matrix A. 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.
On exit, if ifactor = 0, alf contains information on
the L*D*U factorization of the matrix A. 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. If ifactor =
1, then 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 (if ifactor = 1). 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. ldb >=n.
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 nbx 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 nbx solution vectors.
nbx integer*4
On entry, the number of right sides.
On exit, nbz is unchanged.
iparam integer*4
An array of length at least 100, containing the integer
parameters for the expert driver.
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.
nrparam >= 100.
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 DUSKYX.
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(2) = nrparam = 100
IPARAM(6) = iolevel = 0
IPARAM(8) = istore = 1
IPARAM(9) = ifactor = 0
IPARAM(10) = idet = 0
IPARAM(11) = ipvt = 0
IPARAM(13) = itrans = 0
IPARAM(14) = itmax = 5
RPARAM(1) = pvt_sml = 10**(-12)
If idefault = 1, then you must assign values to the
above variables before the call to the DUSKYX 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): ifactor
On entry, defines if the matrix A has already been
factored. If ifactor = 0, the matrix is unfactored
and arrays AUF and ALF are unspecified. If ifactor =
1, the matrix has been factored by a prior call to the
routine DUSKYF, and the arrays AUF and array ALF
contain the L*D*U factorization of A. Default:
ifactor = 0.
On exit, iparam(9) is unchanged.
iparam(10): idet
On entry, defines if the determinant of the matrix A is
to be calculated. If idet = 0, then the determinant
is not calculated; if idet = 1, the determinant is
calculated as det_base * 10**det_pwr. See RPARAM(4)
and RPARAM(5). If ifactor = 1, then IPARAM(10) is
unspecified. Default: idet = 0.
On exit, iparam(10) is unchanged.
iparam(11): ipvt
On entry, defines if the factorization should continue
when a small pivot, defined by RPARAM(1), is
encountered. If ipvt = 0 and the absolute value of the
pivot element is smaller than pvt_sml = RPARAM(1),
then the factorization process is stopped and control
returned to the calling subprogram. If ipvt = 1 and a
pivot smaller than RPARAM(1) in absolute value is
encountered in the factorization, the process
continues. If ipvt = 2 and a pivot smaller than
RPARAM(1) in absolute value is encountered in the
factorization, it is replaced by a predetermined value
pvt_new = RPARAM(2), and the factorization is
continued. If ifactor = 1, then IPARAM(12) is
unspecified. Default: ipvt = 0.
On exit, iparam(11) is unchanged.
iparam(12): ipvt_loc
On entry, an unspecified variable.
On exit, iparam(12) contains the location of the first
pivot element smaller in absolute value than pvt_sml.
The pivot element is returned in pvt_val = RPARAM(3).
If iparam(12) = 0, then no such pivot element exists.
iparam(13): 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(13) is unchanged.
iparam(14): itmax
On entry, defines the maximum number of iterations for
the iterative refinement process. Default: itmax = 5.
On exit, iparam(14) is unchanged.
rparam real*8
An array of length at least 100, containing the real
parameters for the expert driver.
rparam(1): pvt_sml
On entry, defines the value of the pivot element which
is considered to be small. If a pivot element smaller
than pvt_sml, in absolute value, is encountered in the
factorization process, then, depending on the value of
ipvt = IPARAM(11), the process either stops,
continues or continues after the pivot is set equal to
pvt_new = RPARAM(2). If ifactor = 1, then RPARAM (1)
is unspecified. pvt_sml > 0. Recommended value:
10**(-15) <= pvt_sml <= 1. Default: pvt_sml = 10**(-
12).
On exit, rparam(1) is unchanged.
rparam(2): pvt_new
On entry, defines the value to which the pivot element
must be set if ipvt = 2 and the pivot element is less
than pvt_sml in absolute value. pvt_sml should be
large enough to avoid overflow when calculating the
reciprocal of the pivot element. If ifactor = 1, the
RPARAM(2) is unspecified.
On exit, rparam(2) is unchanged.
rparam(3): pvt_val
On entry, an unspecified variable.
On exit, rparam(3) contains the value of the first
pivot element smaller than pvt_sml in absolute value.
The location of this element is returned in ipvt_loc =
IPARAM(12). If ifactor = 1, then the RPARAM(3) is
unspecified.
rparam(4): det_base
On entry, an unspecified variable.
On exit, defines the base for the determinant of the
matrix A. If idet = 1, the determinant is calculated
as det_base * 10**(det_pwr). If ifactor = 1, then
RPARAM(4) is unspecified. 1.0 <= det_base <= 10.0.
rparam(5): det_pwr
On entry, an unspecified variable.
On exit, defines the power for the determinant of the
matrix A. If idet = 1, the determinant is calculated
as det_base * 10**det_pwr. If ifactor = 1, then
RPARAM(5) is unspecified.
rparam(6): anorm
On entry, an unspecified variable.
On exit, rparam(6) contains the 1-norm or the
Infinity-norm of the matrix A.
rparam(7): ainorm
On entry, an unspecified variable.
On exit, rparam(7) contains the estimate of the 1-norm
or the Infinity-norm of inverse(A)).
rparam(8): rcond
On entry, an unspecified variable.
On exit, rparam(8) contains the reciprocal of the
estimate of the 1-norm or the Infinity-norm condition
number of the matrix A.
iwrk integer*4
On entry, an array of length at least 5n used for
integer workspace. If ifactor = 1, then the first 4n
elements of the array IWRK contain information
generated by the routine DUSKYF. If ifactor = 0, then
this information is unspecified.
On exit, the first 4n elements of the array IWRK
contain information generated by the routine DUSKYF.
This information is used by the routines DUSKYS and
DUSKYR, and should therefore remain unchanged between
the call to the routine DUSKYX and any subsequent call
to the routines DUSKYS and DUSKYR.
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 DUSKYX.
Description
DUSKYX is an expert driver routine that :
• Obtains the L*D*U factorization of the matrix A via a call to the
routine DUSKYF.
• If the factorization is successful, obtains the 1-norm or Infinity-norm
condition number estimate of the matrix A by a call to the routine
DUSKYC.
• If the reciprocal of the condition number estimate is greater than the
machine precision, DUSKYX uses the factorization to solve the system
A X = B
or
trans(A)X = B
using the routine DUSKYS.
• Improves the solution X via iterative refinement and obtains the error
bounds using the routine DUSKYR.
DUSKYX first obtains the factorization of the symmetric matrix A as:
A = L D U
where L is a unit triangular matrix, D is a diagonal matrix and U is a unit
upper triangular matrix. The matrix A is stored in a skyline form, using
either the profile-in storage mode or the diagonal-out storage mode or the
structurally symmetric profile-in storage mode. If the matrix is already
factored, as indicated by ifactor, then this step is skipped.
The routine DUSKYF does not perform any pivoting to preserve the numerical
stability of the L*D*U factorization. It is therefore primarily intended
for the solution of systems that do not require pivoting for numerical
stability, such as diagonally dominant systems. Caution is urged when
using this routine for problems that require pivoting.
If a small pivot, in absolute value, pvt_sml, is encountered in the
process of factorization, you have the option of either stopping the
factorization process and returning to the calling subprogram, continuing
the factorization process with the small value of the pivot, or continuing
after setting the pivot equal to some predetermined value, pvt_new. The
location of the first occurrence of a small pivot is returned in ipvt_loc
and its value in pvt_val.
In addition to the L*D*U factorization, the routine DUSKYF can be used to
obtain the determinant of A. If factorization process is stopped at row i
due to a small pivot, then the determinant are evaluated for rows 1 through
(i-1).
The routine DUSKYX does not allow a partial factorization of the matrix A.
If a partial factorization of A is required, the routine DUSKYF is
recommended.
DUSKYC obtains the reciprocal of the estimate of the condition number of
the unsymmetric matrix A as:
rcond(A) = 1 / (||A||*||inverse(A)||)
If the system being solved is
A * X = B
the reciprocal of the 1-norm condition number estimate is calculated. If
the system being solved is
trans(A) * X = B
the reciprocal if the Infinity-norm of condition number estimate is
calculated.
The 1-norm of inverse(A) or inv_transp(A) is obtained using the LAPACK
routine DLACON, which uses Higham's modification [Higham 1988] of Hager's
method [Hager 1984]. If the reciprocal of the condition number estimate is
larger than the machine precision, the routine DUSKYX solves the system via
a call to the routine DUSKYS and then improves on the solution 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
In addition to the iterative refinement of the solution vectors, the
routine DUSKYX 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 reduces by at least a factor of 2 during the previous iteration.
• berr is larger than the machine precision.
The first 4n elements of the integer workspace array IWRK generated by
DUSKYF, contain information for use by the routines DUSKYS and DUSKYR and
therefore must remain unchanged between the calls to the routine DUSKYX and
any subsequent calls to the routine DUSKYS and DUSKYR.
CXML Home Page Index of CXML Routines