pdgesvx.f

Go to the documentation of this file.

      SUBROUTINE pdgesvx( FACT, TRANS, N, NRHS, A, IA, JA, DESCA, AF,
     $                    IAF, JAF, DESCAF, IPIV, EQUED, R, C, B, IB,
     $                    JB, DESCB, X, IX, JX, DESCX, RCOND, FERR,
     $                    BERR, WORK, LWORK, IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      CHARACTER          EQUED, FACT, TRANS
      INTEGER            IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
     $                   LWORK, N, NRHS
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCAF( * ), DESCB( * ),
     $                   DESCX( * ), IPIV( * ), IWORK( * )
      DOUBLE PRECISION   A( * ), AF( * ), B( * ), BERR( * ), C( * ),
     $                   ferr( * ), r( * ), work( * ), x( * )
*     ..
*
*  Purpose
*  =======
*
*  PDGESVX uses the LU factorization to compute the solution to a real
*  system of linear equations
*
*        A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
*
*  where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
*  B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
*
*  Error bounds on the solution and a condition estimate are also
*  provided.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Description
*  ===========
*
*  In the following description, A denotes A(IA:IA+N-1,JA:JA+N-1),
*  B denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
*  X(IX:IX+N-1,JX:JX+NRHS-1).
*
*  The following steps are performed:
*
*  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')
*     or diag(C)*B (if TRANS = 'T' or 'C').
*
*  2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
*     matrix A (after equilibration if FACT = 'E') as
*        A = P * L * U,
*     where P is a permutation matrix, L is a unit lower triangular
*     matrix, and U is upper triangular.
*
*  3. 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, steps 4-6 are skipped.
*
*  4. The system of equations is solved for X using the factored 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 FACT = 'E' and equilibration was used, the matrix X is
*     premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if
*     TRANS = 'T' or 'C') so that it solves the original system
*     before equilibration.
*
*  Arguments
*  =========
*
*  FACT    (global input) CHARACTER
*          Specifies whether or not the factored form of the matrix
*          A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
*          whether the matrix A(IA:IA+N-1,JA:JA+N-1) should be
*          equilibrated before it is factored.
*          = 'F':  On entry, AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
*                  tain the factored form of A(IA:IA+N-1,JA:JA+N-1).
*                  If EQUED is not 'N', the matrix
*                  A(IA:IA+N-1,JA:JA+N-1) has been equilibrated with
*                  scaling factors given by R and C.
*                  A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),
*                  and IPIV are not modified.
*          = 'N':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be copied to
*                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
*          = 'E':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be equili-
*                  brated if necessary, then copied to
*                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
*
*  TRANS   (global input) CHARACTER
*          Specifies the form of the system of equations:
*          = 'N':  A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)     (No transpose)
*          = 'T':  A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
*          = 'C':  A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on, i.e. the
*          order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
*          N >= 0.
*
*  NRHS    (global input) INTEGER
*          The number of right-hand sides, i.e., the number of columns
*          of the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
*          X(IX:IX+N-1,JX:JX+NRHS-1).  NRHS >= 0.
*
*  A       (local input/local output) DOUBLE PRECISION pointer into
*          the local memory to an array of local dimension
*          (LLD_A,LOCc(JA+N-1)).  On entry, the N-by-N matrix
*          A(IA:IA+N-1,JA:JA+N-1).  If FACT = 'F' and EQUED is not 'N',
*          then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
*          the scaling factors in R and/or C.  A(IA:IA+N-1,JA:JA+N-1) is
*          not modified if FACT = 'F' or  'N', or if FACT = 'E' and
*          EQUED = 'N' on exit.
*
*          On exit, if EQUED .ne. 'N', A(IA:IA+N-1,JA:JA+N-1) is scaled
*          as follows:
*          EQUED = 'R':  A(IA:IA+N-1,JA:JA+N-1) :=
*                                      diag(R) * A(IA:IA+N-1,JA:JA+N-1)
*          EQUED = 'C':  A(IA:IA+N-1,JA:JA+N-1) :=
*                                      A(IA:IA+N-1,JA:JA+N-1) * diag(C)
*          EQUED = 'B':  A(IA:IA+N-1,JA:JA+N-1) :=
*                            diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  AF      (local input or local output) DOUBLE PRECISION pointer
*          into the local memory to an array of local dimension
*          (LLD_AF,LOCc(JA+N-1)).  If FACT = 'F', then
*          AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and on
*          entry contains the factors L and U from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by PDGETRF.
*          If EQUED .ne. 'N', then AF is the factored form of the
*          equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'N', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
*          argument and on exit returns the factors L and U from the
*          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
*          matrix A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'E', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
*          argument and on exit returns the factors L and U from the
*          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
*          brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
*          A(IA:IA+N-1,JA:JA+N-1) for the form of the equilibrated
*          matrix).
*
*  IAF     (global input) INTEGER
*          The row index in the global array AF indicating the first
*          row of sub( AF ).
*
*  JAF     (global input) INTEGER
*          The column index in the global array AF indicating the
*          first column of sub( AF ).
*
*  DESCAF  (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix AF.
*
*  IPIV    (local input or local output) INTEGER array, dimension
*          LOCr(M_A)+MB_A. If FACT = 'F', then IPIV is an input argu-
*          ment and on entry contains the pivot indices from the fac-
*          torization A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by
*          PDGETRF; IPIV(i) -> The global row local row i was
*          swapped with.  This array must be aligned with
*          A( IA:IA+N-1, * ).
*
*          If FACT = 'N', then IPIV is an output argument and on exit
*          contains the pivot indices from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
*          A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'E', then IPIV is an output argument and on exit
*          contains the pivot indices from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
*          A(IA:IA+N-1,JA:JA+N-1).
*
*  EQUED   (global input or global output) CHARACTER
*          Specifies the form of equilibration that was done.
*          = 'N':  No equilibration (always true if FACT = 'N').
*          = 'R':  Row equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1) has
*                  been premultiplied by diag(R).
*          = 'C':  Column equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1)
*                  has been postmultiplied by diag(C).
*          = 'B':  Both row and column equilibration, i.e.,
*                  A(IA:IA+N-1,JA:JA+N-1) has been replaced by
*                  diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
*          EQUED is an input variable if FACT = 'F'; otherwise, it is an
*          output variable.
*
*  R       (local input or local output) DOUBLE PRECISION array,
*                                                  dimension LOCr(M_A).
*          The row scale factors for A(IA:IA+N-1,JA:JA+N-1).
*          If EQUED = 'R' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
*          on the left by diag(R); if EQUED='N' or 'C', R is not acces-
*          sed.  R is an input variable if FACT = 'F'; otherwise, R is
*          an output variable.
*          If FACT = 'F' and EQUED = 'R' or 'B', each element of R must
*          be positive.
*          R is replicated in every process column, and is aligned
*          with the distributed matrix A.
*
*  C       (local input or local output) DOUBLE PRECISION array,
*                                                  dimension LOCc(N_A).
*          The column scale factors for A(IA:IA+N-1,JA:JA+N-1).
*          If EQUED = 'C' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
*          on the right by diag(C); if EQUED = 'N' or 'R', C is not
*          accessed.  C is an input variable if FACT = 'F'; otherwise,
*          C is an output variable.  If FACT = 'F' and EQUED = 'C' or
*          'B', each element of C must be positive.
*          C is replicated in every process row, and is aligned with
*          the distributed matrix A.
*
*  B       (local input/local output) DOUBLE PRECISION pointer
*          into the local memory to an array of local dimension
*          (LLD_B,LOCc(JB+NRHS-1) ).  On entry, the N-by-NRHS right-hand
*          side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
*          EQUED = 'N', B(IB:IB+N-1,JB:JB+NRHS-1) is not modified; if
*          TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
*          diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = 'T' or 'C'
*          and EQUED = 'C' or 'B', B(IB:IB+N-1,JB:JB+NRHS-1) is over-
*          written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).
*
*  IB      (global input) INTEGER
*          The row index in the global array B indicating the first
*          row of sub( B ).
*
*  JB      (global input) INTEGER
*          The column index in the global array B indicating the
*          first column of sub( B ).
*
*  DESCB   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix B.
*
*  X       (local input/local output) DOUBLE PRECISION pointer
*          into the local memory to an array of local dimension
*          (LLD_X, LOCc(JX+NRHS-1)).  If INFO = 0, the N-by-NRHS
*          solution matrix X(IX:IX+N-1,JX:JX+NRHS-1) to the original
*          system of equations.  Note that A(IA:IA+N-1,JA:JA+N-1) and
*          B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit if
*          EQUED .ne. 'N', and the solution to the equilibrated system
*          is inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'N'
*          and EQUED = 'C' or 'B', or
*          inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'T' or 'C'
*          and EQUED = 'R' or 'B'.
*
*  IX      (global input) INTEGER
*          The row index in the global array X indicating the first
*          row of sub( X ).
*
*  JX      (global input) INTEGER
*          The column index in the global array X indicating the
*          first column of sub( X ).
*
*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix X.
*
*  RCOND   (global output) DOUBLE PRECISION
*          The estimate of the reciprocal condition number of the matrix
*          A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done).  If
*          RCOND is less than the machine precision (in particular, if
*          RCOND = 0), the matrix is singular to working precision.
*          This condition is indicated by a return code of INFO > 0.
*
*  FERR    (local output) DOUBLE PRECISION array, dimension LOCc(N_B)
*          The estimated forward error bounds for each solution vector
*          X(j) (the j-th column of the solution matrix
*          X(IX:IX+N-1,JX:JX+NRHS-1). If XTRUE is the true solution,
*          FERR(j) bounds the magnitude of the largest entry in
*          (X(j) - XTRUE) divided by the magnitude of the largest entry
*          in X(j).  The estimate is as reliable as the estimate for
*          RCOND, and is almost always a slight overestimate of the
*          true error.  FERR is replicated in every process row, and is
*          aligned with the matrices B and X.
*
*  BERR    (local output) DOUBLE PRECISION array, dimension LOCc(N_B).
*          The componentwise relative backward error of each solution
*          vector X(j) (i.e., the smallest relative change in
*          any entry of A(IA:IA+N-1,JA:JA+N-1) or
*          B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an exact solution).
*          BERR is replicated in every process row, and is aligned
*          with the matrices B and X.
*
*  WORK    (local workspace/local output) DOUBLE PRECISION array,
*                                                    dimension (LWORK)
*          On exit, WORK(1) returns the minimal and optimal LWORK.
*
*  LWORK   (local or global input) INTEGER
*          The dimension of the array WORK.
*          LWORK is local input and must be at least
*          LWORK = MAX( PDGECON( LWORK ), PDGERFS( LWORK ) )
*                  + LOCr( N_A ).
*
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IWORK   (local workspace/local output) INTEGER array,
*                                                  dimension (LIWORK)
*          On exit, IWORK(1) returns the minimal and optimal LIWORK.
*
*  LIWORK  (local or global input) INTEGER
*          The dimension of the array IWORK.
*          LIWORK is local input and must be at least
*          LIWORK = LOCr(N_A).
*
*          If LIWORK = -1, then LIWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, and i is
*                <= N:  U(IA+I-1,IA+I-1) 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.
*                = N+1: RCOND is less than machine precision.  The
*                       factorization has been completed, but the
*                       matrix is singular to working precision, and
*                       the solution and error bounds have not been
*                       computed.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
     $                     ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                     rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0d+0, zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLEQU, EQUIL, LQUERY, NOFACT, NOTRAN, ROWEQU
      CHARACTER          NORM
      INTEGER            CONWRK, I, IACOL, IAROW, IAFROW, IBROW, IBCOL,
     $                   icoffa, icoffb, icoffx, ictxt, idumm,
     $                   iia, iib, iix,
     $                   infequ, iroffa, iroffaf, iroffb,
     $                   iroffx, ixcol, ixrow, j, jja, jjb, jjx,
     $                   lcm, lcmq,
     $                   liwmin, lwmin, mycol, myrow, np, npcol, nprow,
     $                   nq, nqb, nrhsq, rfswrk
      DOUBLE PRECISION   AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
     $                   ROWCND, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            CDESC( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, descset, pchk2mat,
     $                   dgebr2d, dgebs2d, dgamn2d,
     $                   dgamx2d, infog2l, pdcopy, pdgecon,
     $                   pdgeequ, pdgerfs, pdgetrf, pdgetrs,
     $                   pdlacpy, pdlaqge, pxerbla
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL, ILCM, INDXG2P, NUMROC
      DOUBLE PRECISION   PDLAMCH, PDLANGE
      EXTERNAL           iceil, ilcm, indxg2p, lsame, numroc, pdlange,
     $                   pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, ichar, max, min, mod
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ictxt = desca( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
*     Test the input parameters
*
      info = 0
      IF( nprow.EQ.-1 ) THEN
         info = -(800+ctxt_)
      ELSE
         CALL chk1mat( n, 3, n, 3, ia, ja, desca, 8, info )
         IF( lsame( fact, 'F' ) )
     $      CALL chk1mat( n, 3, n, 3, iaf, jaf, descaf, 12, info )
         CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 20, info )
         CALL chk1mat( n, 3, nrhs, 4, ix, jx, descx, 24, info )
         nofact = lsame( fact, 'N' )
         equil = lsame( fact, 'E' )
         notran = lsame( trans, 'N' )
         IF( nofact .OR. equil ) THEN
            equed = 'N'
            rowequ = .false.
            colequ = .false.
         ELSE
            rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
            colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
            smlnum = pdlamch( ictxt, 'Safe minimum' )
            bignum = one / smlnum
         END IF
         IF( info.EQ.0 ) THEN
            iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
     $                       nprow )
            iafrow = indxg2p( iaf, descaf( mb_ ), myrow,
     $                        descaf( rsrc_ ), nprow )
            ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
     $                       nprow )
            ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
     $                       nprow )
            iroffa = mod( ia-1, desca( mb_ ) )
            iroffaf = mod( iaf-1, descaf( mb_ ) )
            icoffa = mod( ja-1, desca( nb_ ) )
            iroffb = mod( ib-1, descb( mb_ ) )
            icoffb = mod( jb-1, descb( nb_ ) )
            iroffx = mod( ix-1, descx( mb_ ) )
            icoffx = mod( jx-1, descx( nb_ ) )
            CALL infog2l( ia, ja, desca, nprow, npcol, myrow,
     $                    mycol, iia, jja, iarow, iacol )
            np = numroc( n+iroffa, desca( mb_ ), myrow, iarow,
     $                   nprow )
            IF( myrow.EQ.iarow )
     $         np = np-iroffa
            nq = numroc( n+icoffa, desca( nb_ ), mycol, iacol,
     $                   npcol )
            IF( mycol.EQ.iacol )
     $         nq = nq-icoffa
            nqb = iceil( n+iroffa, desca( nb_ )*npcol )
            lcm = ilcm( nprow, npcol )
            lcmq = lcm / npcol
            conwrk = 2*np + 2*nq + max( 2, max( desca( nb_ )*
     $               max( 1, iceil( nprow-1, npcol ) ), nq +
     $               desca( nb_ )*
     $               max( 1, iceil( npcol-1, nprow ) ) ) )
            rfswrk = 3*np
            IF( lsame( trans, 'N' ) ) THEN
               rfswrk = rfswrk + np + nq +
     $                 iceil( nqb, lcmq )*desca( nb_ )
            ELSE IF( lsame( trans, 'T' ).OR.lsame( trans, 'C' ) ) THEN
               rfswrk = rfswrk + np + nq
            END IF
            lwmin = max( conwrk, rfswrk )
            work( 1 ) = dble( lwmin )
            liwmin = np
            iwork( 1 ) = liwmin
            IF( .NOT.nofact .AND. .NOT.equil .AND.
     $          .NOT.lsame( fact, 'F' ) ) THEN
               info = -1
            ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND.
     $         .NOT. lsame( trans, 'C' ) ) THEN
               info = -2
            ELSE IF( iroffa.NE.0 ) THEN
               info = -6
            ELSE IF( icoffa.NE.0 .OR. iroffa.NE.icoffa ) THEN
               info = -7
            ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
               info = -(800+nb_)
            ELSE IF( iafrow.NE.iarow ) THEN
               info = -10
            ELSE IF( iroffaf.NE.0 ) THEN
               info = -10
            ELSE IF( ictxt.NE.descaf( ctxt_ ) ) THEN
               info = -(1200+ctxt_)
            ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
     $         ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
               info = -13
            ELSE
               IF( rowequ ) THEN
                  rcmin = bignum
                  rcmax = zero
                  DO 10 j = iia, iia + np - 1
                     rcmin = min( rcmin, r( j ) )
                     rcmax = max( rcmax, r( j ) )
   10             CONTINUE
                  CALL dgamn2d( ictxt, 'Columnwise', ' ', 1, 1, rcmin,
     $                          1, idumm, idumm, -1, -1, mycol )
                  CALL dgamx2d( ictxt, 'Columnwise', ' ', 1, 1, rcmax,
     $                          1, idumm, idumm, -1, -1, mycol )
                  IF( rcmin.LE.zero ) THEN
                     info = -14
                  ELSE IF( n.GT.0 ) THEN
                     rowcnd = max( rcmin, smlnum ) /
     $                        min( rcmax, bignum )
                  ELSE
                     rowcnd = one
                  END IF
               END IF
               IF( colequ .AND. info.EQ.0 ) THEN
                  rcmin = bignum
                  rcmax = zero
                  DO 20 j = jja, jja+nq-1
                     rcmin = min( rcmin, c( j ) )
                     rcmax = max( rcmax, c( j ) )
   20             CONTINUE
                  CALL dgamn2d( ictxt, 'Rowwise', ' ', 1, 1, rcmin,
     $                          1, idumm, idumm, -1, -1, mycol )
                  CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, rcmax,
     $                          1, idumm, idumm, -1, -1, mycol )
                  IF( rcmin.LE.zero ) THEN
                     info = -15
                  ELSE IF( n.GT.0 ) THEN
                     colcnd = max( rcmin, smlnum ) /
     $                        min( rcmax, bignum )
                  ELSE
                     colcnd = one
                  END IF
               END IF
            END IF
         END IF
