pdposvx.f

Go to the documentation of this file.

      SUBROUTINE pdposvx( FACT, UPLO, N, NRHS, A, IA, JA, DESCA, AF,
     $                    IAF, JAF, DESCAF, EQUED, SR, SC, 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, UPLO
      INTEGER            IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
     $                   LWORK, N, NRHS
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCAF( * ), DESCB( * ),
     $                   DESCX( * ), IWORK( * )
      DOUBLE PRECISION   A( * ), AF( * ),
     $                   b( * ), berr( * ), ferr( * ),
     $                   sc( * ), sr( * ), work( * ), x( * )
*     ..
*
*  Purpose
*  =======
*
*  PDPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T 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.  In the following comments Y denotes Y(IY:IY+M-1,JY:JY+K-1)
*  a M-by-K matrix where Y can be A, AF, B and X.
*
*  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
*  ===========
*
*  The following steps are performed:
*
*  1. If FACT = 'E', real scaling factors are computed to equilibrate
*     the system:
*        diag(SR) * A * diag(SC) * inv(diag(SC)) * X = diag(SR) * 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(SR)*A*diag(SC) and B by diag(SR)*B.
*
*  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
*     factor the matrix A (after equilibration if FACT = 'E') as
*        A = U**T* U,  if UPLO = 'U', or
*        A = L * L**T,  if UPLO = 'L',
*     where U is an upper triangular matrix and L is a lower triangular
*     matrix.
*
*  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 equilibration was used, the matrix X is premultiplied by
*     diag(SR) 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 is
*          supplied on entry, and if not, whether the matrix A should be
*          equilibrated before it is factored.
*          = 'F':  On entry, AF contains the factored form of A.
*                  If EQUED = 'Y', the matrix A has been equilibrated
*                  with scaling factors given by S.  A and AF will not
*                  be modified.
*          = 'N':  The matrix A will be copied to AF and factored.
*          = 'E':  The matrix A will be equilibrated if necessary, then
*                  copied to AF and factored.
*
*  UPLO    (global input) CHARACTER
*          = 'U':  Upper triangle of A is stored;
*          = 'L':  Lower triangle of A is stored.
*
*  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 and X.  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 symmetric matrix A, except if FACT = 'F' and
*          EQUED = 'Y', then A must contain the equilibrated matrix
*          diag(SR)*A*diag(SC).  If UPLO = 'U', the leading
*          N-by-N upper triangular part of A contains the upper
*          triangular part of the matrix A, and the strictly lower
*          triangular part of A is not referenced.  If UPLO = 'L', the
*          leading N-by-N lower triangular part of A contains the lower
*          triangular part of the matrix A, and the strictly upper
*          triangular part of A is not referenced.  A is not modified if
*          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
*
*          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
*          diag(SR)*A*diag(SC).
*
*  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 is an input argument and on entry
*          contains the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T, in the same storage
*          format as A.  If EQUED .ne. 'N', then AF is the factored form
*          of the equilibrated matrix diag(SR)*A*diag(SC).
*
*          If FACT = 'N', then AF is an output argument and on exit
*          returns the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T of the original
*          matrix A.
*
*          If FACT = 'E', then AF is an output argument and on exit
*          returns the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T of the equilibrated
*          matrix A (see the description of A 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.
*
*  EQUED   (global input/global output) CHARACTER
*          Specifies the form of equilibration that was done.
*          = 'N':  No equilibration (always true if FACT = 'N').
*          = 'Y':  Equilibration was done, i.e., A has been replaced by
*                  diag(SR) * A * diag(SC).
*          EQUED is an input variable if FACT = 'F'; otherwise, it is an
*          output variable.
*
*  SR      (local input/local output) DOUBLE PRECISION array,
*                                                    dimension (LLD_A)
*          The scale factors for A distributed across process rows;
*          not accessed if EQUED = 'N'.  SR is an input variable if
*          FACT = 'F'; otherwise, SR is an output variable.
*          If FACT = 'F' and EQUED = 'Y', each element of SR must be
*          positive.
*
*  SC      (local input/local output) DOUBLE PRECISION array,
*                                              dimension (LOC(N_A))
*          The scale factors for A distributed across
*          process columns; not accessed if EQUED = 'N'. SC is an input
*          variable if FACT = 'F'; otherwise, SC is an output variable.
*          If FACT = 'F' and EQUED = 'Y', each element of SC must be
*          positive.
*
*  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.
*          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.
*
*  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 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(SC))*X if TRANS = 'N' and EQUED = 'C' or
*          'B', or inv(diag(SR))*X 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 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, and the solution and
*          error bounds are not computed.
*
*  FERR    (local output) DOUBLE PRECISION array, dimension (LOC(N_B))
*          The estimated forward error bounds for each solution vector
*          X(j) (the j-th column of the solution matrix X).
*          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 quality of
*          the error bound depends on the quality of the estimate of
*          norm(inv(A)) computed in the code; if the estimate of
*          norm(inv(A)) is accurate, the error bound is guaranteed.
*
*  BERR    (local output) DOUBLE PRECISION array, dimension (LOC(N_B))
*          The componentwise relative backward error of each solution
*          vector X(j) (i.e., the smallest relative change in
*          any entry of A or B that makes X(j) an exact solution).
*
*  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( PDPOCON( LWORK ), PDPORFS( LWORK ) )
*                  + LOCr( N_A ).
*          LWORK = 3*DESCA( LLD_ )
*
*          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 = DESCA( LLD_ )
*          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: if INFO = i, the leading minor of order i of A
*                     is not positive definite, so the factorization
*                     could not be completed, and 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            EQUIL, LQUERY, NOFACT, RCEQU
      INTEGER            I, IACOL, IAROW, IAFROW, IBROW, IBCOL, ICOFF,
     $                   ICOFFA, ICTXT, IDUMM, IIA, IIB, IIX, INFEQU,
     $                   iroff, iroffa, iroffaf, iroffb, iroffx, ixcol,
     $                   ixrow, j, jja, jjb, jjx, ldb, ldx, liwmin,
     $                   lwmin, mycol, myrow, np, npcol, nprow, nrhsq,
     $                   nq
      DOUBLE PRECISION   AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            IDUM1( 5 ), IDUM2( 5 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, pchk2mat,
     $                   dgamn2d, dgamx2d, infog2l,
     $                   pdpocon, pdpoequ, pdporfs,
     $                   pdpotrf,
     $                   pdpotrs, pdlacpy, pdlaqsy, pxerbla
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2P, NUMROC
      DOUBLE PRECISION   PDLAMCH, PDLANSY
      EXTERNAL           pdlamch, indxg2p, lsame, numroc, pdlansy
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          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 )
         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_ ) )
            iroffx = mod( ix-1, descx( mb_ ) )
            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
            lwmin = 3*desca( lld_ )
            liwmin = np
            nofact = lsame( fact, 'N' )
            equil = lsame( fact, 'E' )
            IF( nofact .OR. equil ) THEN
               equed = 'N'
               rcequ = .false.
            ELSE
               rcequ = lsame( equed, 'Y' )
               smlnum = pdlamch( ictxt, 'Safe minimum' )
               bignum = one / smlnum
            END IF
            IF( .NOT.nofact .AND. .NOT.equil .AND.
     $          .NOT.lsame( fact, 'F' ) ) THEN
               info = -1
            ELSE IF( .NOT.lsame( uplo, 'U' ) .AND.
     $               .NOT.lsame( uplo, 'L' ) ) 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.
     $         ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
               info = -13
            ELSE
               IF( rcequ ) THEN
