Return to dsprfs.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

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