*
         work( 1 ) = dble( lwmin )
         iwork( 1 ) = liwmin
         lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
         IF( info.EQ.0 ) THEN
            IF( ibrow.NE.iarow ) THEN
               info = -18
            ELSE IF( ixrow.NE.ibrow ) THEN
               info = -22
            ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
               info = -(2000+nb_)
            ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
               info = -(2000+ctxt_)
            ELSE IF( descx( mb_ ).NE.desca( nb_ ) ) THEN
               info = -(2400+nb_)
            ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
               info = -(2400+ctxt_)
            ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
               info = -29
            ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
               info = -31
            END IF
            idum1( 1 ) = ichar( fact )
            idum2( 1 ) = 1
            idum1( 2 ) = ichar( trans )
            idum2( 2 ) = 2
            IF( lsame( fact, 'F' ) ) THEN
               idum1( 3 ) = ichar( equed )
               idum2( 3 ) = 14
               IF( lwork.EQ.-1 ) THEN
                  idum1( 4 ) = -1
               ELSE
                  idum1( 4 ) = 1
               END IF
               idum2( 4 ) = 29
               IF( liwork.EQ.-1 ) THEN
                  idum1( 5 ) = -1
               ELSE
                  idum1( 5 ) = 1
               END IF
               idum2( 5 ) = 31
               CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3,
     $                        nrhs, 4, ib, jb, descb, 20, 5, idum1,
     $                        idum2, info )
            ELSE
               IF( lwork.EQ.-1 ) THEN
                  idum1( 3 ) = -1
               ELSE
                  idum1( 3 ) = 1
               END IF
               idum2( 3 ) = 29
               IF( liwork.EQ.-1 ) THEN
                  idum1( 4 ) = -1
               ELSE
                  idum1( 4 ) = 1
               END IF
               idum2( 4 ) = 31
               CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3,
     $                        nrhs, 4, ib, jb, descb, 20, 4, idum1,
     $                        idum2, info )
            END IF
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL pxerbla( ictxt, 'PDGESVX', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
      IF( equil ) THEN
*
*        Compute row and column scalings to equilibrate the matrix A.
*
         CALL pdgeequ( n, n, a, ia, ja, desca, r, c, rowcnd, colcnd,
     $                 amax, infequ )
         IF( infequ.EQ.0 ) THEN
*
*           Equilibrate the matrix.
*
            CALL pdlaqge( n, n, a, ia, ja, desca, r, c, rowcnd, colcnd,
     $                    amax, equed )
            rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
            colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
         END IF
      END IF
*
*     Scale the right-hand side.
*
      CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol, iib,
     $              jjb, ibrow, ibcol )
      np = numroc( n+iroffb, descb( mb_ ), myrow, ibrow, nprow )
      nrhsq = numroc( nrhs+icoffb, descb( nb_ ), mycol, ibcol, npcol )
      IF( myrow.EQ.ibrow )
     $   np = np-iroffb
      IF( mycol.EQ.ibcol )
     $   nrhsq = nrhsq-icoffb
