ENH: Improving the travis dashboard name
[platform/upstream/lapack.git] / 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