LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sptsvx.f
Go to the documentation of this file.
00001 *> \brief \b SPTSVX
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SPTSVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sptsvx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sptsvx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sptsvx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
00022 *                          RCOND, FERR, BERR, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          FACT
00026 *       INTEGER            INFO, LDB, LDX, N, NRHS
00027 *       REAL               RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       REAL               B( LDB, * ), BERR( * ), D( * ), DF( * ),
00031 *      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
00032 *      $                   X( LDX, * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> SPTSVX uses the factorization A = L*D*L**T to compute the solution
00042 *> to a real system of linear equations A*X = B, where A is an N-by-N
00043 *> symmetric positive definite tridiagonal matrix and X and B are
00044 *> N-by-NRHS matrices.
00045 *>
00046 *> Error bounds on the solution and a condition estimate are also
00047 *> provided.
00048 *> \endverbatim
00049 *
00050 *> \par Description:
00051 *  =================
00052 *>
00053 *> \verbatim
00054 *>
00055 *> The following steps are performed:
00056 *>
00057 *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L
00058 *>    is a unit lower bidiagonal matrix and D is diagonal.  The
00059 *>    factorization can also be regarded as having the form
00060 *>    A = U**T*D*U.
00061 *>
00062 *> 2. If the leading i-by-i principal minor is not positive definite,
00063 *>    then the routine returns with INFO = i. Otherwise, the factored
00064 *>    form of A is used to estimate the condition number of the matrix
00065 *>    A.  If the reciprocal of the condition number is less than machine
00066 *>    precision, INFO = N+1 is returned as a warning, but the routine
00067 *>    still goes on to solve for X and compute error bounds as
00068 *>    described below.
00069 *>
00070 *> 3. The system of equations is solved for X using the factored form
00071 *>    of A.
00072 *>
00073 *> 4. Iterative refinement is applied to improve the computed solution
00074 *>    matrix and calculate error bounds and backward error estimates
00075 *>    for it.
00076 *> \endverbatim
00077 *
00078 *  Arguments:
00079 *  ==========
00080 *
00081 *> \param[in] FACT
00082 *> \verbatim
00083 *>          FACT is CHARACTER*1
00084 *>          Specifies whether or not the factored form of A has been
00085 *>          supplied on entry.
00086 *>          = 'F':  On entry, DF and EF contain the factored form of A.
00087 *>                  D, E, DF, and EF will not be modified.
00088 *>          = 'N':  The matrix A will be copied to DF and EF and
00089 *>                  factored.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] N
00093 *> \verbatim
00094 *>          N is INTEGER
00095 *>          The order of the matrix A.  N >= 0.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] NRHS
00099 *> \verbatim
00100 *>          NRHS is INTEGER
00101 *>          The number of right hand sides, i.e., the number of columns
00102 *>          of the matrices B and X.  NRHS >= 0.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] D
00106 *> \verbatim
00107 *>          D is REAL array, dimension (N)
00108 *>          The n diagonal elements of the tridiagonal matrix A.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] E
00112 *> \verbatim
00113 *>          E is REAL array, dimension (N-1)
00114 *>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
00115 *> \endverbatim
00116 *>
00117 *> \param[in,out] DF
00118 *> \verbatim
00119 *>          DF is REAL array, dimension (N)
00120 *>          If FACT = 'F', then DF is an input argument and on entry
00121 *>          contains the n diagonal elements of the diagonal matrix D
00122 *>          from the L*D*L**T factorization of A.
00123 *>          If FACT = 'N', then DF is an output argument and on exit
00124 *>          contains the n diagonal elements of the diagonal matrix D
00125 *>          from the L*D*L**T factorization of A.
00126 *> \endverbatim
00127 *>
00128 *> \param[in,out] EF
00129 *> \verbatim
00130 *>          EF is REAL array, dimension (N-1)
00131 *>          If FACT = 'F', then EF is an input argument and on entry
00132 *>          contains the (n-1) subdiagonal elements of the unit
00133 *>          bidiagonal factor L from the L*D*L**T factorization of A.
00134 *>          If FACT = 'N', then EF is an output argument and on exit
00135 *>          contains the (n-1) subdiagonal elements of the unit
00136 *>          bidiagonal factor L from the L*D*L**T factorization of A.
00137 *> \endverbatim
00138 *>
00139 *> \param[in] B
00140 *> \verbatim
00141 *>          B is REAL array, dimension (LDB,NRHS)
00142 *>          The N-by-NRHS right hand side matrix B.
00143 *> \endverbatim
00144 *>
00145 *> \param[in] LDB
00146 *> \verbatim
00147 *>          LDB is INTEGER
00148 *>          The leading dimension of the array B.  LDB >= max(1,N).
00149 *> \endverbatim
00150 *>
00151 *> \param[out] X
00152 *> \verbatim
00153 *>          X is REAL array, dimension (LDX,NRHS)
00154 *>          If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X.
00155 *> \endverbatim
00156 *>
00157 *> \param[in] LDX
00158 *> \verbatim
00159 *>          LDX is INTEGER
00160 *>          The leading dimension of the array X.  LDX >= max(1,N).
00161 *> \endverbatim
00162 *>
00163 *> \param[out] RCOND
00164 *> \verbatim
00165 *>          RCOND is REAL
00166 *>          The reciprocal condition number of the matrix A.  If RCOND
00167 *>          is less than the machine precision (in particular, if
00168 *>          RCOND = 0), the matrix is singular to working precision.
00169 *>          This condition is indicated by a return code of INFO > 0.
00170 *> \endverbatim
00171 *>
00172 *> \param[out] FERR
00173 *> \verbatim
00174 *>          FERR is REAL array, dimension (NRHS)
00175 *>          The forward error bound for each solution vector
00176 *>          X(j) (the j-th column of the solution matrix X).
00177 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
00178 *>          is an estimated upper bound for the magnitude of the largest
00179 *>          element in (X(j) - XTRUE) divided by the magnitude of the
00180 *>          largest element in X(j).
00181 *> \endverbatim
00182 *>
00183 *> \param[out] BERR
00184 *> \verbatim
00185 *>          BERR is REAL array, dimension (NRHS)
00186 *>          The componentwise relative backward error of each solution
00187 *>          vector X(j) (i.e., the smallest relative change in any
00188 *>          element of A or B that makes X(j) an exact solution).
00189 *> \endverbatim
00190 *>
00191 *> \param[out] WORK
00192 *> \verbatim
00193 *>          WORK is REAL array, dimension (2*N)
00194 *> \endverbatim
00195 *>
00196 *> \param[out] INFO
00197 *> \verbatim
00198 *>          INFO is INTEGER
00199 *>          = 0:  successful exit
00200 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00201 *>          > 0:  if INFO = i, and i is
00202 *>                <= N:  the leading minor of order i of A is
00203 *>                       not positive definite, so the factorization
00204 *>                       could not be completed, and the solution has not
00205 *>                       been computed. RCOND = 0 is returned.
00206 *>                = N+1: U is nonsingular, but RCOND is less than machine
00207 *>                       precision, meaning that the matrix is singular
00208 *>                       to working precision.  Nevertheless, the
00209 *>                       solution and error bounds are computed because
00210 *>                       there are a number of situations where the
00211 *>                       computed solution can be more accurate than the
00212 *>                       value of RCOND would suggest.
00213 *> \endverbatim
00214 *
00215 *  Authors:
00216 *  ========
00217 *
00218 *> \author Univ. of Tennessee 
00219 *> \author Univ. of California Berkeley 
00220 *> \author Univ. of Colorado Denver 
00221 *> \author NAG Ltd. 
00222 *
00223 *> \date April 2012
00224 *
00225 *> \ingroup realOTHERcomputational
00226 *
00227 *  =====================================================================
00228       SUBROUTINE SPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
00229      $                   RCOND, FERR, BERR, WORK, INFO )
00230 *
00231 *  -- LAPACK computational routine (version 3.4.1) --
00232 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00233 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00234 *     April 2012
00235 *
00236 *     .. Scalar Arguments ..
00237       CHARACTER          FACT
00238       INTEGER            INFO, LDB, LDX, N, NRHS
00239       REAL               RCOND
00240 *     ..
00241 *     .. Array Arguments ..
00242       REAL               B( LDB, * ), BERR( * ), D( * ), DF( * ),
00243      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
00244      $                   X( LDX, * )
00245 *     ..
00246 *
00247 *  =====================================================================
00248 *
00249 *     .. Parameters ..
00250       REAL               ZERO
00251       PARAMETER          ( ZERO = 0.0E+0 )
00252 *     ..
00253 *     .. Local Scalars ..
00254       LOGICAL            NOFACT
00255       REAL               ANORM
00256 *     ..
00257 *     .. External Functions ..
00258       LOGICAL            LSAME
00259       REAL               SLAMCH, SLANST
00260       EXTERNAL           LSAME, SLAMCH, SLANST
00261 *     ..
00262 *     .. External Subroutines ..
00263       EXTERNAL           SCOPY, SLACPY, SPTCON, SPTRFS, SPTTRF, SPTTRS,
00264      $                   XERBLA
00265 *     ..
00266 *     .. Intrinsic Functions ..
00267       INTRINSIC          MAX
00268 *     ..
00269 *     .. Executable Statements ..
00270 *
00271 *     Test the input parameters.
00272 *
00273       INFO = 0
00274       NOFACT = LSAME( FACT, 'N' )
00275       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
00276          INFO = -1
00277       ELSE IF( N.LT.0 ) THEN
00278          INFO = -2
00279       ELSE IF( NRHS.LT.0 ) THEN
00280          INFO = -3
00281       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00282          INFO = -9
00283       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00284          INFO = -11
00285       END IF
00286       IF( INFO.NE.0 ) THEN
00287          CALL XERBLA( 'SPTSVX', -INFO )
00288          RETURN
00289       END IF
00290 *
00291       IF( NOFACT ) THEN
00292 *
00293 *        Compute the L*D*L**T (or U**T*D*U) factorization of A.
00294 *
00295          CALL SCOPY( N, D, 1, DF, 1 )
00296          IF( N.GT.1 )
00297      $      CALL SCOPY( N-1, E, 1, EF, 1 )
00298          CALL SPTTRF( N, DF, EF, INFO )
00299 *
00300 *        Return if INFO is non-zero.
00301 *
00302          IF( INFO.GT.0 )THEN
00303             RCOND = ZERO
00304             RETURN
00305          END IF
00306       END IF
00307 *
00308 *     Compute the norm of the matrix A.
00309 *
00310       ANORM = SLANST( '1', N, D, E )
00311 *
00312 *     Compute the reciprocal of the condition number of A.
00313 *
00314       CALL SPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO )
00315 *
00316 *     Compute the solution vectors X.
00317 *
00318       CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00319       CALL SPTTRS( N, NRHS, DF, EF, X, LDX, INFO )
00320 *
00321 *     Use iterative refinement to improve the computed solutions and
00322 *     compute error bounds and backward error estimates for them.
00323 *
00324       CALL SPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR,
00325      $             WORK, INFO )
00326 *
00327 *     Set INFO = N+1 if the matrix is singular to working precision.
00328 *
00329       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00330      $   INFO = N + 1
00331 *
00332       RETURN
00333 *
00334 *     End of SPTSVX
00335 *
00336       END
 All Files Functions