SRC/dptsvx.f

   1 *> \brief <b> DPTSVX computes the solution to system of linear equations A * X = B for PT matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DPTSVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptsvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptsvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptsvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
  22 *                          RCOND, FERR, BERR, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          FACT
  26 *       INTEGER            INFO, LDB, LDX, N, NRHS
  27 *       DOUBLE PRECISION   RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
  31 *      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
  32 *      $                   X( LDX, * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> DPTSVX uses the factorization A = L*D*L**T to compute the solution
  42 *> to a real system of linear equations A*X = B, where A is an N-by-N
  43 *> symmetric positive definite tridiagonal matrix and X and B are
  44 *> N-by-NRHS matrices.
  45 *>
  46 *> Error bounds on the solution and a condition estimate are also
  47 *> provided.
  48 *> \endverbatim
  49 *
  50 *> \par Description:
  51 *  =================
  52 *>
  53 *> \verbatim
  54 *>
  55 *> The following steps are performed:
  56 *>
  57 *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L
  58 *>    is a unit lower bidiagonal matrix and D is diagonal.  The
  59 *>    factorization can also be regarded as having the form
  60 *>    A = U**T*D*U.
  61 *>
  62 *> 2. If the leading i-by-i principal minor is not positive definite,
  63 *>    then the routine returns with INFO = i. Otherwise, the factored
  64 *>    form of A is used to estimate the condition number of the matrix
  65 *>    A.  If the reciprocal of the condition number is less than machine
  66 *>    precision, INFO = N+1 is returned as a warning, but the routine
  67 *>    still goes on to solve for X and compute error bounds as
  68 *>    described below.
  69 *>
  70 *> 3. The system of equations is solved for X using the factored form
  71 *>    of A.
  72 *>
  73 *> 4. Iterative refinement is applied to improve the computed solution
  74 *>    matrix and calculate error bounds and backward error estimates
  75 *>    for it.
  76 *> \endverbatim
  77 *
  78 *  Arguments:
  79 *  ==========
  80 *
  81 *> \param[in] FACT
  82 *> \verbatim
  83 *>          FACT is CHARACTER*1
  84 *>          Specifies whether or not the factored form of A has been
  85 *>          supplied on entry.
  86 *>          = 'F':  On entry, DF and EF contain the factored form of A.
  87 *>                  D, E, DF, and EF will not be modified.
  88 *>          = 'N':  The matrix A will be copied to DF and EF and
  89 *>                  factored.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] N
  93 *> \verbatim
  94 *>          N is INTEGER
  95 *>          The order of the matrix A.  N >= 0.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] NRHS
  99 *> \verbatim
 100 *>          NRHS is INTEGER
 101 *>          The number of right hand sides, i.e., the number of columns
 102 *>          of the matrices B and X.  NRHS >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] D
 106 *> \verbatim
 107 *>          D is DOUBLE PRECISION array, dimension (N)
 108 *>          The n diagonal elements of the tridiagonal matrix A.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] E
 112 *> \verbatim
 113 *>          E is DOUBLE PRECISION array, dimension (N-1)
 114 *>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
 115 *> \endverbatim
 116 *>
 117 *> \param[in,out] DF
 118 *> \verbatim
 119 *>          DF is DOUBLE PRECISION array, dimension (N)
 120 *>          If FACT = 'F', then DF is an input argument and on entry
 121 *>          contains the n diagonal elements of the diagonal matrix D
 122 *>          from the L*D*L**T factorization of A.
 123 *>          If FACT = 'N', then DF is an output argument and on exit
 124 *>          contains the n diagonal elements of the diagonal matrix D
 125 *>          from the L*D*L**T factorization of A.
 126 *> \endverbatim
 127 *>
 128 *> \param[in,out] EF
 129 *> \verbatim
 130 *>          EF is DOUBLE PRECISION array, dimension (N-1)
 131 *>          If FACT = 'F', then EF is an input argument and on entry
 132 *>          contains the (n-1) subdiagonal elements of the unit
 133 *>          bidiagonal factor L from the L*D*L**T factorization of A.
 134 *>          If FACT = 'N', then EF is an output argument and on exit
 135 *>          contains the (n-1) subdiagonal elements of the unit
 136 *>          bidiagonal factor L from the L*D*L**T factorization of A.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] B
 140 *> \verbatim
 141 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
 142 *>          The N-by-NRHS right hand side matrix B.
 143 *> \endverbatim
 144 *>
 145 *> \param[in] LDB
 146 *> \verbatim
 147 *>          LDB is INTEGER
 148 *>          The leading dimension of the array B.  LDB >= max(1,N).
 149 *> \endverbatim
 150 *>
 151 *> \param[out] X
 152 *> \verbatim
 153 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
 154 *>          If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X.
 155 *> \endverbatim
 156 *>
 157 *> \param[in] LDX
 158 *> \verbatim
 159 *>          LDX is INTEGER
 160 *>          The leading dimension of the array X.  LDX >= max(1,N).
 161 *> \endverbatim
 162 *>
 163 *> \param[out] RCOND
 164 *> \verbatim
 165 *>          RCOND is DOUBLE PRECISION
 166 *>          The reciprocal condition number of the matrix A.  If RCOND
 167 *>          is less than the machine precision (in particular, if
 168 *>          RCOND = 0), the matrix is singular to working precision.
 169 *>          This condition is indicated by a return code of INFO > 0.
 170 *> \endverbatim
 171 *>
 172 *> \param[out] FERR
 173 *> \verbatim
 174 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 175 *>          The forward error bound for each solution vector
 176 *>          X(j) (the j-th column of the solution matrix X).
 177 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 178 *>          is an estimated upper bound for the magnitude of the largest
 179 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 180 *>          largest element in X(j).
 181 *> \endverbatim
 182 *>
 183 *> \param[out] BERR
 184 *> \verbatim
 185 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 186 *>          The componentwise relative backward error of each solution
 187 *>          vector X(j) (i.e., the smallest relative change in any
 188 *>          element of A or B that makes X(j) an exact solution).
 189 *> \endverbatim
 190 *>
 191 *> \param[out] WORK
 192 *> \verbatim
 193 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
 194 *> \endverbatim
 195 *>
 196 *> \param[out] INFO
 197 *> \verbatim
 198 *>          INFO is INTEGER
 199 *>          = 0:  successful exit
 200 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 201 *>          > 0:  if INFO = i, and i is
 202 *>                <= N:  the leading minor of order i of A is
 203 *>                       not positive definite, so the factorization
 204 *>                       could not be completed, and the solution has not
 205 *>                       been computed. RCOND = 0 is returned.
 206 *>                = N+1: U is nonsingular, but RCOND is less than machine
 207 *>                       precision, meaning that the matrix is singular
 208 *>                       to working precision.  Nevertheless, the
 209 *>                       solution and error bounds are computed because
 210 *>                       there are a number of situations where the
 211 *>                       computed solution can be more accurate than the
 212 *>                       value of RCOND would suggest.
 213 *> \endverbatim
 214 *
 215 *  Authors:
 216 *  ========
 217 *
 218 *> \author Univ. of Tennessee
 219 *> \author Univ. of California Berkeley
 220 *> \author Univ. of Colorado Denver
 221 *> \author NAG Ltd.
 222 *
 223 *> \date September 2012
 224 *
 225 *> \ingroup doublePTsolve
 226 *
 227 *  =====================================================================
 228       SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
 229      $                   RCOND, FERR, BERR, WORK, INFO )
 230 *
 231 *  -- LAPACK driver routine (version 3.4.2) --
 232 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 233 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 234 *     September 2012
 235 *
 236 *     .. Scalar Arguments ..
 237       CHARACTER          FACT
 238       INTEGER            INFO, LDB, LDX, N, NRHS
 239       DOUBLE PRECISION   RCOND
 240 *     ..
 241 *     .. Array Arguments ..
 242       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
 243      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
 244      $                   X( LDX, * )
 245 *     ..
 246 *
 247 *  =====================================================================
 248 *
 249 *     .. Parameters ..
 250       DOUBLE PRECISION   ZERO
 251       PARAMETER          ( ZERO = 0.0D+0 )
 252 *     ..
 253 *     .. Local Scalars ..
 254       LOGICAL            NOFACT
 255       DOUBLE PRECISION   ANORM
 256 *     ..
 257 *     .. External Functions ..
 258       LOGICAL            LSAME
 259       DOUBLE PRECISION   DLAMCH, DLANST
 260       EXTERNAL           LSAME, DLAMCH, DLANST
 261 *     ..
 262 *     .. External Subroutines ..
 263       EXTERNAL           DCOPY, DLACPY, DPTCON, DPTRFS, DPTTRF, DPTTRS,
 264      $                   XERBLA
 265 *     ..
 266 *     .. Intrinsic Functions ..
 267       INTRINSIC          MAX
 268 *     ..
 269 *     .. Executable Statements ..
 270 *
 271 *     Test the input parameters.
 272 *
 273       INFO = 0
 274       NOFACT = LSAME( FACT, 'N' )
 275       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
 276          INFO = -1
 277       ELSE IF( N.LT.0 ) THEN
 278          INFO = -2
 279       ELSE IF( NRHS.LT.0 ) THEN
 280          INFO = -3
 281       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 282          INFO = -9
 283       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 284          INFO = -11
 285       END IF
 286       IF( INFO.NE.0 ) THEN
 287          CALL XERBLA( 'DPTSVX', -INFO )
 288          RETURN
 289       END IF
 290 *
 291       IF( NOFACT ) THEN
 292 *
 293 *        Compute the L*D*L**T (or U**T*D*U) factorization of A.
 294 *
 295          CALL DCOPY( N, D, 1, DF, 1 )
 296          IF( N.GT.1 )
 297      $      CALL DCOPY( N-1, E, 1, EF, 1 )
 298          CALL DPTTRF( N, DF, EF, INFO )
 299 *
 300 *        Return if INFO is non-zero.
 301 *
 302          IF( INFO.GT.0 )THEN
 303             RCOND = ZERO
 304             RETURN
 305          END IF
 306       END IF
 307 *
 308 *     Compute the norm of the matrix A.
 309 *
 310       ANORM = DLANST( '1', N, D, E )
 311 *
 312 *     Compute the reciprocal of the condition number of A.
 313 *
 314       CALL DPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO )
 315 *
 316 *     Compute the solution vectors X.
 317 *
 318       CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
 319       CALL DPTTRS( N, NRHS, DF, EF, X, LDX, INFO )
 320 *
 321 *     Use iterative refinement to improve the computed solutions and
 322 *     compute error bounds and backward error estimates for them.
 323 *
 324       CALL DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR,
 325      $             WORK, INFO )
 326 *
 327 *     Set INFO = N+1 if the matrix is singular to working precision.
 328 *
 329       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
 330      $   INFO = N + 1
 331 *
 332       RETURN
 333 *
 334 *     End of DPTSVX
 335 *
 336       END