Digital

Name

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

ninteger∗4
On entry, the order of the matrix A.
On exit, n is unchanged. 

aureal∗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. 

aufreal∗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. 

iaudiaginteger∗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. 

nauinteger∗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. 

alreal∗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. 

alfreal∗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. 

ialdiaginteger∗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. 

nalinteger∗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. 

breal∗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. 

ldbinteger∗4
On entry, the leading dimension of  array B.
On exit, ldb is unchanged. 

xreal∗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. 

ldxinteger∗4
On entry, the leading dimension of  array X. ldx >=n.
On exit, ldx is unchanged. 

ferrreal∗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. 

berrreal∗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. 

nbxinteger∗4
On entry, the number of right sides.
On exit, nbx is unchanged. 

iparaminteger∗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. 

rparamreal∗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. 

iwrkinteger∗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. 

ierrorinteger∗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.