Annotation of rpl/lapack/lapack/dgesvxx.f, revision 1.16

1.5       bertrand    1: *> \brief <b> DGESVXX computes the solution to system of linear equations A * X = B for GE matrices</b>
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.12      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.5       bertrand    7: *
                      8: *> \htmlonly
1.12      bertrand    9: *> Download DGESVXX + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvxx.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvxx.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvxx.f">
1.5       bertrand   15: *> [TXT]</a>
1.12      bertrand   16: *> \endhtmlonly
1.5       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGESVXX( 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, IWORK,
                     25: *                           INFO )
1.12      bertrand   26: *
1.5       bertrand   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( * ), IWORK( * )
                     35: *       DOUBLE PRECISION   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, * )
                     40: *       ..
1.12      bertrand   41: *
1.5       bertrand   42: *
                     43: *> \par Purpose:
                     44: *  =============
                     45: *>
                     46: *> \verbatim
                     47: *>
                     48: *>    DGESVXX uses the LU factorization to compute the solution to a
                     49: *>    double precision 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. DGESVXX 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: *>    DGESVXX 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: *>    DGESVXX 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 DGESVXX 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 = 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 DOUBLE PRECISION 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
1.7       bertrand  186: *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
1.5       bertrand  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 DGETRF.  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
1.7       bertrand  210: *>          IPIV is INTEGER array, dimension (N)
1.5       bertrand  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 DGETRF; 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
1.7       bertrand  227: *>          EQUED is CHARACTER*1
1.5       bertrand  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
1.7       bertrand  242: *>          R is DOUBLE PRECISION array, dimension (N)
1.5       bertrand  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
1.7       bertrand  260: *>          C is DOUBLE PRECISION array, dimension (N)
1.5       bertrand  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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DGESVX, 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
1.15      bertrand  414: *>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
1.5       bertrand  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
1.15      bertrand  450: *>     Specifies the number of parameters set in PARAMS.  If <= 0, the
1.5       bertrand  451: *>     PARAMS array is never referenced and default values are used.
                    452: *> \endverbatim
                    453: *>
                    454: *> \param[in,out] PARAMS
                    455: *> \verbatim
1.7       bertrand  456: *>          PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
1.15      bertrand  457: *>     Specifies algorithm parameters.  If an entry is < 0.0, then
1.5       bertrand  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
1.15      bertrand  465: *>            = 0.0:  No refinement is performed, and no error bounds are
1.5       bertrand  466: *>                    computed.
1.15      bertrand  467: *>            = 1.0:  Use the extra-precise refinement algorithm.
1.5       bertrand  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 DOUBLE PRECISION array, dimension (4*N)
                    490: *> \endverbatim
                    491: *>
                    492: *> \param[out] IWORK
                    493: *> \verbatim
                    494: *>          IWORK is INTEGER array, dimension (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: *
1.12      bertrand  527: *> \author Univ. of Tennessee
                    528: *> \author Univ. of California Berkeley
                    529: *> \author Univ. of Colorado Denver
                    530: *> \author NAG Ltd.
1.5       bertrand  531: *
                    532: *> \ingroup doubleGEsolve
                    533: *
                    534: *  =====================================================================
1.1       bertrand  535:       SUBROUTINE DGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
                    536:      $                    EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
                    537:      $                    BERR, N_ERR_BNDS, ERR_BNDS_NORM,
                    538:      $                    ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK,
                    539:      $                    INFO )
                    540: *
1.16    ! bertrand  541: *  -- LAPACK driver routine --
1.5       bertrand  542: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    543: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.1       bertrand  544: *
                    545: *     .. Scalar Arguments ..
                    546:       CHARACTER          EQUED, FACT, TRANS
                    547:       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
                    548:      $                   N_ERR_BNDS
                    549:       DOUBLE PRECISION   RCOND, RPVGRW
                    550: *     ..
                    551: *     .. Array Arguments ..
                    552:       INTEGER            IPIV( * ), IWORK( * )
                    553:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
                    554:      $                   X( LDX , * ),WORK( * )
                    555:       DOUBLE PRECISION   R( * ), C( * ), PARAMS( * ), BERR( * ),
                    556:      $                   ERR_BNDS_NORM( NRHS, * ),
                    557:      $                   ERR_BNDS_COMP( NRHS, * )
                    558: *     ..
                    559: *