*
      IF( notran ) THEN
         IF( rowequ ) THEN
            DO 40 j = jjb, jjb+nrhsq-1
               DO 30 i = iib, iib+np-1
                  b( i+( j-1 )*descb( lld_ ) ) = r( i )*
     $                                     b( i+( j-1 )*descb( lld_ ) )
   30          CONTINUE
   40       CONTINUE
         END IF
      ELSE IF( colequ ) THEN
*
*        Transpose the Column scale factors
*
         CALL descset( cdesc, 1, n+icoffa, 1, desca( nb_ ), myrow,
     $                 iacol, ictxt, 1 )
         CALL pdcopy( n, c, 1, ja, cdesc, cdesc( lld_ ), work, ib, jb,
     $                descb, 1 )
         IF( mycol.EQ.ibcol ) THEN
            CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( iib ),
     $                    descb( lld_ ) )
         ELSE
            CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( iib ),
     $                    descb( lld_ ), myrow, ibcol )
         END IF
         DO 60 j = jjb, jjb+nrhsq-1
            DO 50 i = iib, iib+np-1
               b( i+( j-1 )*descb( lld_ ) ) = work( i )*
     $                                     b( i+( j-1 )*descb( lld_ ) )
   50       CONTINUE
   60    CONTINUE
      END IF
