Return to dsyrfs.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

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