1.5       bertrand  560: *  =====================================================================
1.1       bertrand  561: *
                    562: *     .. Parameters ..
                    563:       DOUBLE PRECISION   ZERO, ONE
                    564:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
                    565:       INTEGER            FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
                    566:       INTEGER            RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
                    567:       INTEGER            CMP_ERR_I, PIV_GROWTH_I
                    568:       PARAMETER          ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
                    569:      $                   BERR_I = 3 )
                    570:       PARAMETER          ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
                    571:       PARAMETER          ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
                    572:      $                   PIV_GROWTH_I = 9 )
                    573: *     ..
                    574: *     .. Local Scalars ..
                    575:       LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
                    576:       INTEGER            INFEQU, J
                    577:       DOUBLE PRECISION   AMAX, BIGNUM, COLCND, RCMAX, RCMIN, ROWCND,
                    578:      $                   SMLNUM
                    579: *     ..
                    580: *     .. External Functions ..
1.7       bertrand  581:       EXTERNAL           LSAME, DLAMCH, DLA_GERPVGRW
1.1       bertrand  582:       LOGICAL            LSAME
1.7       bertrand  583:       DOUBLE PRECISION   DLAMCH, DLA_GERPVGRW
1.1       bertrand  584: *     ..
                    585: *     .. External Subroutines ..
                    586:       EXTERNAL           DGEEQUB, DGETRF, DGETRS, DLACPY, DLAQGE,
                    587:      $                   XERBLA, DLASCL2, DGERFSX
                    588: *     ..
                    589: *     .. Intrinsic Functions ..
                    590:       INTRINSIC          MAX, MIN
                    591: *     ..
                    592: *     .. Executable Statements ..
                    593: *
                    594:       INFO = 0
                    595:       NOFACT = LSAME( FACT, 'N' )
                    596:       EQUIL = LSAME( FACT, 'E' )
                    597:       NOTRAN = LSAME( TRANS, 'N' )
                    598:       SMLNUM = DLAMCH( 'Safe minimum' )
                    599:       BIGNUM = ONE / SMLNUM
                    600:       IF( NOFACT .OR. EQUIL ) THEN
                    601:          EQUED = 'N'
                    602:          ROWEQU = .FALSE.
                    603:          COLEQU = .FALSE.
                    604:       ELSE
                    605:          ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
                    606:          COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
                    607:       END IF
                    608: *
                    609: *     Default is failure.  If an input parameter is wrong or
                    610: *     factorization fails, make everything look horrible.  Only the
                    611: *     pivot growth is set here, the rest is initialized in DGERFSX.
                    612: *
                    613:       RPVGRW = ZERO
                    614: *
                    615: *     Test the input parameters.  PARAMS is not tested until DGERFSX.
                    616: *
                    617:       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.
                    618:      $     LSAME( FACT, 'F' ) ) THEN
                    619:          INFO = -1
                    620:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
                    621:      $        LSAME( TRANS, 'C' ) ) THEN
                    622:          INFO = -2
                    623:       ELSE IF( N.LT.0 ) THEN
                    624:          INFO = -3
                    625:       ELSE IF( NRHS.LT.0 ) THEN
                    626:          INFO = -4
                    627:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    628:          INFO = -6
                    629:       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
                    630:          INFO = -8
                    631:       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
                    632:      $        ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
                    633:          INFO = -10
                    634:       ELSE
                    635:          IF( ROWEQU ) THEN
                    636:             RCMIN = BIGNUM
                    637:             RCMAX = ZERO
                    638:             DO 10 J = 1, N
                    639:                RCMIN = MIN( RCMIN, R( J ) )
                    640:                RCMAX = MAX( RCMAX, R( J ) )
                    641:  10         CONTINUE
                    642:             IF( RCMIN.LE.ZERO ) THEN
                    643:                INFO = -11
                    644:             ELSE IF( N.GT.0 ) THEN
                    645:                ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
                    646:             ELSE
                    647:                ROWCND = ONE
                    648:             END IF
                    649:          END IF
                    650:          IF( COLEQU .AND. INFO.EQ.0 ) THEN
                    651:             RCMIN = BIGNUM
                    652:             RCMAX = ZERO
                    653:             DO 20 J = 1, N
                    654:                RCMIN = MIN( RCMIN, C( J ) )
                    655:                RCMAX = MAX( RCMAX, C( J ) )
                    656:  20         CONTINUE
                    657:             IF( RCMIN.LE.ZERO ) THEN
                    658:                INFO = -12
                    659:             ELSE IF( N.GT.0 ) THEN
                    660:                COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
                    661:             ELSE
                    662:                COLCND = ONE
                    663:             END IF
                    664:          END IF
                    665:          IF( INFO.EQ.0 ) THEN
                    666:             IF( LDB.LT.MAX( 1, N ) ) THEN
                    667:                INFO = -14
                    668:             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
                    669:                INFO = -16
                    670:             END IF
                    671:          END IF
                    672:       END IF
                    673: *
                    674:       IF( INFO.NE.0 ) THEN
                    675:          CALL XERBLA( 'DGESVXX', -INFO )
                    676:          RETURN
                    677:       END IF
                    678: *
                    679:       IF( EQUIL ) THEN
                    680: *
                    681: *     Compute row and column scalings to equilibrate the matrix A.
                    682: *
                    683:          CALL DGEEQUB( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
                    684:      $        INFEQU )
                    685:          IF( INFEQU.EQ.0 ) THEN
                    686: *
                    687: *     Equilibrate the matrix.
                    688: *
                    689:             CALL DLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
                    690:      $           EQUED )
                    691:             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
                    692:             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
                    693:          END IF
                    694: *
                    695: *     If the scaling factors are not applied, set them to 1.0.
                    696: *
                    697:          IF ( .NOT.ROWEQU ) THEN
                    698:             DO J = 1, N
                    699:                R( J ) = 1.0D+0
                    700:             END DO
                    701:          END IF
                    702:          IF ( .NOT.COLEQU ) THEN
                    703:             DO J = 1, N
                    704:                C( J ) = 1.0D+0
                    705:             END DO
                    706:          END IF
                    707:       END IF
                    708: *
                    709: *     Scale the right-hand side.
                    710: *
                    711:       IF( NOTRAN ) THEN
                    712:          IF( ROWEQU ) CALL DLASCL2( N, NRHS, R, B, LDB )
                    713:       ELSE
                    714:          IF( COLEQU ) CALL DLASCL2( N, NRHS, C, B, LDB )
                    715:       END IF
                    716: *
                    717:       IF( NOFACT .OR. EQUIL ) THEN
                    718: *
                    719: *        Compute the LU factorization of A.
                    720: *
                    721:          CALL DLACPY( 'Full', N, N, A, LDA, AF, LDAF )
                    722:          CALL DGETRF( N, N, AF, LDAF, IPIV, INFO )
                    723: *
                    724: *        Return if INFO is non-zero.
                    725: *
                    726:          IF( INFO.GT.0 ) THEN
                    727: *
                    728: *           Pivot in column INFO is exactly 0
                    729: *           Compute the reciprocal pivot growth factor of the
                    730: *           leading rank-deficient INFO columns of A.
                    731: *
