CXML
duskyf
Unsymmetric sparse matrix factorization using skyline storage
scheme (Serial and Parallel Versions)
FORMAT
DUSKYF (n, au, iaudiag, nau, al, ialdiag, nal, 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, if istore = 1 or 2, au contains the factors
U and D of the L*D*U factorization of the matrix A. If
istore = 3, then au contains the factors L, U and D of
the L*D*U factorization of the matrix A. au must
remain unchanged between the call to the routine DUSKYF
and any routines that use the factors such as DUSKYS,
DUSKYC, DUSKYR, and DUSKYX.
iaudiag integer*4
On entry, an array containing the pointers to the
locations of the diagonal elements in array AU.
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, if istore = 1 or 2, al contains the factor L
of the L*D*U factorization of the matrix A. If
istore = 3, then al is undefined. al must remain
unchanged between the call to the routine DUSKYF and
any routines that use the factors such as DUSKYS,
DUSKYC, DUSKYR, and DUSKYX.
ialdiag integer*4
On entry, an array containing the pointers to the
locations of the diagonal elements in array AL.
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.
iparam integer*4
An array of length at least 100, containing the integer
parameters for the L*D*U factorization.
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 IPARAM.
nrparam >= 100.
On exit, iparam(2) is unchanged.
iparam(3): niwrk
On entry, defines the size of the integer work array,
IWRK. niwrk >= 4n.
On exit, iparam(3) is unchanged.
iparam(4): nrwrk
On entry, defines the size of the real work array,
RWRK. As the real work array is not used at present,
nrwrk may be unspecified.
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 DUSKYF.
The I/O unit must be opened in the calling subprogram.
If iounit <= 0, then 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 : fatal errors, warnings and minimal
information
iolevel = 2 : detailed information and statistics
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) = ibeg = 0
IPARAM(10) = idet = 0
IPARAM(11) = ipvt = 0
RPARAM(1) = pvt_sml = 10**(-12)
If idefault = 1, then you must assign values to the
above variables before the call to the DUSKYF 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): ibeg
On entry, defines if full or partial factorization is
to be performed. If ibeg = 0, then a full
factorization is performed for rows and columns 1
through n. If ibeg > 0, then a partial factorization
is performed starting from rows and columns ibeg+1
through n, that is, rows and columns from 1 through
ibeg have already been factorized. Default: ibeg = 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.
rparam real*8
An array of length at least 100, containing the real
parameters for the L*D*U factorization.
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. pvt_new must be large enough to avoid
overflow when calculating the reciprocal of the pivot
element.
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).
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). 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).
iwrk integer*4
On entry, an array of length at least 4n used for
integer workspace.
On exit, iwrk contains information for use by routines
that use the factorization such as DUSKYS, DUSKYC,
DUSKYR and DUSKYX. The first 4n elements of iwrk
should therefore be passed unchanged to these routines.
rwrk real*8
On entry, an array used for real workspace.
On exit, rwrk is unchanged. rwrk is not used by the
routine DUSKYF at present but is reserved for future
use. It can be a dummy variable.
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 DUSKYF.
Description
DUSKYF obtains the factorization of the unsymmetric matrix A as:
A = L D U
where L is a unit lower 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, the diagonal-out storage mode, or
the structurally symmetric profile-in storage mode.
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 program, continuing the factorization
process, 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. A partial factorization can also be obtained
by appropriately setting the value of ibeg. If ibeg > 0, then
factorization begins at row and column ibeg + 1; the rows and columns from
1 to ibeg are assumed to have been already factorized. When ibeg > 0, the
determinant of A and the statistics on the matrix are calculated from rows
and columns ibeg + 1 through n. If the factorization process is stopped at
row i due to a small pivot, then the determinant and the statistics on the
matrix are evaluated for rows ibeg + 1 through (i-1).
The data in the first 4n elements of the integer workspace array, IWRK, are
used in routines that use the L*D*U factorization, such as DUSKYS, DUSKYC,
and DUSKYR. This data must therefore remain unchanged between the call to
DUSKYF and any one of these routines. The real workspace array, RWRK, is
not used at present.
This routine is available in both serial and parallel versions. The routine
names and parameter list are identical for both versions. For information
about linking to the serial or to the parallel library, refer to the CXML
Reference Manual.
CXML Home Page Index of CXML Routines