CXML
dsskyx
Symmetric sparse expert driver using skyline storage scheme
FORMAT
DSSKYX
(n, au, auf, iaudiag, nau, 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 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.
auf real*8
On entry, if RPARAM(9) = ifactor = 0, auf is an
unspecified array of length at least nau. If ifactor
= 1, auf is an array of length at least nau, containing
the transp(U)*D*U factorization of 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.
On exit, if ifactor = 0, auf contains the
transp(U)*D*U factorization of the matrix A stored in
the skyline storage scheme, using either the profile-in
or the diagonal-out storage mode. If ifactor = 1, then
auf is unchanged.
iaudiag integer*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 (if
ifactor =1).
On exit, iaudiag is unchanged.
nau integer*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.
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.
On exit, x contains the 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
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, nbx 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 >=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 DSSKYX.
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) = inertia = 0
IPARAM(17) = 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 DSSKYX 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): ifactor
On entry, defines if the matrix A has already been
factored on input to the routine DSSKYX. If ifactor =
0, the matrix is unfactored and array AUF is
unspecified. If ifactor = 1, the matrix has been
factored by a prior call to the routine DSSKYF, and
array AUF contains the transp(U)*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). 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. 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.
If ifactor = 1, then IPARAM (12) is unspecified.
iparam(13): inertia
On entry, defines if the inertia of the matrix A should
be calculated during factorization. The inertia of the
symmetric matrix A is the triplet of integers (ipeigen,
ineigen, izeigen), consisting of the number of
positive, negative and zero eigenvalues, respectively.
If inertia = 0, then the inertia is not calculated; if
inertia = 1, then the number of positive and negative
eigenvalues are returned in ipeigen = IPARAM(14) and
ineigen = IPARAM(15), respectively. An indication of
the existence of zero eigenvalues is returned in
izeigen = IPARAM(16). Default: inertia = 0.
On exit, iparam(13) is unchanged.
iparam(14): ipeigen
On entry, an unspecified variable.
On exit, if inertia = 1, iparam(14) contains the
number of positive eigenvalues of the matrix A. If
ifactor = 0, then IPARAM(14) is unspecified.
iparam(15): ineigen
On entry, an unspecified variable.
On exit, if inertia = 1, iparam(15) contains the
number of negative eigenvalues of the matrix A. If
ifactor = 0, then IPARAM(15) is unspecified.
iparam(16): izeigen
On entry, an unspecified variable.
On exit, if inertia = 1, iparam(16) indicates if the
matrix A has any zero eigenvalues. If izeigen = 0,
then the matrix A does not have a zero eigenvalue; if
izeigen = 1, then the matrix A has at least one zero
eigenvalue. If ifactor = 0, the IPARAM(16) is
unspecified.
iparam(17): itmax
On entry, defines the maximum number of iterations for
the iterative refinement process. Default: itmax = 5.
On exit, iparam(17) 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). 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.
This element occurs at the location returned in
IPARAM(12). If no such pivot element is found, the
value of pvt_val is unspecified. If ifactor = 1, then
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, the
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, the
RPARAM(5) is unspecified.
rparam(6): anorm
On entry, an unspecified variable.
On exit, rparam(6) contains the 1-norm of the matrix A.
rparam(7): ainorm
On entry, an unspecified variable.
On exit, rparam(7) contains the estimate of the 1-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 condition number of the matrix
A.
iwrk integer*4
On entry, an array of length at least 3n used for
integer workspace. If ifactor = 1, then the first 2n
elements of the array IWRK contain information
generated by the routine DSSKYF. If ifactor = 0, then
this information is unspecified.
On exit, the first 2n elements of the array IWRK
contain information generated by the routine DSSKYF.
This information is used by the routines DSSKYS and
DSSKYR, and should therefore remain unchanged between
the call to the routine DSSKYX and any subsequent call
to the routines DSSKYS and DSSKYR.
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 DSSKYX.
Description
DSSKYX is an expert driver routine that:
• Obtains the L*D*U factorization of the matrix A via a call to the
routine DSSKYF.
• If the factorization is successful, obtains the 1-norm condition number
estimate of the matix A by a call to the routine DSSKYC.
• If the reciprocal of the condition number estimate is greater than the
machine precision, DSSKYX uses the factorization to solve the system
A X = B
using the routine DSSKYS.
• Improves the solution X via iterative refinement and obtains the error
bounds using the routine DSSKYR.
DSSKYX first obtains the factorization of the symmetric matrix A as:
A = transp(U)*D*U
where 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. If the matrix is already factored,
as indicated by ifactor, then this step is skipped.
The routine DSSKYF does not perform any pivoting to preserve the numerical
stability of the transp(U)*D*U factorization. It is therefore primarily
intended for the solution of symmetric positive (or negative) definite
systems as they do not require pivoting for numerical stability. Caution
is urged when using this routine for symmetric indefinite systems.
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 transp(U)*D*U factorization, you can also obtain the
determinant of A, the number of positive and negative eigenvalues of the
matrix A, and an indication of the existence of zero eigenvalues. If the
factorization process is stopped at row i due to a small pivot, then the
inertia and determinant are evaluated for rows 1 through (i-1).
The routine DSSKYX does not allow a partial factorization of the matrix A.
If a partial factorization of A is required, the routine DSSKYF is
recommended.
DUSKYC obtains the reciprocal of the estimate of the condition number of
the symmetric matrix A as:
rcond(A) = 1 / (||A|| * || inverse(A)||))
The 1-norm of inverse(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 DSSKYX solves the system via a call to the routine
DSSKYS 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 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
In addition to the iterative refinement of the solution vectors, the
routine DSSKYX 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
DSSKYF, contain information for use by the routines DSSKYS and DSSKYR.
They must therefore remain unchanged between the calls to the routine
DSSKYX and any subsequent calls to the routines DSSKYS and DSSKYR.
CXML Home Page Index of CXML Routines