1.7       bertrand  732:             RPVGRW = DLA_GERPVGRW( N, INFO, A, LDA, AF, LDAF )
1.1       bertrand  733:             RETURN
                    734:          END IF
                    735:       END IF
                    736: *
                    737: *     Compute the reciprocal pivot growth factor RPVGRW.
                    738: *
1.7       bertrand  739:       RPVGRW = DLA_GERPVGRW( N, N, A, LDA, AF, LDAF )
1.1       bertrand  740: *
                    741: *     Compute the solution matrix X.
                    742: *
                    743:       CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
                    744:       CALL DGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
                    745: *
                    746: *     Use iterative refinement to improve the computed solution and
                    747: *     compute error bounds and backward error estimates for it.
                    748: *
                    749:       CALL DGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF,
                    750:      $     IPIV, R, C, B, LDB, X, LDX, RCOND, BERR,
                    751:      $     N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
                    752:      $     WORK, IWORK, INFO )
                    753: *
                    754: *     Scale solutions.
                    755: *
                    756:       IF ( COLEQU .AND. NOTRAN ) THEN
                    757:          CALL DLASCL2 ( N, NRHS, C, X, LDX )
                    758:       ELSE IF ( ROWEQU .AND. .NOT.NOTRAN ) THEN
                    759:          CALL DLASCL2 ( N, NRHS, R, X, LDX )
                    760:       END IF
                    761: *
                    762:       RETURN
                    763: *
                    764: *     End of DGESVXX
                    765: 
                    766:       END

CVSweb interface <joel.bertrand@systella.fr>