Annotation of rpl/lapack/lapack/zgtsvx.f, revision 1.8

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

CVSweb interface <joel.bertrand@systella.fr>