*
      IF( nofact.OR.equil ) THEN
*
*        Compute the LU factorization of A.
*
         CALL pdlacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
     $                 descaf )
         CALL pdgetrf( n, n, af, iaf, jaf, descaf, ipiv, info )
*
*        Return if INFO is non-zero.
*
         IF( info.NE.0 ) THEN
            IF( info.GT.0 )
     $         rcond = zero
            RETURN
         END IF
      END IF
*
*     Compute the norm of the matrix A.
*
      IF( notran ) THEN
         norm = '1'
      ELSE
         norm = 'I'
      END IF
      anorm = pdlange( norm, n, n, a, ia, ja, desca, work )
*
*     Compute the reciprocal of the condition number of A.
*
      CALL pdgecon( norm, n, af, iaf, jaf, descaf, anorm, rcond, work,
     $              lwork, iwork, liwork, info )
*
*     Return if the matrix is singular to working precision.
*
      IF( rcond.LT.pdlamch( ictxt, 'Epsilon' ) ) THEN
         info = ia + n
         RETURN
      END IF
*
*     Compute the solution matrix X.
*
      CALL pdlacpy( 'Full', n, nrhs, b, ib, jb, descb, x, ix, jx,
     $              descx )
      CALL pdgetrs( trans, n, nrhs, af, iaf, jaf, descaf, ipiv, x, ix,
     $              jx, descx, info )
