Annotation of rpl/lapack/lapack/zherfs.f, revision 1.9

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

CVSweb interface <joel.bertrand@systella.fr>