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

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

CVSweb interface <joel.bertrand@systella.fr>