SRC/dgesvx.f

   1 *> \brief <b> DGESVX computes the solution to system of linear equations A * X = B for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGESVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
  22 *                          EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR,
  23 *                          WORK, IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          EQUED, FACT, TRANS
  27 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
  28 *       DOUBLE PRECISION   RCOND
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IPIV( * ), IWORK( * )
  32 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  33 *      $                   BERR( * ), C( * ), FERR( * ), R( * ),
  34 *      $                   WORK( * ), X( LDX, * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> DGESVX uses the LU factorization to compute the solution to a real
  44 *> system of linear equations
  45 *>    A * X = B,
  46 *> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
  47 *>
  48 *> Error bounds on the solution and a condition estimate are also
  49 *> provided.
  50 *> \endverbatim
  51 *
  52 *> \par Description:
  53 *  =================
  54 *>
  55 *> \verbatim
  56 *>
  57 *> The following steps are performed:
  58 *>
  59 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
  60 *>    the system:
  61 *>       TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
  62 *>       TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
  63 *>       TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
  64 *>    Whether or not the system will be equilibrated depends on the
  65 *>    scaling of the matrix A, but if equilibration is used, A is
  66 *>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
  67 *>    or diag(C)*B (if TRANS = 'T' or 'C').
  68 *>
  69 *> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
  70 *>    matrix A (after equilibration if FACT = 'E') as
  71 *>       A = P * L * U,
  72 *>    where P is a permutation matrix, L is a unit lower triangular
  73 *>    matrix, and U is upper triangular.
  74 *>
  75 *> 3. If some U(i,i)=0, so that U is exactly singular, then the routine
  76 *>    returns with INFO = i. Otherwise, the factored form of A is used
  77 *>    to estimate the condition number of the matrix A.  If the
  78 *>    reciprocal of the condition number is less than machine precision,
  79 *>    INFO = N+1 is returned as a warning, but the routine still goes on
  80 *>    to solve for X and compute error bounds as described below.
  81 *>
  82 *> 4. The system of equations is solved for X using the factored form
  83 *>    of A.
  84 *>
  85 *> 5. Iterative refinement is applied to improve the computed solution
  86 *>    matrix and calculate error bounds and backward error estimates
  87 *>    for it.
  88 *>
  89 *> 6. If equilibration was used, the matrix X is premultiplied by
  90 *>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
  91 *>    that it solves the original system before equilibration.
  92 *> \endverbatim
  93 *
  94 *  Arguments:
  95 *  ==========
  96 *
  97 *> \param[in] FACT
  98 *> \verbatim
  99 *>          FACT is CHARACTER*1
 100 *>          Specifies whether or not the factored form of the matrix A is
 101 *>          supplied on entry, and if not, whether the matrix A should be
 102 *>          equilibrated before it is factored.
 103 *>          = 'F':  On entry, AF and IPIV contain the factored form of A.
 104 *>                  If EQUED is not 'N', the matrix A has been
 105 *>                  equilibrated with scaling factors given by R and C.
 106 *>                  A, AF, and IPIV are not modified.
 107 *>          = 'N':  The matrix A will be copied to AF and factored.
 108 *>          = 'E':  The matrix A will be equilibrated if necessary, then
 109 *>                  copied to AF and factored.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] TRANS
 113 *> \verbatim
 114 *>          TRANS is CHARACTER*1
 115 *>          Specifies the form of the system of equations:
 116 *>          = 'N':  A * X = B     (No transpose)
 117 *>          = 'T':  A**T * X = B  (Transpose)
 118 *>          = 'C':  A**H * X = B  (Transpose)
 119 *> \endverbatim
 120 *>
 121 *> \param[in] N
 122 *> \verbatim
 123 *>          N is INTEGER
 124 *>          The number of linear equations, i.e., the order of the
 125 *>          matrix A.  N >= 0.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] NRHS
 129 *> \verbatim
 130 *>          NRHS is INTEGER
 131 *>          The number of right hand sides, i.e., the number of columns
 132 *>          of the matrices B and X.  NRHS >= 0.
 133 *> \endverbatim
 134 *>
 135 *> \param[in,out] A
 136 *> \verbatim
 137 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 138 *>          On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
 139 *>          not 'N', then A must have been equilibrated by the scaling
 140 *>          factors in R and/or C.  A is not modified if FACT = 'F' or
 141 *>          'N', or if FACT = 'E' and EQUED = 'N' on exit.
 142 *>
 143 *>          On exit, if EQUED .ne. 'N', A is scaled as follows:
 144 *>          EQUED = 'R':  A := diag(R) * A
 145 *>          EQUED = 'C':  A := A * diag(C)
 146 *>          EQUED = 'B':  A := diag(R) * A * diag(C).
 147 *> \endverbatim
 148 *>
 149 *> \param[in] LDA
 150 *> \verbatim
 151 *>          LDA is INTEGER
 152 *>          The leading dimension of the array A.  LDA >= max(1,N).
 153 *> \endverbatim
 154 *>
 155 *> \param[in,out] AF
 156 *> \verbatim
 157 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
 158 *>          If FACT = 'F', then AF is an input argument and on entry
 159 *>          contains the factors L and U from the factorization
 160 *>          A = P*L*U as computed by DGETRF.  If EQUED .ne. 'N', then
 161 *>          AF is the factored form of the equilibrated matrix A.
 162 *>
 163 *>          If FACT = 'N', then AF is an output argument and on exit
 164 *>          returns the factors L and U from the factorization A = P*L*U
 165 *>          of the original matrix A.
 166 *>
 167 *>          If FACT = 'E', then AF is an output argument and on exit
 168 *>          returns the factors L and U from the factorization A = P*L*U
 169 *>          of the equilibrated matrix A (see the description of A for
 170 *>          the form of the equilibrated matrix).
 171 *> \endverbatim
 172 *>
 173 *> \param[in] LDAF
 174 *> \verbatim
 175 *>          LDAF is INTEGER
 176 *>          The leading dimension of the array AF.  LDAF >= max(1,N).
 177 *> \endverbatim
 178 *>
 179 *> \param[in,out] IPIV
 180 *> \verbatim
 181 *>          IPIV is INTEGER array, dimension (N)
 182 *>          If FACT = 'F', then IPIV is an input argument and on entry
 183 *>          contains the pivot indices from the factorization A = P*L*U
 184 *>          as computed by DGETRF; row i of the matrix was interchanged
 185 *>          with row IPIV(i).
 186 *>
 187 *>          If FACT = 'N', then IPIV is an output argument and on exit
 188 *>          contains the pivot indices from the factorization A = P*L*U
 189 *>          of the original matrix A.
 190 *>
 191 *>          If FACT = 'E', then IPIV is an output argument and on exit
 192 *>          contains the pivot indices from the factorization A = P*L*U
 193 *>          of the equilibrated matrix A.
 194 *> \endverbatim
 195 *>
 196 *> \param[in,out] EQUED
 197 *> \verbatim
 198 *>          EQUED is CHARACTER*1
 199 *>          Specifies the form of equilibration that was done.
 200 *>          = 'N':  No equilibration (always true if FACT = 'N').
 201 *>          = 'R':  Row equilibration, i.e., A has been premultiplied by
 202 *>                  diag(R).
 203 *>          = 'C':  Column equilibration, i.e., A has been postmultiplied
 204 *>                  by diag(C).
 205 *>          = 'B':  Both row and column equilibration, i.e., A has been
 206 *>                  replaced by diag(R) * A * diag(C).
 207 *>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
 208 *>          output argument.
 209 *> \endverbatim
 210 *>
 211 *> \param[in,out] R
 212 *> \verbatim
 213 *>          R is DOUBLE PRECISION array, dimension (N)
 214 *>          The row scale factors for A.  If EQUED = 'R' or 'B', A is
 215 *>          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
 216 *>          is not accessed.  R is an input argument if FACT = 'F';
 217 *>          otherwise, R is an output argument.  If FACT = 'F' and
 218 *>          EQUED = 'R' or 'B', each element of R must be positive.
 219 *> \endverbatim
 220 *>
 221 *> \param[in,out] C
 222 *> \verbatim
 223 *>          C is DOUBLE PRECISION array, dimension (N)
 224 *>          The column scale factors for A.  If EQUED = 'C' or 'B', A is
 225 *>          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
 226 *>          is not accessed.  C is an input argument if FACT = 'F';
 227 *>          otherwise, C is an output argument.  If FACT = 'F' and
 228 *>          EQUED = 'C' or 'B', each element of C must be positive.
 229 *> \endverbatim
 230 *>
 231 *> \param[in,out] B
 232 *> \verbatim
 233 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
 234 *>          On entry, the N-by-NRHS right hand side matrix B.
 235 *>          On exit,
 236 *>          if EQUED = 'N', B is not modified;
 237 *>          if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
 238 *>          diag(R)*B;
 239 *>          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
 240 *>          overwritten by diag(C)*B.
 241 *> \endverbatim
 242 *>
 243 *> \param[in] LDB
 244 *> \verbatim
 245 *>          LDB is INTEGER
 246 *>          The leading dimension of the array B.  LDB >= max(1,N).
 247 *> \endverbatim
 248 *>
 249 *> \param[out] X
 250 *> \verbatim
 251 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
 252 *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
 253 *>          to the original system of equations.  Note that A and B are
 254 *>          modified on exit if EQUED .ne. 'N', and the solution to the
 255 *>          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
 256 *>          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
 257 *>          and EQUED = 'R' or 'B'.
 258 *> \endverbatim
 259 *>
 260 *> \param[in] LDX
 261 *> \verbatim
 262 *>          LDX is INTEGER
 263 *>          The leading dimension of the array X.  LDX >= max(1,N).
 264 *> \endverbatim
 265 *>
 266 *> \param[out] RCOND
 267 *> \verbatim
 268 *>          RCOND is DOUBLE PRECISION
 269 *>          The estimate of the reciprocal condition number of the matrix
 270 *>          A after equilibration (if done).  If RCOND is less than the
 271 *>          machine precision (in particular, if RCOND = 0), the matrix
 272 *>          is singular to working precision.  This condition is
 273 *>          indicated by a return code of INFO > 0.
 274 *> \endverbatim
 275 *>
 276 *> \param[out] FERR
 277 *> \verbatim
 278 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 279 *>          The estimated forward error bound for each solution vector
 280 *>          X(j) (the j-th column of the solution matrix X).
 281 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 282 *>          is an estimated upper bound for the magnitude of the largest
 283 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 284 *>          largest element in X(j).  The estimate is as reliable as
 285 *>          the estimate for RCOND, and is almost always a slight
 286 *>          overestimate of the true error.
 287 *> \endverbatim
 288 *>
 289 *> \param[out] BERR
 290 *> \verbatim
 291 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 292 *>          The componentwise relative backward error of each solution
 293 *>          vector X(j) (i.e., the smallest relative change in
 294 *>          any element of A or B that makes X(j) an exact solution).
 295 *> \endverbatim
 296 *>
 297 *> \param[out] WORK
 298 *> \verbatim
 299 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
 300 *>          On exit, WORK(1) contains the reciprocal pivot growth
 301 *>          factor norm(A)/norm(U). The "max absolute element" norm is
 302 *>          used. If WORK(1) is much less than 1, then the stability
 303 *>          of the LU factorization of the (equilibrated) matrix A
 304 *>          could be poor. This also means that the solution X, condition
 305 *>          estimator RCOND, and forward error bound FERR could be
 306 *>          unreliable. If factorization fails with 0<INFO<=N, then
 307 *>          WORK(1) contains the reciprocal pivot growth factor for the
 308 *>          leading INFO columns of A.
 309 *> \endverbatim
 310 *>
 311 *> \param[out] IWORK
 312 *> \verbatim
 313 *>          IWORK is INTEGER array, dimension (N)
 314 *> \endverbatim
 315 *>
 316 *> \param[out] INFO
 317 *> \verbatim
 318 *>          INFO is INTEGER
 319 *>          = 0:  successful exit
 320 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 321 *>          > 0:  if INFO = i, and i is
 322 *>                <= N:  U(i,i) is exactly zero.  The factorization has
 323 *>                       been completed, but the factor U is exactly
 324 *>                       singular, so the solution and error bounds
 325 *>                       could not be computed. RCOND = 0 is returned.
 326 *>                = N+1: U is nonsingular, but RCOND is less than machine
 327 *>                       precision, meaning that the matrix is singular
 328 *>                       to working precision.  Nevertheless, the
 329 *>                       solution and error bounds are computed because
 330 *>                       there are a number of situations where the
 331 *>                       computed solution can be more accurate than the
 332 *>                       value of RCOND would suggest.
 333 *> \endverbatim
 334 *
 335 *  Authors:
 336 *  ========
 337 *
 338 *> \author Univ. of Tennessee
 339 *> \author Univ. of California Berkeley
 340 *> \author Univ. of Colorado Denver
 341 *> \author NAG Ltd.
 342 *
 343 *> \date April 2012
 344 *
 345 *> \ingroup doubleGEsolve
 346 *
 347 *  =====================================================================
 348       SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
 349      $                   EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR,
 350      $                   WORK, IWORK, INFO )
 351 *
 352 *  -- LAPACK driver routine (version 3.4.1) --
 353 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 354 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 355 *     April 2012
 356 *
 357 *     .. Scalar Arguments ..
 358       CHARACTER          EQUED, FACT, TRANS
 359       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
 360       DOUBLE PRECISION   RCOND
 361 *     ..
 362 *     .. Array Arguments ..
 363       INTEGER            IPIV( * ), IWORK( * )
 364       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
 365      $                   BERR( * ), C( * ), FERR( * ), R( * ),
 366      $                   WORK( * ), X( LDX, * )
 367 *     ..
 368 *
 369 *  =====================================================================
 370 *
 371 *     .. Parameters ..
 372       DOUBLE PRECISION   ZERO, ONE
 373       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 374 *     ..
 375 *     .. Local Scalars ..
 376       LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
 377       CHARACTER          NORM
 378       INTEGER            I, INFEQU, J
 379       DOUBLE PRECISION   AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
 380      $                   ROWCND, RPVGRW, SMLNUM
 381 *     ..
 382 *     .. External Functions ..
 383       LOGICAL            LSAME
 384       DOUBLE PRECISION   DLAMCH, DLANGE, DLANTR
 385       EXTERNAL           LSAME, DLAMCH, DLANGE, DLANTR
 386 *     ..
 387 *     .. External Subroutines ..
 388       EXTERNAL           DGECON, DGEEQU, DGERFS, DGETRF, DGETRS, DLACPY,
 389      $                   DLAQGE, XERBLA
 390 *     ..
 391 *     .. Intrinsic Functions ..
 392       INTRINSIC          MAX, MIN
 393 *     ..
 394 *     .. Executable Statements ..
 395 *
 396       INFO = 0
 397       NOFACT = LSAME( FACT, 'N' )
 398       EQUIL = LSAME( FACT, 'E' )
 399       NOTRAN = LSAME( TRANS, 'N' )
 400       IF( NOFACT .OR. EQUIL ) THEN
 401          EQUED = 'N'
 402          ROWEQU = .FALSE.
 403          COLEQU = .FALSE.
 404       ELSE
 405          ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 406          COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 407          SMLNUM = DLAMCH( 'Safe minimum' )
 408          BIGNUM = ONE / SMLNUM
 409       END IF
 410 *
 411 *     Test the input parameters.
 412 *
 413       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
 414      $     THEN
 415          INFO = -1
 416       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 417      $         LSAME( TRANS, 'C' ) ) THEN
 418          INFO = -2
 419       ELSE IF( N.LT.0 ) THEN
 420          INFO = -3
 421       ELSE IF( NRHS.LT.0 ) THEN
 422          INFO = -4
 423       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 424          INFO = -6
 425       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 426          INFO = -8
 427       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
 428      $         ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
 429          INFO = -10
 430       ELSE
 431          IF( ROWEQU ) THEN
 432             RCMIN = BIGNUM
 433             RCMAX = ZERO
 434             DO 10 J = 1, N
 435                RCMIN = MIN( RCMIN, R( J ) )
 436                RCMAX = MAX( RCMAX, R( J ) )
 437    10       CONTINUE
 438             IF( RCMIN.LE.ZERO ) THEN
 439                INFO = -11
 440             ELSE IF( N.GT.0 ) THEN
 441                ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 442             ELSE
 443                ROWCND = ONE
 444             END IF
 445          END IF
 446          IF( COLEQU .AND. INFO.EQ.0 ) THEN
 447             RCMIN = BIGNUM
 448             RCMAX = ZERO
 449             DO 20 J = 1, N
 450                RCMIN = MIN( RCMIN, C( J ) )
 451                RCMAX = MAX( RCMAX, C( J ) )
 452    20       CONTINUE
 453             IF( RCMIN.LE.ZERO ) THEN
 454                INFO = -12
 455             ELSE IF( N.GT.0 ) THEN
 456                COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 457             ELSE
 458                COLCND = ONE
 459             END IF
 460          END IF
 461          IF( INFO.EQ.0 ) THEN
 462             IF( LDB.LT.MAX( 1, N ) ) THEN
 463                INFO = -14
 464             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 465                INFO = -16
 466             END IF
 467          END IF
 468       END IF
 469 *
 470       IF( INFO.NE.0 ) THEN
 471          CALL XERBLA( 'DGESVX', -INFO )
 472          RETURN
 473       END IF
 474 *
 475       IF( EQUIL ) THEN
 476 *
 477 *        Compute row and column scalings to equilibrate the matrix A.
 478 *
 479          CALL DGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU )
 480          IF( INFEQU.EQ.0 ) THEN
 481 *
 482 *           Equilibrate the matrix.
 483 *
 484             CALL DLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 485      $                   EQUED )
 486             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 487             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 488          END IF
 489       END IF
 490 *
 491 *     Scale the right hand side.
 492 *
 493       IF( NOTRAN ) THEN
 494          IF( ROWEQU ) THEN
 495             DO 40 J = 1, NRHS
 496                DO 30 I = 1, N
 497                   B( I, J ) = R( I )*B( I, J )
 498    30          CONTINUE
 499    40       CONTINUE
 500          END IF
 501       ELSE IF( COLEQU ) THEN
 502          DO 60 J = 1, NRHS
 503             DO 50 I = 1, N
 504                B( I, J ) = C( I )*B( I, J )
 505    50       CONTINUE
 506    60    CONTINUE
 507       END IF
 508 *
 509       IF( NOFACT .OR. EQUIL ) THEN
 510 *
 511 *        Compute the LU factorization of A.
 512 *
 513          CALL DLACPY( 'Full', N, N, A, LDA, AF, LDAF )
 514          CALL DGETRF( N, N, AF, LDAF, IPIV, INFO )
 515 *
 516 *        Return if INFO is non-zero.
 517 *
 518          IF( INFO.GT.0 ) THEN
 519 *
 520 *           Compute the reciprocal pivot growth factor of the
 521 *           leading rank-deficient INFO columns of A.
 522 *
 523             RPVGRW = DLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF,
 524      $               WORK )
 525             IF( RPVGRW.EQ.ZERO ) THEN
 526                RPVGRW = ONE
 527             ELSE
 528                RPVGRW = DLANGE( 'M', N, INFO, A, LDA, WORK ) / RPVGRW
 529             END IF
 530             WORK( 1 ) = RPVGRW
 531             RCOND = ZERO
 532             RETURN
 533          END IF
 534       END IF
 535 *
 536 *     Compute the norm of the matrix A and the
 537 *     reciprocal pivot growth factor RPVGRW.
 538 *
 539       IF( NOTRAN ) THEN
 540          NORM = '1'
 541       ELSE
 542          NORM = 'I'
 543       END IF
 544       ANORM = DLANGE( NORM, N, N, A, LDA, WORK )
 545       RPVGRW = DLANTR( 'M', 'U', 'N', N, N, AF, LDAF, WORK )
 546       IF( RPVGRW.EQ.ZERO ) THEN
 547          RPVGRW = ONE
 548       ELSE
 549          RPVGRW = DLANGE( 'M', N, N, A, LDA, WORK ) / RPVGRW
 550       END IF
 551 *
 552 *     Compute the reciprocal of the condition number of A.
 553 *
 554       CALL DGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO )
 555 *
 556 *     Compute the solution matrix X.
 557 *
 558       CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
 559       CALL DGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
 560 *
 561 *     Use iterative refinement to improve the computed solution and
 562 *     compute error bounds and backward error estimates for it.
 563 *
 564       CALL DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
 565      $             LDX, FERR, BERR, WORK, IWORK, INFO )
 566 *
 567 *     Transform the solution matrix X to a solution of the original
 568 *     system.
 569 *
 570       IF( NOTRAN ) THEN
 571          IF( COLEQU ) THEN
 572             DO 80 J = 1, NRHS
 573                DO 70 I = 1, N
 574                   X( I, J ) = C( I )*X( I, J )
 575    70          CONTINUE
 576    80       CONTINUE
 577             DO 90 J = 1, NRHS
 578                FERR( J ) = FERR( J ) / COLCND
 579    90       CONTINUE
 580          END IF
 581       ELSE IF( ROWEQU ) THEN
 582          DO 110 J = 1, NRHS
 583             DO 100 I = 1, N
 584                X( I, J ) = R( I )*X( I, J )
 585   100       CONTINUE
 586   110    CONTINUE
 587          DO 120 J = 1, NRHS
 588             FERR( J ) = FERR( J ) / ROWCND
 589   120    CONTINUE
 590       END IF
 591 *
 592       WORK( 1 ) = RPVGRW
 593 *
 594 *     Set INFO = N+1 if the matrix is singular to working precision.
 595 *
 596       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
 597      $   INFO = N + 1
 598       RETURN
 599 *
 600 *     End of DGESVX
 601 *
 602       END