DGBSVX  - use the LU factorization to compute the solution
       to a real system of linear equations A * X = B, A**T * X =
       B, or A**H * X = B,


SYNOPSIS

       SUBROUTINE DGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB,
                          AFB, LDAFB, IPIV, EQUED, R, C, B,  LDB,
                          X, LDX, RCOND, FERR, BERR, WORK, IWORK,
                          INFO )

           CHARACTER      EQUED, FACT, TRANS

           INTEGER        INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N,
                          NRHS

           DOUBLE         PRECISION RCOND

           INTEGER        IPIV( * ), IWORK( * )

           DOUBLE         PRECISION  AB( LDAB, * ), AFB( LDAFB, *
                          ), B( LDB, * ), BERR(  *  ),  C(  *  ),
                          FERR( * ), R( * ), WORK( * ), X( LDX, *
                          )


PURPOSE

       DGBSVX uses the LU factorization to compute  the  solution
       to a real system of linear equations A * X = B, A**T * X =
       B, or A**H * X = B, where A is a band matrix  of  order  N
       with  KL  subdiagonals  and KU superdiagonals, and X and B
       are N-by-NRHS matrices.

       Error bounds on the solution and a condition estimate  are
       also provided.


DESCRIPTION

       The following steps are performed by this subroutine:

       1.  If  FACT  =  'E', real scaling factors are computed to
       equilibrate
          the system:
             TRANS = 'N':  diag(R)*A*diag(C)      *inv(diag(C))*X
       = diag(R)*B
             TRANS  = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X
       = diag(C)*B
             TRANS = 'C': (diag(R)*A*diag(C))**H  *inv(diag(R))*X
       = diag(C)*B
          Whether  or not the system will be equilibrated depends
       on the
          scaling of the matrix A, but if equilibration is  used,
       A is
          overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if
       TRANS='N')


       2. If FACT = 'N' or 'E', the LU decomposition is  used  to
       factor the
          matrix A (after equilibration if FACT = 'E') as
             A = L * U,
          where L is a product of permutation and unit lower tri­
       angular
          matrices with KL subdiagonals, and U is upper  triangu­
       lar with
          KL+KU superdiagonals.

       3.  If  some U(i,i)=0, so that U is exactly singular, then
       the routine
          returns with INFO = i. Otherwise, the factored form  of
       A is used
          to  estimate  the condition number of the matrix A.  If
       the
          reciprocal of the condition number is less than machine
       precision,
          INFO  =  N+1  is returned as a warning, but the routine
       still goes on
          to solve for X and compute error  bounds  as  described
       below.

       4.  The system of equations is solved for X using the fac­
       tored form
          of A.

       5. Iterative refinement is applied to improve the computed
       solution
          matrix  and  calculate  error bounds and backward error
       estimates
          for it.

       6. If equilibration was used, the matrix  X  is  premulti­
       plied by
          diag(C)  (if TRANS = 'N') or diag(R) (if TRANS = 'T' or
       'C') so
          that it solves the original  system  before  equilibra­
       tion.