*
                  smin = bignum
                  smax = zero
                  DO 10 j = iia, iia + np - 1
                     smin = min( smin, sr( j ) )
                     smax = max( smax, sr( j ) )
   10             CONTINUE
                  CALL dgamn2d( ictxt, 'Columnwise', ' ', 1, 1, smin,
     $                          1, idumm, idumm, -1, -1, mycol )
                  CALL dgamx2d( ictxt, 'Columnwise', ' ', 1, 1, smax,
     $                          1, idumm, idumm, -1, -1, mycol )
                  IF( smin.LE.zero ) THEN
                     info = -14
                  ELSE IF( n.GT.0 ) THEN
                     scond = max( smin, smlnum ) / min( smax, bignum )
                  ELSE
                     scond = 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( ictxt.NE.descx( ctxt_ ) ) THEN
               info = -(2400+ctxt_)
            ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
               info = -28
            ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
               info = -30
            END IF
            idum1( 1 ) = ichar( fact )
            idum2( 1 ) = 1
            idum1( 2 ) = ichar( uplo )
            idum2( 2 ) = 2
            IF( lsame( fact, 'F' ) ) THEN
               idum1( 3 ) = ichar( equed )
               idum2( 3 ) = 13
               IF( lwork.EQ.-1 ) THEN
                  idum1( 4 ) = -1
               ELSE
                  idum1( 4 ) = 1
               END IF
               idum2( 4 ) = 28
               IF( liwork.EQ.-1 ) THEN
                  idum1( 5 ) = -1
               ELSE
                  idum1( 5 ) = 1
               END IF
               idum2( 5 ) = 30
               CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
     $                        4, ib, jb, descb, 19, 5, idum1, idum2,
     $                        info )
            ELSE
               IF( lwork.EQ.-1 ) THEN
                  idum1( 3 ) = -1
               ELSE
                  idum1( 3 ) = 1
               END IF
               idum2( 3 ) = 28
               IF( liwork.EQ.-1 ) THEN
                  idum1( 4 ) = -1
               ELSE
                  idum1( 4 ) = 1
               END IF
               idum2( 4 ) = 30
               CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
     $                        4, ib, jb, descb, 19, 4, idum1, idum2,
     $                        info )
            END IF
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL pxerbla( ictxt, 'PDPOSVX', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
      IF( equil ) THEN
*
*        Compute row and column scalings to equilibrate the matrix A.
*
         CALL pdpoequ( n, a, ia, ja, desca, sr, sc, scond, amax,
     $                 infequ )
*
         IF( infequ.EQ.0 ) THEN
*
*           Equilibrate the matrix.
*
            CALL pdlaqsy( uplo, n, a, ia, ja, desca, sr, sc, scond,
     $                    amax, equed )
            rcequ = lsame( equed, 'Y' )
         END IF
      END IF
*
*     Scale the right-hand side.
*
      CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol, iib,
     $              jjb, ibrow, ibcol )
      ldb =  descb( lld_ )
      iroff = mod( ib-1, descb( mb_ ) )
      icoff = mod( jb-1, descb( nb_ ) )
      np = numroc( n+iroff, descb( mb_ ), myrow, ibrow, nprow )
      nrhsq = numroc( nrhs+icoff, descb( nb_ ), mycol, ibcol, npcol )
      IF( myrow.EQ.ibrow ) np = np-iroff
      IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
