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

1.9     ! bertrand    1: *> \brief \b DPTSVX
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DPTSVX + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptsvx.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptsvx.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptsvx.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
        !            22: *                          RCOND, FERR, BERR, WORK, INFO )
        !            23: * 
        !            24: *       .. Scalar Arguments ..
        !            25: *       CHARACTER          FACT
        !            26: *       INTEGER            INFO, LDB, LDX, N, NRHS
        !            27: *       DOUBLE PRECISION   RCOND
        !            28: *       ..
        !            29: *       .. Array Arguments ..
        !            30: *       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
        !            31: *      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
        !            32: *      $                   X( LDX, * )
        !            33: *       ..
        !            34: *  
        !            35: *
        !            36: *> \par Purpose:
        !            37: *  =============
        !            38: *>
        !            39: *> \verbatim
        !            40: *>
        !            41: *> DPTSVX uses the factorization A = L*D*L**T to compute the solution
        !            42: *> to a real system of linear equations A*X = B, where A is an N-by-N
        !            43: *> symmetric positive definite tridiagonal matrix and X and B are
        !            44: *> N-by-NRHS matrices.
        !            45: *>
        !            46: *> Error bounds on the solution and a condition estimate are also
        !            47: *> provided.
        !            48: *> \endverbatim
        !            49: *
        !            50: *> \par Description:
        !            51: *  =================
        !            52: *>
        !            53: *> \verbatim
        !            54: *>
        !            55: *> The following steps are performed:
        !            56: *>
        !            57: *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L
        !            58: *>    is a unit lower bidiagonal matrix and D is diagonal.  The
        !            59: *>    factorization can also be regarded as having the form
        !            60: *>    A = U**T*D*U.
        !            61: *>
        !            62: *> 2. If the leading i-by-i principal minor is not positive definite,
        !            63: *>    then the routine returns with INFO = i. Otherwise, the factored
        !            64: *>    form of A is used to estimate the condition number of the matrix
        !            65: *>    A.  If the reciprocal of the condition number is less than machine
        !            66: *>    precision, INFO = N+1 is returned as a warning, but the routine
        !            67: *>    still goes on to solve for X and compute error bounds as
        !            68: *>    described below.
        !            69: *>
        !            70: *> 3. The system of equations is solved for X using the factored form
        !            71: *>    of A.
        !            72: *>
        !            73: *> 4. Iterative refinement is applied to improve the computed solution
        !            74: *>    matrix and calculate error bounds and backward error estimates
        !            75: *>    for it.
        !            76: *> \endverbatim
        !            77: *
        !            78: *  Arguments:
        !            79: *  ==========
        !            80: *
        !            81: *> \param[in] FACT
        !            82: *> \verbatim
        !            83: *>          FACT is CHARACTER*1
        !            84: *>          Specifies whether or not the factored form of A has been
        !            85: *>          supplied on entry.
        !            86: *>          = 'F':  On entry, DF and EF contain the factored form of A.
        !            87: *>                  D, E, DF, and EF will not be modified.
        !            88: *>          = 'N':  The matrix A will be copied to DF and EF and
        !            89: *>                  factored.
        !            90: *> \endverbatim
        !            91: *>
        !            92: *> \param[in] N
        !            93: *> \verbatim
        !            94: *>          N is INTEGER
        !            95: *>          The order of the matrix A.  N >= 0.
        !            96: *> \endverbatim
        !            97: *>
        !            98: *> \param[in] NRHS
        !            99: *> \verbatim
        !           100: *>          NRHS is INTEGER
        !           101: *>          The number of right hand sides, i.e., the number of columns
        !           102: *>          of the matrices B and X.  NRHS >= 0.
        !           103: *> \endverbatim
        !           104: *>
        !           105: *> \param[in] D
        !           106: *> \verbatim
        !           107: *>          D is DOUBLE PRECISION array, dimension (N)
        !           108: *>          The n diagonal elements of the tridiagonal matrix A.
        !           109: *> \endverbatim
        !           110: *>
        !           111: *> \param[in] E
        !           112: *> \verbatim
        !           113: *>          E is DOUBLE PRECISION array, dimension (N-1)
        !           114: *>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
        !           115: *> \endverbatim
        !           116: *>
        !           117: *> \param[in,out] DF
        !           118: *> \verbatim
        !           119: *>          DF is or output) DOUBLE PRECISION array, dimension (N)
        !           120: *>          If FACT = 'F', then DF is an input argument and on entry
        !           121: *>          contains the n diagonal elements of the diagonal matrix D
        !           122: *>          from the L*D*L**T factorization of A.
        !           123: *>          If FACT = 'N', then DF is an output argument and on exit
        !           124: *>          contains the n diagonal elements of the diagonal matrix D
        !           125: *>          from the L*D*L**T factorization of A.
        !           126: *> \endverbatim
        !           127: *>
        !           128: *> \param[in,out] EF
        !           129: *> \verbatim
        !           130: *>          EF is or output) DOUBLE PRECISION array, dimension (N-1)
        !           131: *>          If FACT = 'F', then EF is an input argument and on entry
        !           132: *>          contains the (n-1) subdiagonal elements of the unit
        !           133: *>          bidiagonal factor L from the L*D*L**T factorization of A.
        !           134: *>          If FACT = 'N', then EF is an output argument and on exit
        !           135: *>          contains the (n-1) subdiagonal elements of the unit
        !           136: *>          bidiagonal factor L from the L*D*L**T factorization of A.
        !           137: *> \endverbatim
        !           138: *>
        !           139: *> \param[in] B
        !           140: *> \verbatim
        !           141: *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
        !           142: *>          The N-by-NRHS right hand side matrix B.
        !           143: *> \endverbatim
        !           144: *>
        !           145: *> \param[in] LDB
        !           146: *> \verbatim
        !           147: *>          LDB is INTEGER
        !           148: *>          The leading dimension of the array B.  LDB >= max(1,N).
        !           149: *> \endverbatim
        !           150: *>
        !           151: *> \param[out] X
        !           152: *> \verbatim
        !           153: *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
        !           154: *>          If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X.
        !           155: *> \endverbatim
        !           156: *>
        !           157: *> \param[in] LDX
        !           158: *> \verbatim
        !           159: *>          LDX is INTEGER
        !           160: *>          The leading dimension of the array X.  LDX >= max(1,N).
        !           161: *> \endverbatim
        !           162: *>
        !           163: *> \param[out] RCOND
        !           164: *> \verbatim
        !           165: *>          RCOND is DOUBLE PRECISION
        !           166: *>          The reciprocal condition number of the matrix A.  If RCOND
        !           167: *>          is less than the machine precision (in particular, if
        !           168: *>          RCOND = 0), the matrix is singular to working precision.
        !           169: *>          This condition is indicated by a return code of INFO > 0.
        !           170: *> \endverbatim
        !           171: *>
        !           172: *> \param[out] FERR
        !           173: *> \verbatim
        !           174: *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
        !           175: *>          The forward error bound for each solution vector
        !           176: *>          X(j) (the j-th column of the solution matrix X).
        !           177: *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
        !           178: *>          is an estimated upper bound for the magnitude of the largest
        !           179: *>          element in (X(j) - XTRUE) divided by the magnitude of the
        !           180: *>          largest element in X(j).
        !           181: *> \endverbatim
        !           182: *>
        !           183: *> \param[out] BERR
        !           184: *> \verbatim
        !           185: *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
        !           186: *>          The componentwise relative backward error of each solution
        !           187: *>          vector X(j) (i.e., the smallest relative change in any
        !           188: *>          element of A or B that makes X(j) an exact solution).
        !           189: *> \endverbatim
        !           190: *>
        !           191: *> \param[out] WORK
        !           192: *> \verbatim
        !           193: *>          WORK is DOUBLE PRECISION array, dimension (2*N)
        !           194: *> \endverbatim
        !           195: *>
        !           196: *> \param[out] INFO
        !           197: *> \verbatim
        !           198: *>          INFO is INTEGER
        !           199: *>          = 0:  successful exit
        !           200: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
        !           201: *>          > 0:  if INFO = i, and i is
        !           202: *>                <= N:  the leading minor of order i of A is
        !           203: *>                       not positive definite, so the factorization
        !           204: *>                       could not be completed, and the solution has not
        !           205: *>                       been computed. RCOND = 0 is returned.
        !           206: *>                = N+1: U is nonsingular, but RCOND is less than machine
        !           207: *>                       precision, meaning that the matrix is singular
        !           208: *>                       to working precision.  Nevertheless, the
        !           209: *>                       solution and error bounds are computed because
        !           210: *>                       there are a number of situations where the
        !           211: *>                       computed solution can be more accurate than the
        !           212: *>                       value of RCOND would suggest.
        !           213: *> \endverbatim
        !           214: *
        !           215: *  Authors:
        !           216: *  ========
        !           217: *
        !           218: *> \author Univ. of Tennessee 
        !           219: *> \author Univ. of California Berkeley 
        !           220: *> \author Univ. of Colorado Denver 
        !           221: *> \author NAG Ltd. 
        !           222: *
        !           223: *> \date November 2011
        !           224: *
        !           225: *> \ingroup doubleOTHERcomputational
        !           226: *
        !           227: *  =====================================================================
