Return to dpbrfs.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b DPBRFS
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DPBRFS + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbrfs.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbrfs.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbrfs.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B,
   22: *                          LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO
   26: *       INTEGER            INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IWORK( * )
   30: *       DOUBLE PRECISION   AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
   31: *      $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
   32: *       ..
   33: *  
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> DPBRFS improves the computed solution to a system of linear
   41: *> equations when the coefficient matrix is symmetric positive definite
   42: *> and banded, 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] KD
   63: *> \verbatim
   64: *>          KD is INTEGER
   65: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
   66: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
   67: *> \endverbatim
   68: *>
   69: *> \param[in] NRHS
   70: *> \verbatim
   71: *>          NRHS is INTEGER
   72: *>          The number of right hand sides, i.e., the number of columns
   73: *>          of the matrices B and X.  NRHS >= 0.
   74: *> \endverbatim
   75: *>
   76: *> \param[in] AB
   77: *> \verbatim
   78: *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
   79: *>          The upper or lower triangle of the symmetric band matrix A,
   80: *>          stored in the first KD+1 rows of the array.  The j-th column
   81: *>          of A is stored in the j-th column of the array AB as follows:
   82: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
   83: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
   84: *> \endverbatim
   85: *>
   86: *> \param[in] LDAB
   87: *> \verbatim
   88: *>          LDAB is INTEGER
   89: *>          The leading dimension of the array AB.  LDAB >= KD+1.
   90: *> \endverbatim
   91: *>
   92: *> \param[in] AFB
   93: *> \verbatim
   94: *>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
   95: *>          The triangular factor U or L from the Cholesky factorization
   96: *>          A = U**T*U or A = L*L**T of the band matrix A as computed by
   97: *>          DPBTRF, in the same storage format as A (see AB).
   98: *> \endverbatim
   99: *>
  100: *> \param[in] LDAFB
  101: *> \verbatim
  102: *>          LDAFB is INTEGER
  103: *>          The leading dimension of the array AFB.  LDAFB >= KD+1.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] B
  107: *> \verbatim
  108: *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  109: *>          The right hand side matrix B.
  110: *> \endverbatim
  111: *>
  112: *> \param[in] LDB
  113: *> \verbatim
  114: *>          LDB is INTEGER
  115: *>          The leading dimension of the array B.  LDB >= max(1,N).
  116: *> \endverbatim
  117: *>
  118: *> \param[in,out] X
  119: *> \verbatim
  120: *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
  121: *>          On entry, the solution matrix X, as computed by DPBTRS.
  122: *>          On exit, the improved solution matrix X.
  123: *> \endverbatim
  124: *>
  125: *> \param[in] LDX
  126: *> \verbatim
  127: *>          LDX is INTEGER
  128: *>          The leading dimension of the array X.  LDX >= max(1,N).
  129: *> \endverbatim
  130: *>
  131: *> \param[out] FERR
  132: *> \verbatim
  133: *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
  134: *>          The estimated forward error bound for each solution vector
  135: *>          X(j) (the j-th column of the solution matrix X).
  136: *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
  137: *>          is an estimated upper bound for the magnitude of the largest
  138: *>          element in (X(j) - XTRUE) divided by the magnitude of the
  139: *>          largest element in X(j).  The estimate is as reliable as
  140: *>          the estimate for RCOND, and is almost always a slight
  141: *>          overestimate of the true error.
  142: *> \endverbatim
  143: *>
  144: *> \param[out] BERR
  145: *> \verbatim
  146: *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
  147: *>          The componentwise relative backward error of each solution
  148: *>          vector X(j) (i.e., the smallest relative change in
  149: *>          any element of A or B that makes X(j) an exact solution).
  150: *> \endverbatim
  151: *>
  152: *> \param[out] WORK
  153: *> \verbatim
  154: *>          WORK is DOUBLE PRECISION array, dimension (3*N)
  155: *> \endverbatim
  156: *>
  157: *> \param[out] IWORK
  158: *> \verbatim
  159: *>          IWORK is INTEGER array, dimension (N)
  160: *> \endverbatim
  161: *>
  162: *> \param[out] INFO
  163: *> \verbatim
  164: *>          INFO is INTEGER
  165: *>          = 0:  successful exit
  166: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  167: *> \endverbatim
  168: *
  169: *> \par Internal Parameters:
  170: *  =========================
  171: *>
  172: *> \verbatim
  173: *>  ITMAX is the maximum number of steps of iterative refinement.
  174: *> \endverbatim
  175: *
  176: *  Authors:
  177: *  ========
  178: *
  179: *> \author Univ. of Tennessee 
  180: *> \author Univ. of California Berkeley 
  181: *> \author Univ. of Colorado Denver 
  182: *> \author NAG Ltd. 
  183: *
  184: *> \date November 2011
  185: *
  186: *> \ingroup doubleOTHERcomputational
  187: *
  188: *  =====================================================================
  189:       SUBROUTINE DPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B,
  190:      $                   LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
  191: *
  192: *  -- LAPACK computational routine (version 3.4.0) --
  193: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  194: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  195: *     November 2011
  196: *
  197: *     .. Scalar Arguments ..
  198:       CHARACTER          UPLO
  199:       INTEGER            INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
  200: *     ..
  201: *     .. Array Arguments ..
  202:       INTEGER            IWORK( * )
  203:       DOUBLE PRECISION   AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
  204:      $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
  205: *     ..
  206: *
  207: *  =====================================================================
  208: *
  209: *     .. Parameters ..
  210:       INTEGER            ITMAX
  211:       PARAMETER          ( ITMAX = 5 )
  212:       DOUBLE PRECISION   ZERO
  213:       PARAMETER          ( ZERO = 0.0D+0 )
  214:       DOUBLE PRECISION   ONE
  215:       PARAMETER          ( ONE = 1.0D+0 )
  216:       DOUBLE PRECISION   TWO
  217:       PARAMETER          ( TWO = 2.0D+0 )
  218:       DOUBLE PRECISION   THREE
  219:       PARAMETER          ( THREE = 3.0D+0 )
  220: *     ..
  221: *     .. Local Scalars ..
  222:       LOGICAL            UPPER
  223:       INTEGER            COUNT, I, J, K, KASE, L, NZ
  224:       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
  225: *     ..
  226: *     .. Local Arrays ..
  227:       INTEGER            ISAVE( 3 )
  228: *     ..
  229: *     .. External Subroutines ..
  230:       EXTERNAL           DAXPY, DCOPY, DLACN2, DPBTRS, DSBMV, XERBLA
  231: *     ..
  232: *     .. Intrinsic Functions ..
  233:       INTRINSIC          ABS, MAX, MIN
  234: *     ..
  235: *     .. External Functions ..
  236:       LOGICAL            LSAME
  237:       DOUBLE PRECISION   DLAMCH
  238:       EXTERNAL           LSAME, DLAMCH
  239: *     ..
  240: *     .. Executable Statements ..
  241: *
  242: *     Test the input parameters.
  243: *
  244:       INFO = 0
  245:       UPPER = LSAME( UPLO, 'U' )
  246:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  247:          INFO = -1
  248:       ELSE IF( N.LT.0 ) THEN
  249:          INFO = -2
  250:       ELSE IF( KD.LT.0 ) THEN
  251:          INFO = -3
  252:       ELSE IF( NRHS.LT.0 ) THEN
  253:          INFO = -4
  254:       ELSE IF( LDAB.LT.KD+1 ) THEN
  255:          INFO = -6
  256:       ELSE IF( LDAFB.LT.KD+1 ) THEN
  257:          INFO = -8
  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( 'DPBRFS', -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 = MIN( N+1, 2*KD+2 )
  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 DSBMV( UPLO, N, KD, -ONE, AB, LDAB, 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:                L = KD + 1 - K
  322:                DO 40 I = MAX( 1, K-KD ), K - 1
  323:                   WORK( I ) = WORK( I ) + ABS( AB( L+I, K ) )*XK
  324:                   S = S + ABS( AB( L+I, K ) )*ABS( X( I, J ) )
  325:    40          CONTINUE
  326:                WORK( K ) = WORK( K ) + ABS( AB( KD+1, K ) )*XK + S
  327:    50       CONTINUE
  328:          ELSE
  329:             DO 70 K = 1, N
  330:                S = ZERO
  331:                XK = ABS( X( K, J ) )
  332:                WORK( K ) = WORK( K ) + ABS( AB( 1, K ) )*XK
  333:                L = 1 - K
  334:                DO 60 I = K + 1, MIN( N, K+KD )
  335:                   WORK( I ) = WORK( I ) + ABS( AB( L+I, K ) )*XK
  336:                   S = S + ABS( AB( L+I, K ) )*ABS( X( I, J ) )
  337:    60          CONTINUE
  338:                WORK( K ) = WORK( K ) + S
  339:    70       CONTINUE
  340:          END IF
  341:          S = ZERO
  342:          DO 80 I = 1, N
  343:             IF( WORK( I ).GT.SAFE2 ) THEN
  344:                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
  345:             ELSE
  346:                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
  347:      $             ( WORK( I )+SAFE1 ) )
  348:             END IF
  349:    80    CONTINUE
  350:          BERR( J ) = S
  351: *
  352: *        Test stopping criterion. Continue iterating if
  353: *           1) The residual BERR(J) is larger than machine epsilon, and
  354: *           2) BERR(J) decreased by at least a factor of 2 during the
  355: *              last iteration, and
  356: *           3) At most ITMAX iterations tried.
  357: *
  358:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
  359:      $       COUNT.LE.ITMAX ) THEN
  360: *
  361: *           Update solution and try again.
  362: *
  363:             CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
  364:      $                   INFO )
  365:             CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
  366:             LSTRES = BERR( J )
  367:             COUNT = COUNT + 1
  368:             GO TO 20
  369:          END IF
  370: *
  371: *        Bound error from formula
  372: *
  373: *        norm(X - XTRUE) / norm(X) .le. FERR =
  374: *        norm( abs(inv(A))*
  375: *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
  376: *
  377: *        where
  378: *          norm(Z) is the magnitude of the largest component of Z
  379: *          inv(A) is the inverse of A
  380: *          abs(Z) is the componentwise absolute value of the matrix or
  381: *             vector Z
  382: *          NZ is the maximum number of nonzeros in any row of A, plus 1
  383: *          EPS is machine epsilon
  384: *
  385: *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
  386: *        is incremented by SAFE1 if the i-th component of
  387: *        abs(A)*abs(X) + abs(B) is less than SAFE2.
  388: *
  389: *        Use DLACN2 to estimate the infinity-norm of the matrix
  390: *           inv(A) * diag(W),
  391: *        where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
  392: *
  393:          DO 90 I = 1, N
  394:             IF( WORK( I ).GT.SAFE2 ) THEN
  395:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
  396:             ELSE
  397:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
  398:             END IF
  399:    90    CONTINUE
  400: *
  401:          KASE = 0
  402:   100    CONTINUE
  403:          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
  404:      $                KASE, ISAVE )
  405:          IF( KASE.NE.0 ) THEN
  406:             IF( KASE.EQ.1 ) THEN
  407: *
  408: *              Multiply by diag(W)*inv(A**T).
  409: *
  410:                CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
  411:      $                      INFO )
  412:                DO 110 I = 1, N
  413:                   WORK( N+I ) = WORK( N+I )*WORK( I )
  414:   110          CONTINUE
  415:             ELSE IF( KASE.EQ.2 ) THEN
  416: *
  417: *              Multiply by inv(A)*diag(W).
  418: *
  419:                DO 120 I = 1, N
  420:                   WORK( N+I ) = WORK( N+I )*WORK( I )
  421:   120          CONTINUE
  422:                CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
  423:      $                      INFO )
  424:             END IF
  425:             GO TO 100
  426:          END IF
  427: *
  428: *        Normalize error.
  429: *
  430:          LSTRES = ZERO
  431:          DO 130 I = 1, N
  432:             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
  433:   130    CONTINUE
  434:          IF( LSTRES.NE.ZERO )
  435:      $      FERR( J ) = FERR( J ) / LSTRES
  436: *
  437:   140 CONTINUE
  438: *
  439:       RETURN
  440: *
  441: *     End of DPBRFS
  442: *
  443:       END