SRC/zptsvx.f

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