1.1       bertrand  228:       SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
                    229:      $                   RCOND, FERR, BERR, WORK, INFO )
                    230: *
1.9     ! bertrand  231: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  232: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    233: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9     ! bertrand  234: *     November 2011
1.1       bertrand  235: *
                    236: *     .. Scalar Arguments ..
                    237:       CHARACTER          FACT
                    238:       INTEGER            INFO, LDB, LDX, N, NRHS
                    239:       DOUBLE PRECISION   RCOND
                    240: *     ..
                    241: *     .. Array Arguments ..
                    242:       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
                    243:      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
                    244:      $                   X( LDX, * )
                    245: *     ..
                    246: *
                    247: *  =====================================================================
                    248: *
                    249: *     .. Parameters ..
                    250:       DOUBLE PRECISION   ZERO
                    251:       PARAMETER          ( ZERO = 0.0D+0 )
                    252: *     ..
                    253: *     .. Local Scalars ..
                    254:       LOGICAL            NOFACT
                    255:       DOUBLE PRECISION   ANORM
                    256: *     ..
                    257: *     .. External Functions ..
                    258:       LOGICAL            LSAME
                    259:       DOUBLE PRECISION   DLAMCH, DLANST
                    260:       EXTERNAL           LSAME, DLAMCH, DLANST
                    261: *     ..
                    262: *     .. External Subroutines ..
                    263:       EXTERNAL           DCOPY, DLACPY, DPTCON, DPTRFS, DPTTRF, DPTTRS,
                    264:      $                   XERBLA
                    265: *     ..
                    266: *     .. Intrinsic Functions ..
                    267:       INTRINSIC          MAX
                    268: *     ..
                    269: *     .. Executable Statements ..
                    270: *
                    271: *     Test the input parameters.
                    272: *
                    273:       INFO = 0
                    274:       NOFACT = LSAME( FACT, 'N' )
                    275:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
                    276:          INFO = -1
                    277:       ELSE IF( N.LT.0 ) THEN
                    278:          INFO = -2
                    279:       ELSE IF( NRHS.LT.0 ) THEN
                    280:          INFO = -3
                    281:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    282:          INFO = -9
                    283:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
                    284:          INFO = -11
                    285:       END IF
                    286:       IF( INFO.NE.0 ) THEN
                    287:          CALL XERBLA( 'DPTSVX', -INFO )
                    288:          RETURN
                    289:       END IF
                    290: *
                    291:       IF( NOFACT ) THEN
                    292: *
