Annotation of rpl/lapack/lapack/zgbrfsx.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGBRFSX( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB,
! 2: $ LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND,
! 3: $ BERR, N_ERR_BNDS, ERR_BNDS_NORM,
! 4: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
! 5: $ INFO )
! 6: *
! 7: * -- LAPACK routine (version 3.2.2) --
! 8: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
! 9: * -- Jason Riedy of Univ. of California Berkeley. --
! 10: * -- June 2010 --
! 11: *
! 12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 13: * -- Univ. of California Berkeley and NAG Ltd. --
! 14: *
! 15: IMPLICIT NONE
! 16: * ..
! 17: * .. Scalar Arguments ..
! 18: CHARACTER TRANS, EQUED
! 19: INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
! 20: $ NPARAMS, N_ERR_BNDS
! 21: DOUBLE PRECISION RCOND
! 22: * ..
! 23: * .. Array Arguments ..
! 24: INTEGER IPIV( * )
! 25: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
! 26: $ X( LDX , * ),WORK( * )
! 27: DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
! 28: $ ERR_BNDS_NORM( NRHS, * ),
! 29: $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
! 30: * ..
! 31: *
! 32: * Purpose
! 33: * =======
! 34: *
! 35: * ZGBRFSX improves the computed solution to a system of linear
! 36: * equations and provides error bounds and backward error estimates
! 37: * for the solution. In addition to normwise error bound, the code
! 38: * provides maximum componentwise error bound if possible. See
! 39: * comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
! 40: * error bounds.
! 41: *
! 42: * The original system of linear equations may have been equilibrated
! 43: * before calling this routine, as described by arguments EQUED, R
! 44: * and C below. In this case, the solution and error bounds returned
! 45: * are for the original unequilibrated system.
! 46: *
! 47: * Arguments
! 48: * =========
! 49: *
! 50: * Some optional parameters are bundled in the PARAMS array. These
! 51: * settings determine how refinement is performed, but often the
! 52: * defaults are acceptable. If the defaults are acceptable, users
! 53: * can pass NPARAMS = 0 which prevents the source code from accessing
! 54: * the PARAMS argument.
! 55: *
! 56: * TRANS (input) CHARACTER*1
! 57: * Specifies the form of the system of equations:
! 58: * = 'N': A * X = B (No transpose)
! 59: * = 'T': A**T * X = B (Transpose)
! 60: * = 'C': A**H * X = B (Conjugate transpose = Transpose)
! 61: *
! 62: * EQUED (input) CHARACTER*1
! 63: * Specifies the form of equilibration that was done to A
! 64: * before calling this routine. This is needed to compute
! 65: * the solution and error bounds correctly.
! 66: * = 'N': No equilibration
! 67: * = 'R': Row equilibration, i.e., A has been premultiplied by
! 68: * diag(R).
! 69: * = 'C': Column equilibration, i.e., A has been postmultiplied
! 70: * by diag(C).
! 71: * = 'B': Both row and column equilibration, i.e., A has been
! 72: * replaced by diag(R) * A * diag(C).
! 73: * The right hand side B has been changed accordingly.
! 74: *
! 75: * N (input) INTEGER
! 76: * The order of the matrix A. N >= 0.
! 77: *
! 78: * KL (input) INTEGER
! 79: * The number of subdiagonals within the band of A. KL >= 0.
! 80: *
! 81: * KU (input) INTEGER
! 82: * The number of superdiagonals within the band of A. KU >= 0.
! 83: *
! 84: * NRHS (input) INTEGER
! 85: * The number of right hand sides, i.e., the number of columns
! 86: * of the matrices B and X. NRHS >= 0.
! 87: *
! 88: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
! 89: * The original band matrix A, stored in rows 1 to KL+KU+1.
! 90: * The j-th column of A is stored in the j-th column of the
! 91: * array AB as follows:
! 92: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
! 93: *
! 94: * LDAB (input) INTEGER
! 95: * The leading dimension of the array AB. LDAB >= KL+KU+1.
! 96: *
! 97: * AFB (input) DOUBLE PRECISION array, dimension (LDAFB,N)
! 98: * Details of the LU factorization of the band matrix A, as
! 99: * computed by DGBTRF. U is stored as an upper triangular band
! 100: * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
! 101: * the multipliers used during the factorization are stored in
! 102: * rows KL+KU+2 to 2*KL+KU+1.
! 103: *
! 104: * LDAFB (input) INTEGER
! 105: * The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
! 106: *
! 107: * IPIV (input) INTEGER array, dimension (N)
! 108: * The pivot indices from DGETRF; for 1<=i<=N, row i of the
! 109: * matrix was interchanged with row IPIV(i).
! 110: *
! 111: * R (input or output) DOUBLE PRECISION array, dimension (N)
! 112: * The row scale factors for A. If EQUED = 'R' or 'B', A is
! 113: * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
! 114: * is not accessed. R is an input argument if FACT = 'F';
! 115: * otherwise, R is an output argument. If FACT = 'F' and
! 116: * EQUED = 'R' or 'B', each element of R must be positive.
! 117: * If R is output, each element of R is a power of the radix.
! 118: * If R is input, each element of R should be a power of the radix
! 119: * to ensure a reliable solution and error estimates. Scaling by
! 120: * powers of the radix does not cause rounding errors unless the
! 121: * result underflows or overflows. Rounding errors during scaling
! 122: * lead to refining with a matrix that is not equivalent to the
! 123: * input matrix, producing error estimates that may not be
! 124: * reliable.
! 125: *
! 126: * C (input or output) DOUBLE PRECISION array, dimension (N)
! 127: * The column scale factors for A. If EQUED = 'C' or 'B', A is
! 128: * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
! 129: * is not accessed. C is an input argument if FACT = 'F';
! 130: * otherwise, C is an output argument. If FACT = 'F' and
! 131: * EQUED = 'C' or 'B', each element of C must be positive.
! 132: * If C is output, each element of C is a power of the radix.
! 133: * If C is input, each element of C should be a power of the radix
! 134: * to ensure a reliable solution and error estimates. Scaling by
! 135: * powers of the radix does not cause rounding errors unless the
! 136: * result underflows or overflows. Rounding errors during scaling
! 137: * lead to refining with a matrix that is not equivalent to the
! 138: * input matrix, producing error estimates that may not be
! 139: * reliable.
! 140: *
! 141: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
! 142: * The right hand side matrix B.
! 143: *
! 144: * LDB (input) INTEGER
! 145: * The leading dimension of the array B. LDB >= max(1,N).
! 146: *
! 147: * X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
! 148: * On entry, the solution matrix X, as computed by DGETRS.
! 149: * On exit, the improved solution matrix X.
! 150: *
! 151: * LDX (input) INTEGER
! 152: * The leading dimension of the array X. LDX >= max(1,N).
! 153: *
! 154: * RCOND (output) DOUBLE PRECISION
! 155: * Reciprocal scaled condition number. This is an estimate of the
! 156: * reciprocal Skeel condition number of the matrix A after
! 157: * equilibration (if done). If this is less than the machine
! 158: * precision (in particular, if it is zero), the matrix is singular
! 159: * to working precision. Note that the error may still be small even
! 160: * if this number is very small and the matrix appears ill-
! 161: * conditioned.
! 162: *
! 163: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
! 164: * Componentwise relative backward error. This is the
! 165: * componentwise relative backward error of each solution vector X(j)
! 166: * (i.e., the smallest relative change in any element of A or B that
! 167: * makes X(j) an exact solution).
! 168: *
! 169: * N_ERR_BNDS (input) INTEGER
! 170: * Number of error bounds to return for each right hand side
! 171: * and each type (normwise or componentwise). See ERR_BNDS_NORM and
! 172: * ERR_BNDS_COMP below.
! 173: *
! 174: * ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
! 175: * For each right-hand side, this array contains information about
! 176: * various error bounds and condition numbers corresponding to the
! 177: * normwise relative error, which is defined as follows:
! 178: *
! 179: * Normwise relative error in the ith solution vector:
! 180: * max_j (abs(XTRUE(j,i) - X(j,i)))
! 181: * ------------------------------
! 182: * max_j abs(X(j,i))
! 183: *
! 184: * The array is indexed by the type of error information as described
! 185: * below. There currently are up to three pieces of information
! 186: * returned.
! 187: *
! 188: * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
! 189: * right-hand side.
! 190: *
! 191: * The second index in ERR_BNDS_NORM(:,err) contains the following
! 192: * three fields:
! 193: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
! 194: * reciprocal condition number is less than the threshold
! 195: * sqrt(n) * dlamch('Epsilon').
! 196: *
! 197: * err = 2 "Guaranteed" error bound: The estimated forward error,
! 198: * almost certainly within a factor of 10 of the true error
! 199: * so long as the next entry is greater than the threshold
! 200: * sqrt(n) * dlamch('Epsilon'). This error bound should only
! 201: * be trusted if the previous boolean is true.
! 202: *
! 203: * err = 3 Reciprocal condition number: Estimated normwise
! 204: * reciprocal condition number. Compared with the threshold
! 205: * sqrt(n) * dlamch('Epsilon') to determine if the error
! 206: * estimate is "guaranteed". These reciprocal condition
! 207: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
! 208: * appropriately scaled matrix Z.
! 209: * Let Z = S*A, where S scales each row by a power of the
! 210: * radix so all absolute row sums of Z are approximately 1.
! 211: *
! 212: * See Lapack Working Note 165 for further details and extra
! 213: * cautions.
! 214: *
! 215: * ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
! 216: * For each right-hand side, this array contains information about
! 217: * various error bounds and condition numbers corresponding to the
! 218: * componentwise relative error, which is defined as follows:
! 219: *
! 220: * Componentwise relative error in the ith solution vector:
! 221: * abs(XTRUE(j,i) - X(j,i))
! 222: * max_j ----------------------
! 223: * abs(X(j,i))
! 224: *
! 225: * The array is indexed by the right-hand side i (on which the
! 226: * componentwise relative error depends), and the type of error
! 227: * information as described below. There currently are up to three
! 228: * pieces of information returned for each right-hand side. If
! 229: * componentwise accuracy is not requested (PARAMS(3) = 0.0), then
! 230: * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
! 231: * the first (:,N_ERR_BNDS) entries are returned.
! 232: *
! 233: * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
! 234: * right-hand side.
! 235: *
! 236: * The second index in ERR_BNDS_COMP(:,err) contains the following
! 237: * three fields:
! 238: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
! 239: * reciprocal condition number is less than the threshold
! 240: * sqrt(n) * dlamch('Epsilon').
! 241: *
! 242: * err = 2 "Guaranteed" error bound: The estimated forward error,
! 243: * almost certainly within a factor of 10 of the true error
! 244: * so long as the next entry is greater than the threshold
! 245: * sqrt(n) * dlamch('Epsilon'). This error bound should only
! 246: * be trusted if the previous boolean is true.
! 247: *
! 248: * err = 3 Reciprocal condition number: Estimated componentwise
! 249: * reciprocal condition number. Compared with the threshold
! 250: * sqrt(n) * dlamch('Epsilon') to determine if the error
! 251: * estimate is "guaranteed". These reciprocal condition
! 252: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
! 253: * appropriately scaled matrix Z.
! 254: * Let Z = S*(A*diag(x)), where x is the solution for the
! 255: * current right-hand side and S scales each row of
! 256: * A*diag(x) by a power of the radix so all absolute row
! 257: * sums of Z are approximately 1.
! 258: *
! 259: * See Lapack Working Note 165 for further details and extra
! 260: * cautions.
! 261: *
! 262: * NPARAMS (input) INTEGER
! 263: * Specifies the number of parameters set in PARAMS. If .LE. 0, the
! 264: * PARAMS array is never referenced and default values are used.
! 265: *
! 266: * PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
! 267: * Specifies algorithm parameters. If an entry is .LT. 0.0, then
! 268: * that entry will be filled with default value used for that
! 269: * parameter. Only positions up to NPARAMS are accessed; defaults
! 270: * are used for higher-numbered parameters.
! 271: *
! 272: * PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
! 273: * refinement or not.
! 274: * Default: 1.0D+0
! 275: * = 0.0 : No refinement is performed, and no error bounds are
! 276: * computed.
! 277: * = 1.0 : Use the double-precision refinement algorithm,
! 278: * possibly with doubled-single computations if the
! 279: * compilation environment does not support DOUBLE
! 280: * PRECISION.
! 281: * (other values are reserved for future use)
! 282: *
! 283: * PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
! 284: * computations allowed for refinement.
! 285: * Default: 10
! 286: * Aggressive: Set to 100 to permit convergence using approximate
! 287: * factorizations or factorizations other than LU. If
! 288: * the factorization uses a technique other than
! 289: * Gaussian elimination, the guarantees in
! 290: * err_bnds_norm and err_bnds_comp may no longer be
! 291: * trustworthy.
! 292: *
! 293: * PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
! 294: * will attempt to find a solution with small componentwise
! 295: * relative error in the double-precision algorithm. Positive
! 296: * is true, 0.0 is false.
! 297: * Default: 1.0 (attempt componentwise convergence)
! 298: *
! 299: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
! 300: *
! 301: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 302: *
! 303: * INFO (output) INTEGER
! 304: * = 0: Successful exit. The solution to every right-hand side is
! 305: * guaranteed.
! 306: * < 0: If INFO = -i, the i-th argument had an illegal value
! 307: * > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
! 308: * has been completed, but the factor U is exactly singular, so
! 309: * the solution and error bounds could not be computed. RCOND = 0
! 310: * is returned.
! 311: * = N+J: The solution corresponding to the Jth right-hand side is
! 312: * not guaranteed. The solutions corresponding to other right-
! 313: * hand sides K with K > J may not be guaranteed as well, but
! 314: * only the first such right-hand side is reported. If a small
! 315: * componentwise error is not requested (PARAMS(3) = 0.0) then
! 316: * the Jth right-hand side is the first with a normwise error
! 317: * bound that is not guaranteed (the smallest J such
! 318: * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
! 319: * the Jth right-hand side is the first with either a normwise or
! 320: * componentwise error bound that is not guaranteed (the smallest
! 321: * J such that either ERR_BNDS_NORM(J,1) = 0.0 or
! 322: * ERR_BNDS_COMP(J,1) = 0.0). See the definition of
! 323: * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
! 324: * about all of the right-hand sides check ERR_BNDS_NORM or
! 325: * ERR_BNDS_COMP.
! 326: *
! 327: * ==================================================================
! 328: *
! 329: * .. Parameters ..
! 330: DOUBLE PRECISION ZERO, ONE
! 331: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 332: DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
! 333: DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
! 334: DOUBLE PRECISION DZTHRESH_DEFAULT
! 335: PARAMETER ( ITREF_DEFAULT = 1.0D+0 )
! 336: PARAMETER ( ITHRESH_DEFAULT = 10.0D+0 )
! 337: PARAMETER ( COMPONENTWISE_DEFAULT = 1.0D+0 )
! 338: PARAMETER ( RTHRESH_DEFAULT = 0.5D+0 )
! 339: PARAMETER ( DZTHRESH_DEFAULT = 0.25D+0 )
! 340: INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
! 341: $ LA_LINRX_CWISE_I
! 342: PARAMETER ( LA_LINRX_ITREF_I = 1,
! 343: $ LA_LINRX_ITHRESH_I = 2 )
! 344: PARAMETER ( LA_LINRX_CWISE_I = 3 )
! 345: INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
! 346: $ LA_LINRX_RCOND_I
! 347: PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
! 348: PARAMETER ( LA_LINRX_RCOND_I = 3 )
! 349: * ..
! 350: * .. Local Scalars ..
! 351: CHARACTER(1) NORM
! 352: LOGICAL ROWEQU, COLEQU, NOTRAN, IGNORE_CWISE
! 353: INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS,
! 354: $ ITHRESH
! 355: DOUBLE PRECISION ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND,
! 356: $ CWISE_WRONG, RTHRESH, UNSTABLE_THRESH
! 357: * ..
! 358: * .. External Subroutines ..
! 359: EXTERNAL XERBLA, ZGBCON, ZLA_GBRFSX_EXTENDED
! 360: * ..
! 361: * .. Intrinsic Functions ..
! 362: INTRINSIC MAX, SQRT, TRANSFER
! 363: * ..
! 364: * .. External Functions ..
! 365: EXTERNAL LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC
! 366: EXTERNAL DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
! 367: DOUBLE PRECISION DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
! 368: LOGICAL LSAME
! 369: INTEGER BLAS_FPINFO_X
! 370: INTEGER ILATRANS, ILAPREC
! 371: * ..
! 372: * .. Executable Statements ..
! 373: *
! 374: * Check the input parameters.
! 375: *
! 376: INFO = 0
! 377: TRANS_TYPE = ILATRANS( TRANS )
! 378: REF_TYPE = INT( ITREF_DEFAULT )
! 379: IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
! 380: IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
! 381: PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
! 382: ELSE
! 383: REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
! 384: END IF
! 385: END IF
! 386: *
! 387: * Set default parameters.
! 388: *
! 389: ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
! 390: ITHRESH = INT( ITHRESH_DEFAULT )
! 391: RTHRESH = RTHRESH_DEFAULT
! 392: UNSTABLE_THRESH = DZTHRESH_DEFAULT
! 393: IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
! 394: *
! 395: IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
! 396: IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
! 397: PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
! 398: ELSE
! 399: ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
! 400: END IF
! 401: END IF
! 402: IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
! 403: IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
! 404: IF ( IGNORE_CWISE ) THEN
! 405: PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
! 406: ELSE
! 407: PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
! 408: END IF
! 409: ELSE
! 410: IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
! 411: END IF
! 412: END IF
! 413: IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
! 414: N_NORMS = 0
! 415: ELSE IF ( IGNORE_CWISE ) THEN
! 416: N_NORMS = 1
! 417: ELSE
! 418: N_NORMS = 2
! 419: END IF
! 420: *
! 421: NOTRAN = LSAME( TRANS, 'N' )
! 422: ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
! 423: COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
! 424: *
! 425: * Test input parameters.
! 426: *
! 427: IF( TRANS_TYPE.EQ.-1 ) THEN
! 428: INFO = -1
! 429: ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
! 430: $ .NOT.LSAME( EQUED, 'N' ) ) THEN
! 431: INFO = -2
! 432: ELSE IF( N.LT.0 ) THEN
! 433: INFO = -3
! 434: ELSE IF( KL.LT.0 ) THEN
! 435: INFO = -4
! 436: ELSE IF( KU.LT.0 ) THEN
! 437: INFO = -5
! 438: ELSE IF( NRHS.LT.0 ) THEN
! 439: INFO = -6
! 440: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
! 441: INFO = -8
! 442: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
! 443: INFO = -10
! 444: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 445: INFO = -13
! 446: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
! 447: INFO = -15
! 448: END IF
! 449: IF( INFO.NE.0 ) THEN
! 450: CALL XERBLA( 'ZGBRFSX', -INFO )
! 451: RETURN
! 452: END IF
! 453: *
! 454: * Quick return if possible.
! 455: *
! 456: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
! 457: RCOND = 1.0D+0
! 458: DO J = 1, NRHS
! 459: BERR( J ) = 0.0D+0
! 460: IF ( N_ERR_BNDS .GE. 1 ) THEN
! 461: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 462: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 463: END IF
! 464: IF ( N_ERR_BNDS .GE. 2 ) THEN
! 465: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
! 466: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
! 467: END IF
! 468: IF ( N_ERR_BNDS .GE. 3 ) THEN
! 469: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
! 470: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
! 471: END IF
! 472: END DO
! 473: RETURN
! 474: END IF
! 475: *
! 476: * Default to failure.
! 477: *
! 478: RCOND = 0.0D+0
! 479: DO J = 1, NRHS
! 480: BERR( J ) = 1.0D+0
! 481: IF ( N_ERR_BNDS .GE. 1 ) THEN
! 482: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 483: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 484: END IF
! 485: IF ( N_ERR_BNDS .GE. 2 ) THEN
! 486: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
! 487: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
! 488: END IF
! 489: IF ( N_ERR_BNDS .GE. 3 ) THEN
! 490: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
! 491: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
! 492: END IF
! 493: END DO
! 494: *
! 495: * Compute the norm of A and the reciprocal of the condition
! 496: * number of A.
! 497: *
! 498: IF( NOTRAN ) THEN
! 499: NORM = 'I'
! 500: ELSE
! 501: NORM = '1'
! 502: END IF
! 503: ANORM = ZLANGB( NORM, N, KL, KU, AB, LDAB, RWORK )
! 504: CALL ZGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND,
! 505: $ WORK, RWORK, INFO )
! 506: *
! 507: * Perform refinement on each right-hand side
! 508: *
! 509: IF ( REF_TYPE .NE. 0 ) THEN
! 510:
! 511: PREC_TYPE = ILAPREC( 'E' )
! 512:
! 513: IF ( NOTRAN ) THEN
! 514: CALL ZLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
! 515: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B,
! 516: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
! 517: $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
! 518: $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
! 519: $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
! 520: $ INFO )
! 521: ELSE
! 522: CALL ZLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
! 523: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B,
! 524: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
! 525: $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
! 526: $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
! 527: $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
! 528: $ INFO )
! 529: END IF
! 530: END IF
! 531:
! 532: ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
! 533: IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN
! 534: *
! 535: * Compute scaled normwise condition number cond(A*C).
! 536: *
! 537: IF ( COLEQU .AND. NOTRAN ) THEN
! 538: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
! 539: $ LDAFB, IPIV, C, .TRUE., INFO, WORK, RWORK )
! 540: ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
! 541: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
! 542: $ LDAFB, IPIV, R, .TRUE., INFO, WORK, RWORK )
! 543: ELSE
! 544: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
! 545: $ LDAFB, IPIV, C, .FALSE., INFO, WORK, RWORK )
! 546: END IF
! 547: DO J = 1, NRHS
! 548: *
! 549: * Cap the error at 1.0.
! 550: *
! 551: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
! 552: $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0)
! 553: $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
! 554: *
! 555: * Threshold the error (see LAWN).
! 556: *
! 557: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
! 558: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
! 559: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
! 560: IF ( INFO .LE. N ) INFO = N + J
! 561: ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
! 562: $ THEN
! 563: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
! 564: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 565: END IF
! 566: *
! 567: * Save the condition number.
! 568: *
! 569: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
! 570: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
! 571: END IF
! 572:
! 573: END DO
! 574: END IF
! 575:
! 576: IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
! 577: *
! 578: * Compute componentwise condition number cond(A*diag(Y(:,J))) for
! 579: * each right-hand side using the current solution as an estimate of
! 580: * the true solution. If the componentwise error estimate is too
! 581: * large, then the solution is a lousy estimate of truth and the
! 582: * estimated RCOND may be too optimistic. To avoid misleading users,
! 583: * the inverse condition number is set to 0.0 when the estimated
! 584: * cwise error is at least CWISE_WRONG.
! 585: *
! 586: CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
! 587: DO J = 1, NRHS
! 588: IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
! 589: $ THEN
! 590: RCOND_TMP = ZLA_GBRCOND_X( TRANS, N, KL, KU, AB, LDAB,
! 591: $ AFB, LDAFB, IPIV, X( 1, J ), INFO, WORK, RWORK )
! 592: ELSE
! 593: RCOND_TMP = 0.0D+0
! 594: END IF
! 595: *
! 596: * Cap the error at 1.0.
! 597: *
! 598: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
! 599: $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
! 600: $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
! 601: *
! 602: * Threshold the error (see LAWN).
! 603: *
! 604: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
! 605: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
! 606: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
! 607: IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0
! 608: $ .AND. INFO.LT.N + J ) INFO = N + J
! 609: ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
! 610: $ .LT. ERR_LBND ) THEN
! 611: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
! 612: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
! 613: END IF
! 614: *
! 615: * Save the condition number.
! 616: *
! 617: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
! 618: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
! 619: END IF
! 620:
! 621: END DO
! 622: END IF
! 623: *
! 624: RETURN
! 625: *
! 626: * End of ZGBRFSX
! 627: *
! 628: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgbrfsx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack