LAPACK 3.3.0

# dptsvx.f

Go to the documentation of this file.
```00001       SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
00002      \$                   RCOND, FERR, BERR, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          FACT
00011       INTEGER            INFO, LDB, LDX, N, NRHS
00012       DOUBLE PRECISION   RCOND
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
00016      \$                   E( * ), EF( * ), FERR( * ), WORK( * ),
00017      \$                   X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DPTSVX uses the factorization A = L*D*L**T to compute the solution
00024 *  to a real system of linear equations A*X = B, where A is an N-by-N
00025 *  symmetric positive definite tridiagonal matrix and X and B are
00026 *  N-by-NRHS matrices.
00027 *
00028 *  Error bounds on the solution and a condition estimate are also
00029 *  provided.
00030 *
00031 *  Description
00032 *  ===========
00033 *
00034 *  The following steps are performed:
00035 *
00036 *  1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L
00037 *     is a unit lower bidiagonal matrix and D is diagonal.  The
00038 *     factorization can also be regarded as having the form
00039 *     A = U**T*D*U.
00040 *
00041 *  2. If the leading i-by-i principal minor is not positive definite,
00042 *     then the routine returns with INFO = i. Otherwise, the factored
00043 *     form of A is used to estimate the condition number of the matrix
00044 *     A.  If the reciprocal of the condition number is less than machine
00045 *     precision, INFO = N+1 is returned as a warning, but the routine
00046 *     still goes on to solve for X and compute error bounds as
00047 *     described below.
00048 *
00049 *  3. The system of equations is solved for X using the factored form
00050 *     of A.
00051 *
00052 *  4. Iterative refinement is applied to improve the computed solution
00053 *     matrix and calculate error bounds and backward error estimates
00054 *     for it.
00055 *
00056 *  Arguments
00057 *  =========
00058 *
00059 *  FACT    (input) CHARACTER*1
00060 *          Specifies whether or not the factored form of A has been
00061 *          supplied on entry.
00062 *          = 'F':  On entry, DF and EF contain the factored form of A.
00063 *                  D, E, DF, and EF will not be modified.
00064 *          = 'N':  The matrix A will be copied to DF and EF and
00065 *                  factored.
00066 *
00067 *  N       (input) INTEGER
00068 *          The order of the matrix A.  N >= 0.
00069 *
00070 *  NRHS    (input) INTEGER
00071 *          The number of right hand sides, i.e., the number of columns
00072 *          of the matrices B and X.  NRHS >= 0.
00073 *
00074 *  D       (input) DOUBLE PRECISION array, dimension (N)
00075 *          The n diagonal elements of the tridiagonal matrix A.
00076 *
00077 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00078 *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
00079 *
00080 *  DF      (input or output) DOUBLE PRECISION array, dimension (N)
00081 *          If FACT = 'F', then DF is an input argument and on entry
00082 *          contains the n diagonal elements of the diagonal matrix D
00083 *          from the L*D*L**T factorization of A.
00084 *          If FACT = 'N', then DF is an output argument and on exit
00085 *          contains the n diagonal elements of the diagonal matrix D
00086 *          from the L*D*L**T factorization of A.
00087 *
00088 *  EF      (input or output) DOUBLE PRECISION array, dimension (N-1)
00089 *          If FACT = 'F', then EF is an input argument and on entry
00090 *          contains the (n-1) subdiagonal elements of the unit
00091 *          bidiagonal factor L from the L*D*L**T factorization of A.
00092 *          If FACT = 'N', then EF is an output argument and on exit
00093 *          contains the (n-1) subdiagonal elements of the unit
00094 *          bidiagonal factor L from the L*D*L**T factorization of A.
00095 *
00096 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
00097 *          The N-by-NRHS right hand side matrix B.
00098 *
00099 *  LDB     (input) INTEGER
00100 *          The leading dimension of the array B.  LDB >= max(1,N).
00101 *
00102 *  X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
00103 *          If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X.
00104 *
00105 *  LDX     (input) INTEGER
00106 *          The leading dimension of the array X.  LDX >= max(1,N).
00107 *
00108 *  RCOND   (output) DOUBLE PRECISION
00109 *          The reciprocal condition number of the matrix A.  If RCOND
00110 *          is less than the machine precision (in particular, if
00111 *          RCOND = 0), the matrix is singular to working precision.
00112 *          This condition is indicated by a return code of INFO > 0.
00113 *
00114 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00115 *          The forward error bound for each solution vector
00116 *          X(j) (the j-th column of the solution matrix X).
00117 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00118 *          is an estimated upper bound for the magnitude of the largest
00119 *          element in (X(j) - XTRUE) divided by the magnitude of the
00120 *          largest element in X(j).
00121 *
00122 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00123 *          The componentwise relative backward error of each solution
00124 *          vector X(j) (i.e., the smallest relative change in any
00125 *          element of A or B that makes X(j) an exact solution).
00126 *
00127 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00128 *
00129 *  INFO    (output) INTEGER
00130 *          = 0:  successful exit
00131 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00132 *          > 0:  if INFO = i, and i is
00133 *                <= N:  the leading minor of order i of A is
00134 *                       not positive definite, so the factorization
00135 *                       could not be completed, and the solution has not
00136 *                       been computed. RCOND = 0 is returned.
00137 *                = N+1: U is nonsingular, but RCOND is less than machine
00138 *                       precision, meaning that the matrix is singular
00139 *                       to working precision.  Nevertheless, the
00140 *                       solution and error bounds are computed because
00141 *                       there are a number of situations where the
00142 *                       computed solution can be more accurate than the
00143 *                       value of RCOND would suggest.
00144 *
00145 *  =====================================================================
00146 *
00147 *     .. Parameters ..
00148       DOUBLE PRECISION   ZERO
00149       PARAMETER          ( ZERO = 0.0D+0 )
00150 *     ..
00151 *     .. Local Scalars ..
00152       LOGICAL            NOFACT
00153       DOUBLE PRECISION   ANORM
00154 *     ..
00155 *     .. External Functions ..
00156       LOGICAL            LSAME
00157       DOUBLE PRECISION   DLAMCH, DLANST
00158       EXTERNAL           LSAME, DLAMCH, DLANST
00159 *     ..
00160 *     .. External Subroutines ..
00161       EXTERNAL           DCOPY, DLACPY, DPTCON, DPTRFS, DPTTRF, DPTTRS,
00162      \$                   XERBLA
00163 *     ..
00164 *     .. Intrinsic Functions ..
00165       INTRINSIC          MAX
00166 *     ..
00167 *     .. Executable Statements ..
00168 *
00169 *     Test the input parameters.
00170 *
00171       INFO = 0
00172       NOFACT = LSAME( FACT, 'N' )
00173       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
00174          INFO = -1
00175       ELSE IF( N.LT.0 ) THEN
00176          INFO = -2
00177       ELSE IF( NRHS.LT.0 ) THEN
00178          INFO = -3
00179       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00180          INFO = -9
00181       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00182          INFO = -11
00183       END IF
00184       IF( INFO.NE.0 ) THEN
00185          CALL XERBLA( 'DPTSVX', -INFO )
00186          RETURN
00187       END IF
00188 *
00189       IF( NOFACT ) THEN
00190 *
00191 *        Compute the L*D*L' (or U'*D*U) factorization of A.
00192 *
00193          CALL DCOPY( N, D, 1, DF, 1 )
00194          IF( N.GT.1 )
00195      \$      CALL DCOPY( N-1, E, 1, EF, 1 )
00196          CALL DPTTRF( N, DF, EF, INFO )
00197 *
00198 *        Return if INFO is non-zero.
00199 *
00200          IF( INFO.GT.0 )THEN
00201             RCOND = ZERO
00202             RETURN
00203          END IF
00204       END IF
00205 *
00206 *     Compute the norm of the matrix A.
00207 *
00208       ANORM = DLANST( '1', N, D, E )
00209 *
00210 *     Compute the reciprocal of the condition number of A.
00211 *
00212       CALL DPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO )
00213 *
00214 *     Compute the solution vectors X.
00215 *
00216       CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00217       CALL DPTTRS( N, NRHS, DF, EF, X, LDX, INFO )
00218 *
00219 *     Use iterative refinement to improve the computed solutions and
00220 *     compute error bounds and backward error estimates for them.
00221 *
00222       CALL DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR,
00223      \$             WORK, INFO )
00224 *
00225 *     Set INFO = N+1 if the matrix is singular to working precision.
00226 *
00227       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
00228      \$   INFO = N + 1
00229 *
00230       RETURN
00231 *
00232 *     End of DPTSVX
00233 *
00234       END
```