SRC/zgesvxx.f

   1 *> \brief <b> ZGESVXX 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 ZGESVXX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvxx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvxx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvxx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGESVXX( 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 *       DOUBLE PRECISION   RCOND, RPVGRW
  32 *       ..
  33 *       .. Array Arguments ..
  34 *       INTEGER            IPIV( * )
  35 *       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  36 *      $                   X( LDX , * ),WORK( * )
  37 *       DOUBLE PRECISION   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 *>    ZGESVXX uses the LU factorization to compute the solution to a
  49 *>    complex*16 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. ZGESVXX 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 *>    ZGESVXX 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 *>    ZGESVXX 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 ZGESVXX 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', double precision 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*16 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*16 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 ZGETRF.  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 ZGETRF; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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 DOUBLE PRECISION
 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 DOUBLE PRECISION
 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 ZGESVX, this quantity is
 333 *>     returned in WORK(1).
 334 *> \endverbatim
 335 *>
 336 *> \param[out] BERR
 337 *> \verbatim
 338 *>          BERR is DOUBLE PRECISION 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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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.0D+0
 465 *>            = 0.0 : No refinement is performed, and no error bounds are
 466 *>                    computed.
 467 *>            = 1.0 : Use the extra-precise refinement algorithm.
 468 *>              (other values are reserved for future use)
 469 *>
 470 *>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
 471 *>            computations allowed for refinement.
 472 *>         Default: 10
 473 *>         Aggressive: Set to 100 to permit convergence using approximate
 474 *>                     factorizations or factorizations other than LU. If
 475 *>                     the factorization uses a technique other than
 476 *>                     Gaussian elimination, the guarantees in
 477 *>                     err_bnds_norm and err_bnds_comp may no longer be
 478 *>                     trustworthy.
 479 *>
 480 *>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
 481 *>            will attempt to find a solution with small componentwise
 482 *>            relative error in the double-precision algorithm.  Positive
 483 *>            is true, 0.0 is false.
 484 *>         Default: 1.0 (attempt componentwise convergence)
 485 *> \endverbatim
 486 *>
 487 *> \param[out] WORK
 488 *> \verbatim
 489 *>          WORK is COMPLEX*16 array, dimension (2*N)
 490 *> \endverbatim
 491 *>
 492 *> \param[out] RWORK
 493 *> \verbatim
 494 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
 495 *> \endverbatim
 496 *>
 497 *> \param[out] INFO
 498 *> \verbatim
 499 *>          INFO is INTEGER
 500 *>       = 0:  Successful exit. The solution to every right-hand side is
 501 *>         guaranteed.
 502 *>       < 0:  If INFO = -i, the i-th argument had an illegal value
 503 *>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
 504 *>         has been completed, but the factor U is exactly singular, so
 505 *>         the solution and error bounds could not be computed. RCOND = 0
 506 *>         is returned.
 507 *>       = N+J: The solution corresponding to the Jth right-hand side is
 508 *>         not guaranteed. The solutions corresponding to other right-
 509 *>         hand sides K with K > J may not be guaranteed as well, but
 510 *>         only the first such right-hand side is reported. If a small
 511 *>         componentwise error is not requested (PARAMS(3) = 0.0) then
 512 *>         the Jth right-hand side is the first with a normwise error
 513 *>         bound that is not guaranteed (the smallest J such
 514 *>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
 515 *>         the Jth right-hand side is the first with either a normwise or
 516 *>         componentwise error bound that is not guaranteed (the smallest
 517 *>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
 518 *>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
 519 *>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
 520 *>         about all of the right-hand sides check ERR_BNDS_NORM or
 521 *>         ERR_BNDS_COMP.
 522 *> \endverbatim
 523 *
 524 *  Authors:
 525 *  ========
 526 *
 527 *> \author Univ. of Tennessee
 528 *> \author Univ. of California Berkeley
 529 *> \author Univ. of Colorado Denver
 530 *> \author NAG Ltd.
 531 *
 532 *> \date April 2012
 533 *
 534 *> \ingroup complex16GEsolve
 535 *
 536 *  =====================================================================
 537       SUBROUTINE ZGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
 538      $                    EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
 539      $                    BERR, N_ERR_BNDS, ERR_BNDS_NORM,
 540      $                    ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
 541      $                    INFO )
 542 *
 543 *  -- LAPACK driver routine (version 3.4.1) --
 544 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 545 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 546 *     April 2012
 547 *
 548 *     .. Scalar Arguments ..
 549       CHARACTER          EQUED, FACT, TRANS
 550       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
 551      $                   N_ERR_BNDS
 552       DOUBLE PRECISION   RCOND, RPVGRW
 553 *     ..
 554 *     .. Array Arguments ..
 555       INTEGER            IPIV( * )
 556       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
 557      $                   X( LDX , * ),WORK( * )
 558       DOUBLE PRECISION   R( * ), C( * ), PARAMS( * ), BERR( * ),
 559      $                   ERR_BNDS_NORM( NRHS, * ),
 560      $                   ERR_BNDS_COMP( NRHS, * ), RWORK( * )
 561 *     ..
 562 *
 563 *  ==================================================================
 564 *
 565 *     .. Parameters ..
 566       DOUBLE PRECISION   ZERO, ONE
 567       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 568       INTEGER            FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
 569       INTEGER            RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
 570       INTEGER            CMP_ERR_I, PIV_GROWTH_I
 571       PARAMETER          ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
 572      $                   BERR_I = 3 )
 573       PARAMETER          ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
 574       PARAMETER          ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
 575      $                   PIV_GROWTH_I = 9 )
 576 *     ..
 577 *     .. Local Scalars ..
 578       LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
 579       INTEGER            INFEQU, J
 580       DOUBLE PRECISION   AMAX, BIGNUM, COLCND, RCMAX, RCMIN,
 581      $                   ROWCND, SMLNUM
 582 *     ..
 583 *     .. External Functions ..
 584       EXTERNAL           LSAME, DLAMCH, ZLA_GERPVGRW
 585       LOGICAL            LSAME
 586       DOUBLE PRECISION   DLAMCH, ZLA_GERPVGRW
 587 *     ..
 588 *     .. External Subroutines ..
 589       EXTERNAL           ZGEEQUB, ZGETRF, ZGETRS, ZLACPY, ZLAQGE,
 590      $                   XERBLA, ZLASCL2, ZGERFSX
 591 *     ..
 592 *     .. Intrinsic Functions ..
 593       INTRINSIC          MAX, MIN
 594 *     ..
 595 *     .. Executable Statements ..
 596 *
 597       INFO = 0
 598       NOFACT = LSAME( FACT, 'N' )
 599       EQUIL = LSAME( FACT, 'E' )
 600       NOTRAN = LSAME( TRANS, 'N' )
 601       SMLNUM = DLAMCH( 'Safe minimum' )
 602       BIGNUM = ONE / SMLNUM
 603       IF( NOFACT .OR. EQUIL ) THEN
 604          EQUED = 'N'
 605          ROWEQU = .FALSE.
 606          COLEQU = .FALSE.
 607       ELSE
 608          ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 609          COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 610       END IF
 611 *
 612 *     Default is failure.  If an input parameter is wrong or
 613 *     factorization fails, make everything look horrible.  Only the
 614 *     pivot growth is set here, the rest is initialized in ZGERFSX.
 615 *
 616       RPVGRW = ZERO
 617 *
 618 *     Test the input parameters.  PARAMS is not tested until ZGERFSX.
 619 *
 620       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.
 621      $     LSAME( FACT, 'F' ) ) THEN
 622          INFO = -1
 623       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 624      $        LSAME( TRANS, 'C' ) ) THEN
 625          INFO = -2
 626       ELSE IF( N.LT.0 ) THEN
 627          INFO = -3
 628       ELSE IF( NRHS.LT.0 ) THEN
 629          INFO = -4
 630       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 631          INFO = -6
 632       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 633          INFO = -8
 634       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
 635      $        ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
 636          INFO = -10
 637       ELSE
 638          IF( ROWEQU ) THEN
 639             RCMIN = BIGNUM
 640             RCMAX = ZERO
 641             DO 10 J = 1, N
 642                RCMIN = MIN( RCMIN, R( J ) )
 643                RCMAX = MAX( RCMAX, R( J ) )
 644  10         CONTINUE
 645             IF( RCMIN.LE.ZERO ) THEN
 646                INFO = -11
 647             ELSE IF( N.GT.0 ) THEN
 648                ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 649             ELSE
 650                ROWCND = ONE
 651             END IF
 652          END IF
 653          IF( COLEQU .AND. INFO.EQ.0 ) THEN
 654             RCMIN = BIGNUM
 655             RCMAX = ZERO
 656             DO 20 J = 1, N
 657                RCMIN = MIN( RCMIN, C( J ) )
 658                RCMAX = MAX( RCMAX, C( J ) )
 659  20         CONTINUE
 660             IF( RCMIN.LE.ZERO ) THEN
 661                INFO = -12
 662             ELSE IF( N.GT.0 ) THEN
 663                COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 664             ELSE
 665                COLCND = ONE
 666             END IF
 667          END IF
 668          IF( INFO.EQ.0 ) THEN
 669             IF( LDB.LT.MAX( 1, N ) ) THEN
 670                INFO = -14
 671             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 672                INFO = -16
 673             END IF
 674          END IF
 675       END IF
 676 *
 677       IF( INFO.NE.0 ) THEN
 678          CALL XERBLA( 'ZGESVXX', -INFO )
 679          RETURN
 680       END IF
 681 *
 682       IF( EQUIL ) THEN
 683 *
 684 *     Compute row and column scalings to equilibrate the matrix A.
 685 *
 686          CALL ZGEEQUB( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 687      $        INFEQU )
 688          IF( INFEQU.EQ.0 ) THEN
 689 *
 690 *     Equilibrate the matrix.
 691 *
 692             CALL ZLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 693      $           EQUED )
 694             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
 695             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
 696          END IF
 697 *
 698 *     If the scaling factors are not applied, set them to 1.0.
 699 *
 700          IF ( .NOT.ROWEQU ) THEN
 701             DO J = 1, N
 702                R( J ) = 1.0D+0
 703             END DO
 704          END IF
 705          IF ( .NOT.COLEQU ) THEN
 706             DO J = 1, N
 707                C( J ) = 1.0D+0
 708             END DO
 709          END IF
 710       END IF
 711 *
 712 *     Scale the right-hand side.
 713 *
 714       IF( NOTRAN ) THEN
 715          IF( ROWEQU ) CALL ZLASCL2( N, NRHS, R, B, LDB )
 716       ELSE
 717          IF( COLEQU ) CALL ZLASCL2( N, NRHS, C, B, LDB )
 718       END IF
 719 *
 720       IF( NOFACT .OR. EQUIL ) THEN
 721 *
 722 *        Compute the LU factorization of A.
 723 *
 724          CALL ZLACPY( 'Full', N, N, A, LDA, AF, LDAF )
 725          CALL ZGETRF( N, N, AF, LDAF, IPIV, INFO )
 726 *
 727 *        Return if INFO is non-zero.
 728 *
 729          IF( INFO.GT.0 ) THEN
 730 *
 731 *           Pivot in column INFO is exactly 0
 732 *           Compute the reciprocal pivot growth factor of the
 733 *           leading rank-deficient INFO columns of A.
 734 *
 735             RPVGRW = ZLA_GERPVGRW( N, INFO, A, LDA, AF, LDAF )
 736             RETURN
 737          END IF
 738       END IF
 739 *
 740 *     Compute the reciprocal pivot growth factor RPVGRW.
 741 *
 742       RPVGRW = ZLA_GERPVGRW( N, N, A, LDA, AF, LDAF )
 743 *
 744 *     Compute the solution matrix X.
 745 *
 746       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
 747       CALL ZGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
 748 *
 749 *     Use iterative refinement to improve the computed solution and
 750 *     compute error bounds and backward error estimates for it.
 751 *
 752       CALL ZGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF,
 753      $     IPIV, R, C, B, LDB, X, LDX, RCOND, BERR,
 754      $     N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
 755      $     WORK, RWORK, INFO )
 756 *
 757 *     Scale solutions.
 758 *
 759       IF ( COLEQU .AND. NOTRAN ) THEN
 760          CALL ZLASCL2 ( N, NRHS, C, X, LDX )
 761       ELSE IF ( ROWEQU .AND. .NOT.NOTRAN ) THEN
 762          CALL ZLASCL2 ( N, NRHS, R, X, LDX )
 763       END IF
 764 *
 765       RETURN
 766 *
 767 *     End of ZGESVXX
 768 *
 769       END