*
      IF( rcequ ) THEN
         DO 30 j = jjb, jjb+nrhsq-1
            DO 20 i = iib, iib+np-1
               b( i + ( j-1 )*ldb ) = sr( i )*b( i + ( j-1 )*ldb )
   20       CONTINUE
   30    CONTINUE
      END IF
*
      IF( nofact .OR. equil ) THEN
*
*        Compute the Cholesky factorization A = U'*U or A = L*L'.
*
         CALL pdlacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
     $                 descaf )
         CALL pdpotrf( uplo, n, af, iaf, jaf, descaf, 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.
*
      anorm = pdlansy( '1', uplo, n, a, ia, ja, desca, work )
*
*     Compute the reciprocal of the condition number of A.
*
      CALL pdpocon( uplo, 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 pdpotrs( uplo, n, nrhs, af, iaf, jaf, descaf, x, ix, jx,
     $              descx, info )
*
*     Use iterative refinement to improve the computed solution and
*     compute error bounds and backward error estimates for it.
*
      CALL pdporfs( uplo, n, nrhs, a, ia, ja, desca, af, iaf, jaf,
     $              descaf, 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 )
      ldx = descx( lld_ )
      iroff = mod( ix-1, descx( mb_ ) )
      icoff = mod( jx-1, descx( nb_ ) )
      np = numroc( n+iroff, descx( mb_ ), myrow, ixrow, nprow )
      nrhsq = numroc( nrhs+icoff, descx( nb_ ), mycol, ixcol, npcol )
      IF( myrow.EQ.ibrow ) np = np-iroff
      IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
*
      IF( rcequ ) THEN
         DO 50 j = jjx, jjx+nrhsq-1
            DO 40 i = iix, iix+np-1
               x( i + ( j-1 )*ldx ) = sr( i )*x( i + ( j-1 )*ldx )
   40       CONTINUE
   50    CONTINUE
         DO 60 j = jjx, jjx+nrhsq-1
            ferr( j ) = ferr( j ) / scond
   60    CONTINUE
      END IF
*
      work( 1 ) = dble( lwmin )
      iwork( 1 ) = liwmin
      RETURN
*
*     End of PDPOSVX
*
      END