CXML

duskyx 


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