Return to zhprfs.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

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