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

1.16      bertrand    1: *> \brief <b> ZGTSVX computes the solution to system of linear equations A * X = B for GT matrices </b>
1.8       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.16      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.16      bertrand    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">
1.8       bertrand   15: *> [TXT]</a>
1.16      bertrand   16: *> \endhtmlonly
1.8       bertrand   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 )
1.16      bertrand   24: *
1.8       bertrand   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: *       ..
1.16      bertrand   37: *
1.8       bertrand   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
1.10      bertrand  137: *>          DLF is COMPLEX*16 array, dimension (N-1)
1.8       bertrand  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
1.10      bertrand  149: *>          DF is COMPLEX*16 array, dimension (N)
1.8       bertrand  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
1.10      bertrand  161: *>          DUF is COMPLEX*16 array, dimension (N-1)
1.8       bertrand  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
1.10      bertrand  171: *>          DU2 is COMPLEX*16 array, dimension (N-2)
1.8       bertrand  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
1.10      bertrand  183: *>          IPIV is INTEGER array, dimension (N)
1.8       bertrand  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: *
1.16      bertrand  283: *> \author Univ. of Tennessee
                    284: *> \author Univ. of California Berkeley
                    285: *> \author Univ. of Colorado Denver
                    286: *> \author NAG Ltd.
1.8       bertrand  287: *
1.12      bertrand  288: *> \ingroup complex16GTsolve
1.8       bertrand  289: *
                    290: *  =====================================================================
1.1       bertrand  291:       SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
                    292:      $                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
                    293:      $                   WORK, RWORK, INFO )
                    294: *
1.19    ! bertrand  295: *  -- LAPACK driver routine --
1.1       bertrand  296: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    297: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    298: *
                    299: *     .. Scalar Arguments ..
                    300:       CHARACTER          FACT, TRANS
                    301:       INTEGER            INFO, LDB, LDX, N, NRHS
                    302:       DOUBLE PRECISION   RCOND
                    303: *     ..
                    304: *     .. Array Arguments ..
                    305:       INTEGER            IPIV( * )
                    306:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
                    307:       COMPLEX*16         B( LDB, * ), D( * ), DF( * ), DL( * ),
                    308:      $                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
                    309:      $                   WORK( * ), X( LDX, * )
                    310: *     ..
                    311: *
                    312: *  =====================================================================
                    313: *
                    314: *     .. Parameters ..
                    315:       DOUBLE PRECISION   ZERO
                    316:       PARAMETER          ( ZERO = 0.0D+0 )
                    317: *     ..
                    318: *     .. Local Scalars ..
                    319:       LOGICAL            NOFACT, NOTRAN
                    320:       CHARACTER          NORM
                    321:       DOUBLE PRECISION   ANORM
                    322: *     ..
                    323: *     .. External Functions ..
                    324:       LOGICAL            LSAME
                    325:       DOUBLE PRECISION   DLAMCH, ZLANGT
                    326:       EXTERNAL           LSAME, DLAMCH, ZLANGT
                    327: *     ..
                    328: *     .. External Subroutines ..
                    329:       EXTERNAL           XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS,
                    330:      $                   ZLACPY
                    331: *     ..
                    332: *     .. Intrinsic Functions ..
                    333:       INTRINSIC          MAX
                    334: *     ..
                    335: *     .. Executable Statements ..
                    336: *
                    337:       INFO = 0
                    338:       NOFACT = LSAME( FACT, 'N' )
                    339:       NOTRAN = LSAME( TRANS, 'N' )
                    340:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
                    341:          INFO = -1
                    342:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
                    343:      $         LSAME( TRANS, 'C' ) ) THEN
                    344:          INFO = -2
                    345:       ELSE IF( N.LT.0 ) THEN
                    346:          INFO = -3
                    347:       ELSE IF( NRHS.LT.0 ) THEN
                    348:          INFO = -4
                    349:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    350:          INFO = -14
                    351:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
                    352:          INFO = -16
                    353:       END IF
                    354:       IF( INFO.NE.0 ) THEN
                    355:          CALL XERBLA( 'ZGTSVX', -INFO )
                    356:          RETURN
                    357:       END IF
                    358: *
                    359:       IF( NOFACT ) THEN
                    360: *
                    361: *        Compute the LU factorization of A.
                    362: *
                    363:          CALL ZCOPY( N, D, 1, DF, 1 )
                    364:          IF( N.GT.1 ) THEN
                    365:             CALL ZCOPY( N-1, DL, 1, DLF, 1 )
                    366:             CALL ZCOPY( N-1, DU, 1, DUF, 1 )
                    367:          END IF
                    368:          CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
                    369: *
                    370: *        Return if INFO is non-zero.
                    371: *
                    372:          IF( INFO.GT.0 )THEN
                    373:             RCOND = ZERO
                    374:             RETURN
                    375:          END IF
                    376:       END IF
                    377: *
                    378: *     Compute the norm of the matrix A.
                    379: *
                    380:       IF( NOTRAN ) THEN
                    381:          NORM = '1'
                    382:       ELSE
                    383:          NORM = 'I'
                    384:       END IF
                    385:       ANORM = ZLANGT( NORM, N, DL, D, DU )
                    386: *
                    387: *     Compute the reciprocal of the condition number of A.
                    388: *
                    389:       CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
                    390:      $             INFO )
                    391: *
                    392: *     Compute the solution vectors X.
                    393: *
                    394:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
                    395:       CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
                    396:      $             INFO )
                    397: *
                    398: *     Use iterative refinement to improve the computed solutions and
                    399: *     compute error bounds and backward error estimates for them.
                    400: *
                    401:       CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
                    402:      $             B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
                    403: *
                    404: *     Set INFO = N+1 if the matrix is singular to working precision.
                    405: *
                    406:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
                    407:      $   INFO = N + 1
                    408: *
                    409:       RETURN
                    410: *
                    411: *     End of ZGTSVX
                    412: *
                    413:       END

CVSweb interface <joel.bertrand@systella.fr>