1.8       bertrand  293: *        Compute the L*D*L**T (or U**T*D*U) factorization of A.
1.1       bertrand  294: *
                    295:          CALL DCOPY( N, D, 1, DF, 1 )
                    296:          IF( N.GT.1 )
                    297:      $      CALL DCOPY( N-1, E, 1, EF, 1 )
                    298:          CALL DPTTRF( N, DF, EF, INFO )
                    299: *
                    300: *        Return if INFO is non-zero.
                    301: *
                    302:          IF( INFO.GT.0 )THEN
                    303:             RCOND = ZERO
                    304:             RETURN
                    305:          END IF
                    306:       END IF
                    307: *
                    308: *     Compute the norm of the matrix A.
                    309: *
                    310:       ANORM = DLANST( '1', N, D, E )
                    311: *
                    312: *     Compute the reciprocal of the condition number of A.
                    313: *
                    314:       CALL DPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO )
                    315: *
                    316: *     Compute the solution vectors X.
                    317: *
                    318:       CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
                    319:       CALL DPTTRS( N, NRHS, DF, EF, X, LDX, INFO )
                    320: *
                    321: *     Use iterative refinement to improve the computed solutions and
                    322: *     compute error bounds and backward error estimates for them.
                    323: *
                    324:       CALL DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR,
                    325:      $             WORK, INFO )
                    326: *
                    327: *     Set INFO = N+1 if the matrix is singular to working precision.
                    328: *
                    329:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
                    330:      $   INFO = N + 1
                    331: *
                    332:       RETURN
                    333: *
                    334: *     End of DPTSVX
                    335: *
                    336:       END

CVSweb interface <joel.bertrand@systella.fr>