ARGUMENTS

       FACT    (input) CHARACTER*1
               Specifies  whether or not the factored form of the
               matrix A is supplied on entry, and if not, whether
               the  matrix  A should be equilibrated before it is
               factored.  = 'F':  On entry, AFB and IPIV  contain
               the  factored form of A.  If EQUED is not 'N', the
               matrix A has been equilibrated with  scaling  fac­
               tors  given by R and C.  AB, AFB, and IPIV are not
               modified.  = 'N':  The matrix A will be copied  to
               AFB and factored.

               essary, then copied to AFB and factored.

       TRANS   (input) CHARACTER*1
               Specifies the form of the system of equations.   =
               'N':  A * X = B     (No transpose)
               = 'T':  A**T * X = B  (Transpose)
               = 'C':  A**H * X = B  (Transpose)

       N       (input) INTEGER
               The number of linear equations, i.e., the order of
               the matrix A.  N >= 0.

       KL      (input) INTEGER
               The number of subdiagonals within the band  of  A.
               KL >= 0.

       KU      (input) INTEGER
               The number of superdiagonals within the band of A.
               KU >= 0.

       NRHS    (input) INTEGER
               The number of right hand sides, i.e.,  the  number
               of columns of the matrices B and X.  NRHS >= 0.

       AB      (input/output) DOUBLE PRECISION array, dimension
               (LDAB,N)
               On  entry, the matrix A in band storage, in rows 1
               to KL+KU+1.  The j-th column of A is stored in the
               j-th column of the array AB as follows: AB(KU+1+i-
               j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)

               If FACT = 'F' and EQUED is not 'N',  then  A  must
               have been equilibrated by the scaling factors in R
               and/or C.  AB is not modified if  FACT  =  'F'  or
               'N', or if FACT = 'E' and EQUED = 'N' on exit.

               On  exit,  if  EQUED .ne. 'N', A is scaled as fol­
               lows: EQUED = 'R':  A := diag(R) * A
               EQUED = 'C':  A := A * diag(C)
               EQUED = 'B':  A := diag(R) * A * diag(C).

       LDAB    (input) INTEGER
               The leading dimension of the array  AB.   LDAB  >=
               KL+KU+1.

       AFB      (input  or  output)  DOUBLE  PRECISION  array,
               dimension (LDAFB,N)
               If  FACT  = 'F', then AFB is an input argument and
               on entry contains details of the LU  factorization
               of the band matrix A, as computed by DGBTRF.  U is
               stored as an upper  triangular  band  matrix  with
               KL+KU superdiagonals in rows 1 to KL+KU+1, and the
               multipliers  used  during  the  factorization  are

               .ne. 'N', then AFB is the  factored  form  of  the
               equilibrated matrix A.

               If  FACT = 'N', then AFB is an output argument and
               on exit returns details of the LU factorization of
               A.

               If  FACT = 'E', then AFB is an output argument and
               on exit returns details of the LU factorization of
               the  equilibrated matrix A (see the description of
               AB for the form of the equilibrated matrix).

       LDAFB   (input) INTEGER
               The leading dimension of the array AFB.  LDAFB  >=
               2*KL+KU+1.

       IPIV    (input or output) INTEGER array, dimension (N)
               If  FACT = 'F', then IPIV is an input argument and
               on entry contains the pivot indices from the  fac­
               torization A = L*U as computed by DGBTRF; row i of
               the matrix was interchanged with row IPIV(i).

               If FACT = 'N', then IPIV is an output argument and
               on  exit  contains the pivot indices from the fac­
               torization A = L*U of the original matrix A.

               If FACT = 'E', then IPIV is an output argument and
               on  exit  contains the pivot indices from the fac­
               torization A = L*U of the equilibrated matrix A.

       EQUED   (input or output) CHARACTER*1
               Specifies the form of equilibration that was done.
               =  'N':   No  equilibration (always true if FACT =
               'N').
               = 'R':  Row equilibration, i.e., A has  been  pre­
               multiplied  by diag(R).  = 'C':  Column equilibra­
               tion, i.e., A has been postmultiplied by  diag(C).
               = 'B':  Both row and column equilibration, i.e., A
               has been replaced by diag(R) * A * diag(C).  EQUED
               is  an input argument if FACT = 'F'; otherwise, it
               is an output argument.

       R       (input or output) DOUBLE PRECISION array,
               dimension (N)
               The  row  scale  factors for A.  If EQUED = 'R' or
               'B', A is multiplied on the left  by  diag(R);  if
               EQUED  =  'N'  or 'C', R is not accessed.  R is an
               input argument if FACT = 'F'; otherwise, R  is  an
               output argument.  If FACT = 'F' and EQUED = 'R' or
               'B', each element of R must be positive.

       C       (input or output) DOUBLE PRECISION array,
               dimension (N)

               'B', A is multiplied on the right by  diag(C);  if
               EQUED  =  'N'  or 'R', C is not accessed.  C is an
               input argument if FACT = 'F'; otherwise, C  is  an
               output argument.  If FACT = 'F' and EQUED = 'C' or
               'B', each element of C must be positive.

       B        (input/output)   DOUBLE   PRECISION   array,
               dimension (LDB,NRHS)
               On entry, the right hand side matrix B.  On  exit,
               if  EQUED = 'N', B is not modified; if TRANS = 'N'
               and EQUED =  'R'  or  'B',  B  is  overwritten  by
               diag(R)*B;  if  TRANS = 'T' or 'C' and EQUED = 'C'
               or 'B', B is overwritten by diag(C)*B.

       LDB     (input) INTEGER
               The leading dimension of  the  array  B.   LDB  >=
               max(1,N).

       X       (output) DOUBLE PRECISION array, dimension
               (LDX,NRHS)
               If INFO = 0 or INFO = N+1, the N-by-NRHS  solution
               matrix  X  to  the  original  system of equations.
               Note that A and B are modified on  exit  if  EQUED
               .ne.  'N',  and  the  solution to the equilibrated
               system is inv(diag(C))*X if TRANS = 'N' and  EQUED
               =  'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or
               'C' and EQUED = 'R' or 'B'.

       LDX     (input) INTEGER
               The leading dimension of  the  array  X.   LDX  >=
               max(1,N).

       RCOND   (output) DOUBLE PRECISION
               The estimate of the reciprocal condition number of
               the matrix A after equilibration  (if  done).   If
               RCOND  is less than the machine precision (in par­
               ticular, if RCOND = 0), the matrix is singular  to
               working precision.  This condition is indicated by
               a return code of INFO > 0.

       FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
               The estimated forward error bound for  each  solu­
               tion  vector X(j) (the j-th column of the solution
               matrix X).  If XTRUE is the true  solution  corre­
               sponding  to  X(j),  FERR(j) is an estimated upper
               bound for the magnitude of the largest element  in
               (X(j)  -  XTRUE)  divided  by the magnitude of the
               largest element in X(j).  The estimate is as reli­
               able  as  the  estimate  for  RCOND, and is almost
               always a slight overestimate of the true error.

       BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
               The componentwise relative backward error of  each

               change in any element of A or B that makes X(j) an
               exact solution).

       WORK     (workspace/output)  DOUBLE  PRECISION  array,
               dimension (3*N)
               On  exit,  WORK(1)  contains  the reciprocal pivot
               growth factor norm(A)/norm(U). The  "max  absolute
               element"  norm  is  used.  If WORK(1) is much less
               than 1, then the stability of the LU factorization
               of the (equilibrated) matrix A could be poor. This
               also means that the solution X, condition  estima­
               tor  RCOND,  and forward error bound FERR could be
               unreliable. If factorization fails with 0<INFO<=N,
               then  WORK(1) contains the reciprocal pivot growth
               factor for the leading INFO columns of A.

       IWORK   (workspace) INTEGER array, dimension (N)

       INFO    (output) INTEGER
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an ille­
               gal value
               > 0:  if INFO = i, and i is
               <=  N:  U(i,i) is exactly zero.  The factorization
               has been completed, but the factor  U  is  exactly
               singular,  so  the solution and error bounds could
               not be computed. RCOND = 0 is returned.  = N+1:  U
               is  nonsingular,  but  RCOND  is less than machine
               precision, meaning that the matrix is singular  to
               working precision.  Nevertheless, the solution and
               error bounds are computed because there are a num­
               ber  of situations where the computed solution can
               be more accurate than the  value  of  RCOND  would
               suggest.


Man(1) output converted with man2html