*
*     Use iterative refinement to improve the computed solution and
*     compute error bounds and backward error estimates for it.
*
      CALL pdgerfs( trans, n, nrhs, a, ia, ja, desca, af, iaf, jaf,
     $              descaf, ipiv, b, ib, jb, descb, x, ix, jx, descx,
     $              ferr, berr, work, lwork, iwork, liwork, info )
*
*     Transform the solution matrix X to a solution of the original
*     system.
*
      CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix,
     $              jjx, ixrow, ixcol )
      np = numroc( n+iroffx, descx( mb_ ), myrow, ixrow, nprow )
      nrhsq = numroc( nrhs+icoffx, descx( nb_ ), mycol, ixcol, npcol )
      IF( myrow.EQ.ibrow )
     $   np = np-iroffx
      IF( mycol.EQ.ibcol )
     $   nrhsq = nrhsq-icoffx
*
      IF( notran ) THEN
         IF( colequ ) THEN
*
*           Transpose the column scaling factors
*
            CALL descset( cdesc, 1, n+icoffa, 1, desca( nb_ ), myrow,
     $                    iacol, ictxt, 1 )
            CALL pdcopy( n, c, 1, ja, cdesc, cdesc( lld_ ), work, ix,
     $                   jx, descx, 1 )
            IF( mycol.EQ.ibcol ) THEN
                CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1,
     $                        work( iix ), descx( lld_ ) )
            ELSE
                CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1,
     $                        work( iix ), descx( lld_ ), myrow, ibcol )
            END IF
*
            DO 80 j = jjx, jjx+nrhsq-1
               DO 70 i = iix, iix+np-1
                  x( i+( j-1 )*descx( lld_ ) ) = work( i )*
     $                                      x( i+( j-1 )*descx( lld_ ) )
   70          CONTINUE
   80       CONTINUE
            DO 90 j = jjx, jjx+nrhsq-1
               ferr( j ) = ferr( j ) / colcnd
   90       CONTINUE
         END IF
      ELSE IF( rowequ ) THEN
         DO 110 j = jjx, jjx+nrhsq-1
            DO 100 i = iix, iix+np-1
               x( i+( j-1 )*descx( lld_ ) ) = r( i )*
     $                                     x( i+( j-1 )*descx( lld_ ) )
  100       CONTINUE
  110    CONTINUE
         DO 120 j = jjx, jjx+nrhsq-1
            ferr( j ) = ferr( j ) / rowcnd
  120    CONTINUE
      END IF
*
      work( 1 ) = dble( lwmin )
      iwork( 1 ) = liwmin
*
      RETURN
*
*     End of PDGESVX
*
      END