Return to dgerfs.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DGERFS
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DGERFS + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgerfs.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgerfs.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgerfs.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
   22: *                          X, LDX, FERR, BERR, WORK, IWORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          TRANS
   26: *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IPIV( * ), IWORK( * )
   30: *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
   31: *      $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
   32: *       ..
   33: *  
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> DGERFS improves the computed solution to a system of linear
   41: *> equations and provides error bounds and backward error estimates for
   42: *> the solution.
   43: *> \endverbatim
   44: *
   45: *  Arguments:
   46: *  ==========
   47: *
   48: *> \param[in] TRANS
   49: *> \verbatim
   50: *>          TRANS is CHARACTER*1
   51: *>          Specifies the form of the system of equations:
   52: *>          = 'N':  A * X = B     (No transpose)
   53: *>          = 'T':  A**T * X = B  (Transpose)
   54: *>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
   55: *> \endverbatim
   56: *>
   57: *> \param[in] N
   58: *> \verbatim
   59: *>          N is INTEGER
   60: *>          The order of the matrix A.  N >= 0.
   61: *> \endverbatim
   62: *>
   63: *> \param[in] NRHS
   64: *> \verbatim
   65: *>          NRHS is INTEGER
   66: *>          The number of right hand sides, i.e., the number of columns
   67: *>          of the matrices B and X.  NRHS >= 0.
   68: *> \endverbatim
   69: *>
   70: *> \param[in] A
   71: *> \verbatim
   72: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   73: *>          The original N-by-N matrix A.
   74: *> \endverbatim
   75: *>
   76: *> \param[in] LDA
   77: *> \verbatim
   78: *>          LDA is INTEGER
   79: *>          The leading dimension of the array A.  LDA >= max(1,N).
   80: *> \endverbatim
   81: *>
   82: *> \param[in] AF
   83: *> \verbatim
   84: *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
   85: *>          The factors L and U from the factorization A = P*L*U
   86: *>          as computed by DGETRF.
   87: *> \endverbatim
   88: *>
   89: *> \param[in] LDAF
   90: *> \verbatim
   91: *>          LDAF is INTEGER
   92: *>          The leading dimension of the array AF.  LDAF >= max(1,N).
   93: *> \endverbatim
   94: *>
   95: *> \param[in] IPIV
   96: *> \verbatim
   97: *>          IPIV is INTEGER array, dimension (N)
   98: *>          The pivot indices from DGETRF; for 1<=i<=N, row i of the
   99: *>          matrix was interchanged with row IPIV(i).
  100: *> \endverbatim
  101: *>
  102: *> \param[in] B
  103: *> \verbatim
  104: *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  105: *>          The right hand side matrix B.
  106: *> \endverbatim
  107: *>
  108: *> \param[in] LDB
  109: *> \verbatim
  110: *>          LDB is INTEGER
  111: *>          The leading dimension of the array B.  LDB >= max(1,N).
  112: *> \endverbatim
  113: *>
  114: *> \param[in,out] X
  115: *> \verbatim
  116: *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
  117: *>          On entry, the solution matrix X, as computed by DGETRS.
  118: *>          On exit, the improved solution matrix X.
  119: *> \endverbatim
  120: *>
  121: *> \param[in] LDX
  122: *> \verbatim
  123: *>          LDX is INTEGER
  124: *>          The leading dimension of the array X.  LDX >= max(1,N).
  125: *> \endverbatim
  126: *>
  127: *> \param[out] FERR
  128: *> \verbatim
  129: *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
  130: *>          The estimated forward error bound for each solution vector
  131: *>          X(j) (the j-th column of the solution matrix X).
  132: *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
  133: *>          is an estimated upper bound for the magnitude of the largest
  134: *>          element in (X(j) - XTRUE) divided by the magnitude of the
  135: *>          largest element in X(j).  The estimate is as reliable as
  136: *>          the estimate for RCOND, and is almost always a slight
  137: *>          overestimate of the true error.
  138: *> \endverbatim
  139: *>
  140: *> \param[out] BERR
  141: *> \verbatim
  142: *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
  143: *>          The componentwise relative backward error of each solution
  144: *>          vector X(j) (i.e., the smallest relative change in
  145: *>          any element of A or B that makes X(j) an exact solution).
  146: *> \endverbatim
  147: *>
  148: *> \param[out] WORK
  149: *> \verbatim
  150: *>          WORK is DOUBLE PRECISION array, dimension (3*N)
  151: *> \endverbatim
  152: *>
  153: *> \param[out] IWORK
  154: *> \verbatim
  155: *>          IWORK is INTEGER array, dimension (N)
  156: *> \endverbatim
  157: *>
  158: *> \param[out] INFO
  159: *> \verbatim
  160: *>          INFO is INTEGER
  161: *>          = 0:  successful exit
  162: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  163: *> \endverbatim
  164: *
  165: *> \par Internal Parameters:
  166: *  =========================
  167: *>
  168: *> \verbatim
  169: *>  ITMAX is the maximum number of steps of iterative refinement.
  170: *> \endverbatim
  171: *
  172: *  Authors:
  173: *  ========
  174: *
  175: *> \author Univ. of Tennessee 
  176: *> \author Univ. of California Berkeley 
  177: *> \author Univ. of Colorado Denver 
  178: *> \author NAG Ltd. 
  179: *
  180: *> \date November 2011
  181: *
  182: *> \ingroup doubleGEcomputational
  183: *
  184: *  =====================================================================
  185:       SUBROUTINE DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
  186:      $                   X, LDX, FERR, BERR, WORK, IWORK, INFO )
  187: *
  188: *  -- LAPACK computational routine (version 3.4.0) --
  189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  191: *     November 2011
  192: *
  193: *     .. Scalar Arguments ..
  194:       CHARACTER          TRANS
  195:       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
  196: *     ..
  197: *     .. Array Arguments ..
  198:       INTEGER            IPIV( * ), IWORK( * )
  199:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  200:      $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
  201: *     ..
  202: *
  203: *  =====================================================================
  204: *
  205: *     .. Parameters ..
  206:       INTEGER            ITMAX
  207:       PARAMETER          ( ITMAX = 5 )
  208:       DOUBLE PRECISION   ZERO
  209:       PARAMETER          ( ZERO = 0.0D+0 )
  210:       DOUBLE PRECISION   ONE
  211:       PARAMETER          ( ONE = 1.0D+0 )
  212:       DOUBLE PRECISION   TWO
  213:       PARAMETER          ( TWO = 2.0D+0 )
  214:       DOUBLE PRECISION   THREE
  215:       PARAMETER          ( THREE = 3.0D+0 )
  216: *     ..
  217: *     .. Local Scalars ..
  218:       LOGICAL            NOTRAN
  219:       CHARACTER          TRANST
  220:       INTEGER            COUNT, I, J, K, KASE, NZ
  221:       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
  222: *     ..
  223: *     .. Local Arrays ..
  224:       INTEGER            ISAVE( 3 )
  225: *     ..
  226: *     .. External Subroutines ..
  227:       EXTERNAL           DAXPY, DCOPY, DGEMV, DGETRS, DLACN2, XERBLA
  228: *     ..
  229: *     .. Intrinsic Functions ..
  230:       INTRINSIC          ABS, MAX
  231: *     ..
  232: *     .. External Functions ..
  233:       LOGICAL            LSAME
  234:       DOUBLE PRECISION   DLAMCH
  235:       EXTERNAL           LSAME, DLAMCH
  236: *     ..
  237: *     .. Executable Statements ..
  238: *
  239: *     Test the input parameters.
  240: *
  241:       INFO = 0
  242:       NOTRAN = LSAME( TRANS, 'N' )
  243:       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
  244:      $    LSAME( TRANS, 'C' ) ) THEN
  245:          INFO = -1
  246:       ELSE IF( N.LT.0 ) THEN
  247:          INFO = -2
  248:       ELSE IF( NRHS.LT.0 ) THEN
  249:          INFO = -3
  250:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  251:          INFO = -5
  252:       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
  253:          INFO = -7
  254:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  255:          INFO = -10
  256:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  257:          INFO = -12
  258:       END IF
  259:       IF( INFO.NE.0 ) THEN
  260:          CALL XERBLA( 'DGERFS', -INFO )
  261:          RETURN
  262:       END IF
  263: *
  264: *     Quick return if possible
  265: *
  266:       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
  267:          DO 10 J = 1, NRHS
  268:             FERR( J ) = ZERO
  269:             BERR( J ) = ZERO
  270:    10    CONTINUE
  271:          RETURN
  272:       END IF
  273: *
  274:       IF( NOTRAN ) THEN
  275:          TRANST = 'T'
  276:       ELSE
  277:          TRANST = 'N'
  278:       END IF
  279: *
  280: *     NZ = maximum number of nonzero elements in each row of A, plus 1
  281: *
  282:       NZ = N + 1
  283:       EPS = DLAMCH( 'Epsilon' )
  284:       SAFMIN = DLAMCH( 'Safe minimum' )
  285:       SAFE1 = NZ*SAFMIN
  286:       SAFE2 = SAFE1 / EPS
  287: *
  288: *     Do for each right hand side
  289: *
  290:       DO 140 J = 1, NRHS
  291: *
  292:          COUNT = 1
  293:          LSTRES = THREE
  294:    20    CONTINUE
  295: *
  296: *        Loop until stopping criterion is satisfied.
  297: *
  298: *        Compute residual R = B - op(A) * X,
  299: *        where op(A) = A, A**T, or A**H, depending on TRANS.
  300: *
  301:          CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
  302:          CALL DGEMV( TRANS, N, N, -ONE, A, LDA, X( 1, J ), 1, ONE,
  303:      $               WORK( N+1 ), 1 )
  304: *
  305: *        Compute componentwise relative backward error from formula
  306: *
  307: *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
  308: *
  309: *        where abs(Z) is the componentwise absolute value of the matrix
  310: *        or vector Z.  If the i-th component of the denominator is less
  311: *        than SAFE2, then SAFE1 is added to the i-th components of the
  312: *        numerator and denominator before dividing.
  313: *
  314:          DO 30 I = 1, N
  315:             WORK( I ) = ABS( B( I, J ) )
  316:    30    CONTINUE
  317: *
  318: *        Compute abs(op(A))*abs(X) + abs(B).
  319: *
  320:          IF( NOTRAN ) THEN
  321:             DO 50 K = 1, N
  322:                XK = ABS( X( K, J ) )
  323:                DO 40 I = 1, N
  324:                   WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
  325:    40          CONTINUE
  326:    50       CONTINUE
  327:          ELSE
  328:             DO 70 K = 1, N
  329:                S = ZERO
  330:                DO 60 I = 1, N
  331:                   S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
  332:    60          CONTINUE
  333:                WORK( K ) = WORK( K ) + S
  334:    70       CONTINUE
  335:          END IF
  336:          S = ZERO
  337:          DO 80 I = 1, N
  338:             IF( WORK( I ).GT.SAFE2 ) THEN
  339:                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
  340:             ELSE
  341:                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
  342:      $             ( WORK( I )+SAFE1 ) )
  343:             END IF
  344:    80    CONTINUE
  345:          BERR( J ) = S
  346: *
  347: *        Test stopping criterion. Continue iterating if
  348: *           1) The residual BERR(J) is larger than machine epsilon, and
  349: *           2) BERR(J) decreased by at least a factor of 2 during the
  350: *              last iteration, and
  351: *           3) At most ITMAX iterations tried.
  352: *
  353:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
  354:      $       COUNT.LE.ITMAX ) THEN
  355: *
  356: *           Update solution and try again.
  357: *
  358:             CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK( N+1 ), N,
  359:      $                   INFO )
  360:             CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
  361:             LSTRES = BERR( J )
  362:             COUNT = COUNT + 1
  363:             GO TO 20
  364:          END IF
  365: *
  366: *        Bound error from formula
  367: *
  368: *        norm(X - XTRUE) / norm(X) .le. FERR =
  369: *        norm( abs(inv(op(A)))*
  370: *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
  371: *
  372: *        where
  373: *          norm(Z) is the magnitude of the largest component of Z
  374: *          inv(op(A)) is the inverse of op(A)
  375: *          abs(Z) is the componentwise absolute value of the matrix or
  376: *             vector Z
  377: *          NZ is the maximum number of nonzeros in any row of A, plus 1
  378: *          EPS is machine epsilon
  379: *
  380: *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
  381: *        is incremented by SAFE1 if the i-th component of
  382: *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
  383: *
  384: *        Use DLACN2 to estimate the infinity-norm of the matrix
  385: *           inv(op(A)) * diag(W),
  386: *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
  387: *
  388:          DO 90 I = 1, N
  389:             IF( WORK( I ).GT.SAFE2 ) THEN
  390:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
  391:             ELSE
  392:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
  393:             END IF
  394:    90    CONTINUE
  395: *
  396:          KASE = 0
  397:   100    CONTINUE
  398:          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
  399:      $                KASE, ISAVE )
  400:          IF( KASE.NE.0 ) THEN
  401:             IF( KASE.EQ.1 ) THEN
  402: *
  403: *              Multiply by diag(W)*inv(op(A)**T).
  404: *
  405:                CALL DGETRS( TRANST, N, 1, AF, LDAF, IPIV, WORK( N+1 ),
  406:      $                      N, INFO )
  407:                DO 110 I = 1, N
  408:                   WORK( N+I ) = WORK( I )*WORK( N+I )
  409:   110          CONTINUE
  410:             ELSE
  411: *
  412: *              Multiply by inv(op(A))*diag(W).
  413: *
  414:                DO 120 I = 1, N
  415:                   WORK( N+I ) = WORK( I )*WORK( N+I )
  416:   120          CONTINUE
  417:                CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK( N+1 ), N,
  418:      $                      INFO )
  419:             END IF
  420:             GO TO 100
  421:          END IF
  422: *
  423: *        Normalize error.
  424: *
  425:          LSTRES = ZERO
  426:          DO 130 I = 1, N
  427:             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
  428:   130    CONTINUE
  429:          IF( LSTRES.NE.ZERO )
  430:      $      FERR( J ) = FERR( J ) / LSTRES
  431: *
  432:   140 CONTINUE
  433: *
  434:       RETURN
  435: *
  436: *     End of DGERFS
  437: *
  438:       END