SRC/cgesvxx.f

   1 *> \brief <b> CGESVXX 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 CGESVXX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvxx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvxx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvxx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
  22 *                           EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
  23 *                           BERR, N_ERR_BNDS, ERR_BNDS_NORM,
  24 *                           ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
  25 *                           INFO )
  26 *
  27 *       .. Scalar Arguments ..
  28 *       CHARACTER          EQUED, FACT, TRANS
  29 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
  30 *      $                   N_ERR_BNDS
  31 *       REAL               RCOND, RPVGRW
  32 *       ..
  33 *       .. Array Arguments ..
  34 *       INTEGER            IPIV( * )
  35 *       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  36 *      $                   X( LDX , * ),WORK( * )
  37 *       REAL               R( * ), C( * ), PARAMS( * ), BERR( * ),
  38 *      $                   ERR_BNDS_NORM( NRHS, * ),
  39 *      $                   ERR_BNDS_COMP( NRHS, * ), RWORK( * )
  40 *       ..
  41 *
  42 *
  43 *> \par Purpose:
  44 *  =============
  45 *>
  46 *> \verbatim
  47 *>
  48 *>    CGESVXX uses the LU factorization to compute the solution to a
  49 *>    complex system of linear equations  A * X = B,  where A is an
  50 *>    N-by-N matrix and X and B are N-by-NRHS matrices.
  51 *>
  52 *>    If requested, both normwise and maximum componentwise error bounds
  53 *>    are returned. CGESVXX will return a solution with a tiny
  54 *>    guaranteed error (O(eps) where eps is the working machine
  55 *>    precision) unless the matrix is very ill-conditioned, in which
  56 *>    case a warning is returned. Relevant condition numbers also are
  57 *>    calculated and returned.
  58 *>
  59 *>    CGESVXX accepts user-provided factorizations and equilibration
  60 *>    factors; see the definitions of the FACT and EQUED options.
  61 *>    Solving with refinement and using a factorization from a previous
  62 *>    CGESVXX call will also produce a solution with either O(eps)
  63 *>    errors or warnings, but we cannot make that claim for general
  64 *>    user-provided factorizations and equilibration factors if they
  65 *>    differ from what CGESVXX would itself produce.
  66 *> \endverbatim
  67 *
  68 *> \par Description:
  69 *  =================
  70 *>
  71 *> \verbatim
  72 *>
  73 *>    The following steps are performed:
  74 *>
  75 *>    1. If FACT = 'E', real scaling factors are computed to equilibrate
  76 *>    the system:
  77 *>
  78 *>      TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
  79 *>      TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
  80 *>      TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
  81 *>
  82 *>    Whether or not the system will be equilibrated depends on the
  83 *>    scaling of the matrix A, but if equilibration is used, A is
  84 *>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
  85 *>    or diag(C)*B (if TRANS = 'T' or 'C').
  86 *>
  87 *>    2. If FACT = 'N' or 'E', the LU decomposition is used to factor
  88 *>    the matrix A (after equilibration if FACT = 'E') as
  89 *>
  90 *>      A = P * L * U,
  91 *>
  92 *>    where P is a permutation matrix, L is a unit lower triangular
  93 *>    matrix, and U is upper triangular.
  94 *>
  95 *>    3. If some U(i,i)=0, so that U is exactly singular, then the
  96 *>    routine returns with INFO = i. Otherwise, the factored form of A
  97 *>    is used to estimate the condition number of the matrix A (see
  98 *>    argument RCOND). If the reciprocal of the condition number is less
  99 *>    than machine precision, the routine still goes on to solve for X
 100 *>    and compute error bounds as described below.
 101 *>
 102 *>    4. The system of equations is solved for X using the factored form
 103 *>    of A.
 104 *>
 105 *>    5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
 106 *>    the routine will use iterative refinement to try to get a small
 107 *>    error and error bounds.  Refinement calculates the residual to at
 108 *>    least twice the working precision.
 109 *>
 110 *>    6. If equilibration was used, the matrix X is premultiplied by
 111 *>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
 112 *>    that it solves the original system before equilibration.
 113 *> \endverbatim
 114 *
 115 *  Arguments:
 116 *  ==========
 117 *
 118 *> \verbatim
 119 *>     Some optional parameters are bundled in the PARAMS array.  These
 120 *>     settings determine how refinement is performed, but often the
 121 *>     defaults are acceptable.  If the defaults are acceptable, users
 122 *>     can pass NPARAMS = 0 which prevents the source code from accessing
 123 *>     the PARAMS argument.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] FACT
 127 *> \verbatim
 128 *>          FACT is CHARACTER*1
 129 *>     Specifies whether or not the factored form of the matrix A is
 130 *>     supplied on entry, and if not, whether the matrix A should be
 131 *>     equilibrated before it is factored.
 132 *>       = 'F':  On entry, AF and IPIV contain the factored form of A.
 133 *>               If EQUED is not 'N', the matrix A has been
 134 *>               equilibrated with scaling factors given by R and C.
 135 *>               A, AF, and IPIV are not modified.
 136 *>       = 'N':  The matrix A will be copied to AF and factored.
 137 *>       = 'E':  The matrix A will be equilibrated if necessary, then
 138 *>               copied to AF and factored.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] TRANS
 142 *> \verbatim
 143 *>          TRANS is CHARACTER*1
 144 *>     Specifies the form of the system of equations:
 145 *>       = 'N':  A * X = B     (No transpose)
 146 *>       = 'T':  A**T * X = B  (Transpose)
 147 *>       = 'C':  A**H * X = B  (Conjugate Transpose)
 148 *> \endverbatim
 149 *>
 150 *> \param[in] N
 151 *> \verbatim
 152 *>          N is INTEGER
 153 *>     The number of linear equations, i.e., the order of the
 154 *>     matrix A.  N >= 0.
 155 *> \endverbatim
 156 *>
 157 *> \param[in] NRHS
 158 *> \verbatim
 159 *>          NRHS is INTEGER
 160 *>     The number of right hand sides, i.e., the number of columns
 161 *>     of the matrices B and X.  NRHS >= 0.
 162 *> \endverbatim
 163 *>
 164 *> \param[in,out] A
 165 *> \verbatim
 166 *>          A is COMPLEX array, dimension (LDA,N)
 167 *>     On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
 168 *>     not 'N', then A must have been equilibrated by the scaling
 169 *>     factors in R and/or C.  A is not modified if FACT = 'F' or
 170 *>     'N', or if FACT = 'E' and EQUED = 'N' on exit.
 171 *>
 172 *>     On exit, if EQUED .ne. 'N', A is scaled as follows:
 173 *>     EQUED = 'R':  A := diag(R) * A
 174 *>     EQUED = 'C':  A := A * diag(C)
 175 *>     EQUED = 'B':  A := diag(R) * A * diag(C).
 176 *> \endverbatim
 177 *>
 178 *> \param[in] LDA
 179 *> \verbatim
 180 *>          LDA is INTEGER
 181 *>     The leading dimension of the array A.  LDA >= max(1,N).
 182 *> \endverbatim
 183 *>
 184 *> \param[in,out] AF
 185 *> \verbatim
 186 *>          AF is COMPLEX array, dimension (LDAF,N)
 187 *>     If FACT = 'F', then AF is an input argument and on entry
 188 *>     contains the factors L and U from the factorization
 189 *>     A = P*L*U as computed by CGETRF.  If EQUED .ne. 'N', then
 190 *>     AF is the factored form of the equilibrated matrix A.
 191 *>
 192 *>     If FACT = 'N', then AF is an output argument and on exit
 193 *>     returns the factors L and U from the factorization A = P*L*U
 194 *>     of the original matrix A.
 195 *>
 196 *>     If FACT = 'E', then AF is an output argument and on exit
 197 *>     returns the factors L and U from the factorization A = P*L*U
 198 *>     of the equilibrated matrix A (see the description of A for
 199 *>     the form of the equilibrated matrix).
 200 *> \endverbatim
 201 *>
 202 *> \param[in] LDAF
 203 *> \verbatim
 204 *>          LDAF is INTEGER
 205 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
 206 *> \endverbatim
 207 *>
 208 *> \param[in,out] IPIV
 209 *> \verbatim
 210 *>          IPIV is INTEGER array, dimension (N)
 211 *>     If FACT = 'F', then IPIV is an input argument and on entry
 212 *>     contains the pivot indices from the factorization A = P*L*U
 213 *>     as computed by CGETRF; row i of the matrix was interchanged
 214 *>     with row IPIV(i).
 215 *>
 216 *>     If FACT = 'N', then IPIV is an output argument and on exit
 217 *>     contains the pivot indices from the factorization A = P*L*U
 218 *>     of the original matrix A.
 219 *>
 220 *>     If FACT = 'E', then IPIV is an output argument and on exit
 221 *>     contains the pivot indices from the factorization A = P*L*U
 222 *>     of the equilibrated matrix A.
 223 *> \endverbatim
 224 *>
 225 *> \param[in,out] EQUED
 226 *> \verbatim
 227 *>          EQUED is CHARACTER*1
 228 *>     Specifies the form of equilibration that was done.
 229 *>       = 'N':  No equilibration (always true if FACT = 'N').
 230 *>       = 'R':  Row equilibration, i.e., A has been premultiplied by
 231 *>               diag(R).
 232 *>       = 'C':  Column equilibration, i.e., A has been postmultiplied
 233 *>               by diag(C).
 234 *>       = 'B':  Both row and column equilibration, i.e., A has been
 235 *>               replaced by diag(R) * A * diag(C).
 236 *>     EQUED is an input argument if FACT = 'F'; otherwise, it is an
 237 *>     output argument.
 238 *> \endverbatim
 239 *>
 240 *> \param[in,out] R
 241 *> \verbatim
 242 *>          R is REAL array, dimension (N)
 243 *>     The row scale factors for A.  If EQUED = 'R' or 'B', A is
 244 *>     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
 245 *>     is not accessed.  R is an input argument if FACT = 'F';
 246 *>     otherwise, R is an output argument.  If FACT = 'F' and
 247 *>     EQUED = 'R' or 'B', each element of R must be positive.
 248 *>     If R is output, each element of R is a power of the radix.
 249 *>     If R is input, each element of R should be a power of the radix
 250 *>     to ensure a reliable solution and error estimates. Scaling by
 251 *>     powers of the radix does not cause rounding errors unless the
 252 *>     result underflows or overflows. Rounding errors during scaling
 253 *>     lead to refining with a matrix that is not equivalent to the
 254 *>     input matrix, producing error estimates that may not be
 255 *>     reliable.
 256 *> \endverbatim
 257 *>
 258 *> \param[in,out] C
 259 *> \verbatim
 260 *>          C is REAL array, dimension (N)
 261 *>     The column scale factors for A.  If EQUED = 'C' or 'B', A is
 262 *>     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
 263 *>     is not accessed.  C is an input argument if FACT = 'F';
 264 *>     otherwise, C is an output argument.  If FACT = 'F' and
 265 *>     EQUED = 'C' or 'B', each element of C must be positive.
 266 *>     If C is output, each element of C is a power of the radix.
 267 *>     If C is input, each element of C should be a power of the radix
 268 *>     to ensure a reliable solution and error estimates. Scaling by
 269 *>     powers of the radix does not cause rounding errors unless the
 270 *>     result underflows or overflows. Rounding errors during scaling
 271 *>     lead to refining with a matrix that is not equivalent to the
 272 *>     input matrix, producing error estimates that may not be
 273 *>     reliable.
 274 *> \endverbatim
 275 *>
 276 *> \param[in,out] B
 277 *> \verbatim
 278 *>          B is COMPLEX array, dimension (LDB,NRHS)
 279 *>     On entry, the N-by-NRHS right hand side matrix B.
 280 *>     On exit,
 281 *>     if EQUED = 'N', B is not modified;
 282 *>     if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
 283 *>        diag(R)*B;
 284 *>     if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
 285 *>        overwritten by diag(C)*B.
 286 *> \endverbatim
 287 *>
 288 *> \param[in] LDB
 289 *> \verbatim
 290 *>          LDB is INTEGER
 291 *>     The leading dimension of the array B.  LDB >= max(1,N).
 292 *> \endverbatim
 293 *>
 294 *> \param[out] X
 295 *> \verbatim
 296 *>          X is COMPLEX array, dimension (LDX,NRHS)
 297 *>     If INFO = 0, the N-by-NRHS solution matrix X to the original
 298 *>     system of equations.  Note that A and B are modified on exit
 299 *>     if EQUED .ne. 'N', and the solution to the equilibrated system is
 300 *>     inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
 301 *>     inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
 302 *> \endverbatim
 303 *>
 304 *> \param[in] LDX
 305 *> \verbatim
 306 *>          LDX is INTEGER
 307 *>     The leading dimension of the array X.  LDX >= max(1,N).
 308 *> \endverbatim
 309 *>
 310 *> \param[out] RCOND
 311 *> \verbatim
 312 *>          RCOND is REAL
 313 *>     Reciprocal scaled condition number.  This is an estimate of the
 314 *>     reciprocal Skeel condition number of the matrix A after
 315 *>     equilibration (if done).  If this is less than the machine
 316 *>     precision (in particular, if it is zero), the matrix is singular
 317 *>     to working precision.  Note that the error may still be small even
 318 *>     if this number is very small and the matrix appears ill-
 319 *>     conditioned.
 320 *> \endverbatim
 321 *>
 322 *> \param[out] RPVGRW
 323 *> \verbatim
 324 *>          RPVGRW is REAL
 325 *>     Reciprocal pivot growth.  On exit, this contains the reciprocal
 326 *>     pivot growth factor norm(A)/norm(U). The "max absolute element"
 327 *>     norm is used.  If this is much less than 1, then the stability of
 328 *>     the LU factorization of the (equilibrated) matrix A could be poor.
 329 *>     This also means that the solution X, estimated condition numbers,
 330 *>     and error bounds could be unreliable. If factorization fails with
 331 *>     0<INFO<=N, then this contains the reciprocal pivot growth factor
 332 *>     for the leading INFO columns of A.  In CGESVX, this quantity is
 333 *>     returned in WORK(1).
 334 *> \endverbatim
 335 *>
 336 *> \param[out] BERR
 337 *> \verbatim
 338 *>          BERR is REAL array, dimension (NRHS)
 339 *>     Componentwise relative backward error.  This is the
 340 *>     componentwise relative backward error of each solution vector X(j)
 341 *>     (i.e., the smallest relative change in any element of A or B that
 342 *>     makes X(j) an exact solution).
 343 *> \endverbatim
 344 *>
 345 *> \param[in] N_ERR_BNDS
 346 *> \verbatim
 347 *>          N_ERR_BNDS is INTEGER
 348 *>     Number of error bounds to return for each right hand side
 349 *>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
 350 *>     ERR_BNDS_COMP below.
 351 *> \endverbatim
 352 *>
 353 *> \param[out] ERR_BNDS_NORM
 354 *> \verbatim
 355 *>          ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
 356 *>     For each right-hand side, this array contains information about
 357 *>     various error bounds and condition numbers corresponding to the
 358 *>     normwise relative error, which is defined as follows:
 359 *>
 360 *>     Normwise relative error in the ith solution vector:
 361 *>             max_j (abs(XTRUE(j,i) - X(j,i)))
 362 *>            ------------------------------
 363 *>                  max_j abs(X(j,i))
 364 *>
 365 *>     The array is indexed by the type of error information as described
 366 *>     below. There currently are up to three pieces of information
 367 *>     returned.
 368 *>
 369 *>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
 370 *>     right-hand side.
 371 *>
 372 *>     The second index in ERR_BNDS_NORM(:,err) contains the following
 373 *>     three fields:
 374 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
 375 *>              reciprocal condition number is less than the threshold
 376 *>              sqrt(n) * slamch('Epsilon').
 377 *>
 378 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
 379 *>              almost certainly within a factor of 10 of the true error
 380 *>              so long as the next entry is greater than the threshold
 381 *>              sqrt(n) * slamch('Epsilon'). This error bound should only
 382 *>              be trusted if the previous boolean is true.
 383 *>
 384 *>     err = 3  Reciprocal condition number: Estimated normwise
 385 *>              reciprocal condition number.  Compared with the threshold
 386 *>              sqrt(n) * slamch('Epsilon') to determine if the error
 387 *>              estimate is "guaranteed". These reciprocal condition
 388 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
 389 *>              appropriately scaled matrix Z.
 390 *>              Let Z = S*A, where S scales each row by a power of the
 391 *>              radix so all absolute row sums of Z are approximately 1.
 392 *>
 393 *>     See Lapack Working Note 165 for further details and extra
 394 *>     cautions.
 395 *> \endverbatim
 396 *>
 397 *> \param[out] ERR_BNDS_COMP
 398 *> \verbatim
 399 *>          ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
 400 *>     For each right-hand side, this array contains information about
 401 *>     various error bounds and condition numbers corresponding to the
 402 *>     componentwise relative error, which is defined as follows:
 403 *>
 404 *>     Componentwise relative error in the ith solution vector:
 405 *>                    abs(XTRUE(j,i) - X(j,i))
 406 *>             max_j ----------------------
 407 *>                         abs(X(j,i))
 408 *>
 409 *>     The array is indexed by the right-hand side i (on which the
 410 *>     componentwise relative error depends), and the type of error
 411 *>     information as described below. There currently are up to three
 412 *>     pieces of information returned for each right-hand side. If
 413 *>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
 414 *>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
 415 *>     the first (:,N_ERR_BNDS) entries are returned.
 416 *>
 417 *>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
 418 *>     right-hand side.
 419 *>
 420 *>     The second index in ERR_BNDS_COMP(:,err) contains the following
 421 *>     three fields:
 422 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
 423 *>              reciprocal condition number is less than the threshold
 424 *>              sqrt(n) * slamch('Epsilon').
 425 *>
 426 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
 427 *>              almost certainly within a factor of 10 of the true error
 428 *>              so long as the next entry is greater than the threshold
 429 *>              sqrt(n) * slamch('Epsilon'). This error bound should only
 430 *>              be trusted if the previous boolean is true.
 431 *>
 432 *>     err = 3  Reciprocal condition number: Estimated componentwise
 433 *>              reciprocal condition number.  Compared with the threshold
 434 *>              sqrt(n) * slamch('Epsilon') to determine if the error
 435 *>              estimate is "guaranteed". These reciprocal condition
 436 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
 437 *>              appropriately scaled matrix Z.
 438 *>              Let Z = S*(A*diag(x)), where x is the solution for the
 439 *>              current right-hand side and S scales each row of
 440 *>              A*diag(x) by a power of the radix so all absolute row
 441 *>              sums of Z are approximately 1.
 442 *>
 443 *>     See Lapack Working Note 165 for further details and extra
 444 *>     cautions.
 445 *> \endverbatim
 446 *>
 447 *> \param[in] NPARAMS
 448 *> \verbatim
 449 *>          NPARAMS is INTEGER
 450 *>     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
 451 *>     PARAMS array is never referenced and default values are used.
 452 *> \endverbatim
 453 *>
 454 *> \param[in,out] PARAMS
 455 *> \verbatim
 456 *>          PARAMS is REAL array, dimension NPARAMS
 457 *>     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
 458 *>     that entry will be filled with default value used for that
 459 *>     parameter.  Only positions up to NPARAMS are accessed; defaults
 460 *>     are used for higher-numbered parameters.
 461 *>
 462 *>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
 463 *>            refinement or not.
 464 *>         Default: 1.0
 465 *>            = 0.0 : No refinement is performed, and no error bounds are
 466 *>                    computed.
 467 *>            = 1.0 : Use the double-precision refinement algorithm,
 468 *>                    possibly with doubled-single computations if the
 469 *>                    compilation environment does not support DOUBLE
 470 *>                    PRECISION.
 471 *>              (other values are reserved for future use)
 472 *>
 473 *>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
 474 *>            computations allowed for refinement.
 475 *>         Default: 10
 476 *>         Aggressive: Set to 100 to permit convergence using approximate
 477 *>                     factorizations or factorizations other than LU. If
 478 *>                     the factorization uses a technique other than
 479 *>                     Gaussian elimination, the guarantees in
 480 *>                     err_bnds_norm and err_bnds_comp may no longer be
 481 *>                     trustworthy.
 482 *>
 483 *>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
 484 *>            will attempt to find a solution with small componentwise
 485 *>            relative error in the double-precision algorithm.  Positive
 486 *>            is true, 0.0 is false.
 487 *>         Default: 1.0 (attempt componentwise convergence)
 488 *> \endverbatim
 489 *>
 490 *> \param[out] WORK
 491 *> \verbatim
 492 *>          WORK is COMPLEX array, dimension (2*N)
 493 *> \endverbatim
 494 *>
 495 *> \param[out] RWORK
 496 *> \verbatim
 497 *>          RWORK is REAL array, dimension (2*N)
 498 *> \endverbatim
 499 *>
 500 *> \param[out] INFO
 501 *> \verbatim
 502 *>          INFO is INTEGER
 503 *>       = 0:  Successful exit. The solution to every right-hand side is
 504 *>         guaranteed.
 505 *>       < 0:  If INFO = -i, the i-th argument had an illegal value
 506 *>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
 507 *>         has been completed, but the factor U is exactly singular, so
 508 *>         the solution and error bounds could not be computed. RCOND = 0
 509 *>         is returned.
 510 *>       = N+J: The solution corresponding to the Jth right-hand side is
 511 *>         not guaranteed. The solutions corresponding to other right-
 512 *>         hand sides K with K > J may not be guaranteed as well, but
 513 *>         only the first such right-hand side is reported. If a small
 514 *>         componentwise error is not requested (PARAMS(3) = 0.0) then
 515 *>         the Jth right-hand side is the first with a normwise error
 516 *>         bound that is not guaranteed (the smallest J such
 517 *>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
 518 *>         the Jth right-hand side is the first with either a normwise or
 519 *>         componentwise error bound that is not guaranteed (the smallest
 520 *>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
 521 *>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
 522 *>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
 523 *>         about all of the right-hand sides check ERR_BNDS_NORM or
 524 *>         ERR_BNDS_COMP.
 525 *> \endverbatim
 526 *
 527 *  Authors:
 528 *  ========
 529 *
 530 *> \author Univ. of Tennessee
 531 *> \author Univ. of California Berkeley
 532 *> \author Univ. of Colorado Denver
 533 *> \author NAG Ltd.
 534 *
 535 *> \date April 2012
 536 *
 537 *> \ingroup complexGEsolve
 538 *
 539 *  =====================================================================
 540       SUBROUTINE CGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
 541      $                    EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
 542      $                    BERR, N_ERR_BNDS, ERR_BNDS_NORM,
 543      $                    ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
 544      $                    INFO )
 545 *
 546 *  -- LAPACK driver routine (version 3.4.1) --
 547 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 548 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 549 *     April 2012
 550 *
 551 *     .. Scalar Arguments ..
 552       CHARACTER          EQUED, FACT, TRANS
 553       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
 554      $                   N_ERR_BNDS
 555       REAL               RCOND, RPVGRW
 556 *     ..
 557 *     .. Array Arguments ..
 558       INTEGER            IPIV( * )
 559       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
 560      $                   X( LDX , * ),WORK( * )
 561       REAL               R( * ), C( * ), PARAMS( * ), BERR( * ),
 562      $                   ERR_BNDS_NORM( NRHS, * ),
 563      $                   ERR_BNDS_COMP( NRHS, * ), RWORK( * )
 564 *     ..
 565 *
 566 *  ==================================================================
 567 *
 568 *     .. Parameters ..
 569       REAL               ZERO, ONE
 570       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 571       INTEGER            FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
 572       INTEGER            RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
 573       INTEGER            CMP_ERR_I, PIV_GROWTH_I
 574       PARAMETER          ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
 575      $                   BERR_I = 3 )
 576       PARAMETER          ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
 577       PARAMETER          ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
 578      $                   PIV_GROWTH_I = 9 )
 579 *     ..
 580 *     .. Local Scalars ..
 581       LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
 582       INTEGER            INFEQU, J
 583       REAL               AMAX, BIGNUM, COLCND, RCMAX, RCMIN,
 584      $                   ROWCND, SMLNUM
 585 *     ..
 586 *     .. External Functions ..
 587       EXTERNAL           LSAME, SLAMCH, CLA_GERPVGRW
 588       LOGICAL            LSAME
 589       REAL               SLAMCH, CLA_GERPVGRW
 590 *     ..
 591 *     .. External Subroutines ..
 592       EXTERNAL           CGEEQUB, CGETRF, CGETRS, CLACPY, CLAQGE,
 593      $                   XERBLA, CLASCL2, CGERFSX
 594 *     ..
 595 *     .. Intrinsic Functions ..
 596       INTRINSIC          MAX, MIN
 597 *     ..
 598 *     .. Executable Statements ..
 599 *
 600       INFO = 0
 601       NOFACT = LSAME( FACT, 'N' )
 602       EQUIL = LSAME( FACT, 'E' )
 603       NOTRAN = LSAME( TRANS, 'N' )
 604       SMLNUM = SLAMCH( 'Safe minimum' )
 605       BIGNUM = ONE / SMLNUM
 606       IF( NOFACT .OR. EQUIL ) THEN
 607          EQUED = 'N'
 608          ROWEQU = .FALSE.
 609          COLEQU = .FALSE.
 610       ELSE
 611          ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 612          COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 613       END IF
 614 *
 615 *     Default is failure.  If an input parameter is wrong or
 616 *     factorization fails, make everything look horrible.  Only the
 617 *     pivot growth is set here, the rest is initialized in CGERFSX.
 618 *
 619       RPVGRW = ZERO
 620 *
 621 *     Test the input parameters.  PARAMS is not tested until CGERFSX.
 622 *
 623       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.
 624      $     LSAME( FACT, 'F' ) ) THEN
 625          INFO = -1
 626       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 627      $        LSAME( TRANS, 'C' ) ) THEN
 628          INFO = -2
 629       ELSE IF( N.LT.0 ) THEN
 630          INFO = -3
 631       ELSE IF( NRHS.LT.0 ) THEN
 632          INFO = -4
 633       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 634          INFO = -6
 635       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 636          INFO = -8
 637       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
 638      $        ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
 639          INFO = -10
 640       ELSE
 641          IF( ROWEQU ) THEN
 642             RCMIN = BIGNUM
 643             RCMAX = ZERO
 644             DO 10 J = 1, N
 645                RCMIN = MIN( RCMIN, R( J ) )
 646                RCMAX = MAX( RCMAX, R( J ) )
 647  10         CONTINUE
 648             IF( RCMIN.LE.ZERO ) THEN
 649                INFO = -11
 650             ELSE IF( N.GT.0 ) THEN
 651                ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 652             ELSE
 653                ROWCND = ONE
 654             END IF
 655          END IF
 656          IF( COLEQU .AND. INFO.EQ.0 ) THEN
 657             RCMIN = BIGNUM
 658             RCMAX = ZERO
 659             DO 20 J = 1, N
 660                RCMIN = MIN( RCMIN, C( J ) )
 661                RCMAX = MAX( RCMAX, C( J ) )
 662  20         CONTINUE
 663             IF( RCMIN.LE.ZERO ) THEN
 664                INFO = -12
 665             ELSE IF( N.GT.0 ) THEN
 666                COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 667             ELSE
 668                COLCND = ONE
 669             END IF
 670          END IF
 671          IF( INFO.EQ.0 ) THEN
 672             IF( LDB.LT.MAX( 1, N ) ) THEN
 673                INFO = -14
 674             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 675                INFO = -16
 676             END IF
 677          END IF
 678       END IF
 679 *
 680       IF( INFO.NE.0 ) THEN
 681          CALL XERBLA( 'CGESVXX', -INFO )
 682          RETURN
 683       END IF
 684 *
 685       IF( EQUIL ) THEN
 686 *
 687 *     Compute row and column scalings to equilibrate the matrix A.
 688 *
 689          CALL CGEEQUB( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 690      $        INFEQU )
 691          IF( INFEQU.EQ.0 ) THEN
 692 *
 693 *     Equilibrate the matrix.
 694 *
 695             CALL CLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 696      $           EQUED )
 697             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 698             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 699          END IF
 700 *
 701 *     If the scaling factors are not applied, set them to 1.0.
 702 *
 703          IF ( .NOT.ROWEQU ) THEN
 704             DO J = 1, N
 705                R( J ) = 1.0
 706             END DO
 707          END IF
 708          IF ( .NOT.COLEQU ) THEN
 709             DO J = 1, N
 710                C( J ) = 1.0
 711             END DO
 712          END IF
 713       END IF
 714 *
 715 *     Scale the right-hand side.
 716 *
 717       IF( NOTRAN ) THEN
 718          IF( ROWEQU ) CALL CLASCL2( N, NRHS, R, B, LDB )
 719       ELSE
 720          IF( COLEQU ) CALL CLASCL2( N, NRHS, C, B, LDB )
 721       END IF
 722 *
 723       IF( NOFACT .OR. EQUIL ) THEN
 724 *
 725 *        Compute the LU factorization of A.
 726 *
 727          CALL CLACPY( 'Full', N, N, A, LDA, AF, LDAF )
 728          CALL CGETRF( N, N, AF, LDAF, IPIV, INFO )
 729 *
 730 *        Return if INFO is non-zero.
 731 *
 732          IF( INFO.GT.0 ) THEN
 733 *
 734 *           Pivot in column INFO is exactly 0
 735 *           Compute the reciprocal pivot growth factor of the
 736 *           leading rank-deficient INFO columns of A.
 737 *
 738             RPVGRW = CLA_GERPVGRW( N, INFO, A, LDA, AF, LDAF )
 739             RETURN
 740          END IF
 741       END IF
 742 *
 743 *     Compute the reciprocal pivot growth factor RPVGRW.
 744 *
 745       RPVGRW = CLA_GERPVGRW( N, N, A, LDA, AF, LDAF )
 746 *
 747 *     Compute the solution matrix X.
 748 *
 749       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
 750       CALL CGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
 751 *
 752 *     Use iterative refinement to improve the computed solution and
 753 *     compute error bounds and backward error estimates for it.
 754 *
 755       CALL CGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF,
 756      $     IPIV, R, C, B, LDB, X, LDX, RCOND, BERR,
 757      $     N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
 758      $     WORK, RWORK, INFO )
 759 *
 760 *     Scale solutions.
 761 *
 762       IF ( COLEQU .AND. NOTRAN ) THEN
 763          CALL CLASCL2 ( N, NRHS, C, X, LDX )
 764       ELSE IF ( ROWEQU .AND. .NOT.NOTRAN ) THEN
 765          CALL CLASCL2 ( N, NRHS, R, X, LDX )
 766       END IF
 767 *
 768       RETURN
 769 *
 770 *     End of CGESVXX
 771 *
 772       END