Digital

Name

dsskyr − Symmetric sparse iterative refinement using skyline storage scheme

FORMAT

DSSKYR (n, au, auf, iaudiag, nau, 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 of length at least nau, containing the matrix A stored in the skyline storage scheme, using either the profile-in or the diagonal-out storage mode.
On exit,  au is unchanged. 

aufreal∗8
On entry, an array of length at least nau, containing transp(U)∗D∗U factorization of the the matrix A stored in the skyline storage scheme, using either the profile-in or the diagonal-out storage mode. The factorization has been obtained by a prior call to the routine DSSKYF.  auf must remain unchanged between calls to the routines DSSKYF and DSSKYR. 
On exit, auf is unchanged. 

iaudiaginteger∗4
On entry, an array of length at least n for the profile-in storage mode and (n+1) for the diagonal-out storage mode, containing the pointers to the locations of the diagonal elements in arrays AU and AUF.
On exit,  iaudiag is unchanged. 

nauinteger∗4
On entry, the number of elements in array AU. nau is also the envelope size of the symmetric part of the matrix A. For the profile-in storage mode, nau =  IAUDIAG(n).  For the diagonal-out storage mode, nau =  IAUDIAG(n+1) - 1. 
On exit,  nau 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. ldb >=n.
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 DSSKYS.
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 nbx solution vectors. 

berr
On entry, an array BERR of length at least nbz, whose elements are unspecified variables.
On exit, berr contains the component-wise relative backward error for each of the nbz solution vectors. 

nbx
On entry, the number of right sides.
On exit, nbz 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 parameter array is not used at present, nrparam can be unspecified. 
On exit,  iparam(2) is unchanged. 

iparam(3): niwrk
On entry, defines the size of the integer work array, IWRK.  niwrk >= 3n.
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 DSSKYR. 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) = itmax = 5

If idefault = 1, then you must assign values to the above variables before the call to the DSSKYR 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 matrix A is stored using the profile-in storage mode; if istore = 2, the matrix A is stored using the diagonal-out storage mode.  Default: istore = 1. 
On exit,  iparam(8) is unchanged. 

iparam(9): itmax
On entry, defines the maximum number of iterations for the iterative refinement process. Default: itmax = 5. 
On exit,  iparam(9) 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 DSSKYR at present, but is reserved for future use.  It can be a dummy variable. 

iwrkinteger∗4
On entry, an array of length at least 3n used for integer workspace. The first 2n elements of  the array IWRK, generated by the routine DSSKYF, should be passed unchanged to the routine DSSKYR.
On exit, the first 2n 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 DSSKYR. 

Description

DSSKYR obtains an improved solution to the system

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 DSSKYS, and obtaining a new matrix of solutions X(new) as follows:

R = B - A ∗ X_hat

delta_X = inverse(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 transp(U)∗D∗U factorization obtained via the routine DSSKYF. Since this routine overwrites the matrix A by the factorization, a copy of the matrix must be made prior to the call to DSSKYF. Further, both the right sides B and the solution vectors X are required during iterative refinement. Since the solution process in the routine DSSKYS overwrites the right sides with the solution vectors, a copy of the right sides must be made prior to the call to the routine DSSKYS. 

In addition to the iterative refinement of the solution vectors, the routine DSSKYR 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(9) = itmax. 

•berr reduces by at least a factor of 2 during the previous iteration. 

•berr is larger than the machine precision. 

The routine DSSKYR is called after a call to the routine DSSKYF to obtain the transp(U)∗D∗U factorization and a call to the routine DSSKYS to obtain the solution X.  The first 2n elements of the integer workspace array IWRK, generated by DSSKYF, contain information for use by DSSKYR and therefore must remain unchanged between the calls to the routines DSSKYF and DSSKYR. The real work array, RWRK, is not used at present.  The storage scheme used in the routines DSSKYF, DSSKYS, and DSSKYR must be identical.