Return to dporfs.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

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