Annotation of rpl/lapack/lapack/dpprfs.f, revision 1.12

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

CVSweb interface <joel.bertrand@systella.fr>