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