Return to ztprfs.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

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