SRC/zgtsvx.f

   1 *> \brief <b> ZGTSVX computes the solution to system of linear equations A * X = B for GT 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 ZGTSVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgtsvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgtsvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgtsvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
  22 *                          DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
  23 *                          WORK, RWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          FACT, TRANS
  27 *       INTEGER            INFO, LDB, LDX, N, NRHS
  28 *       DOUBLE PRECISION   RCOND
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IPIV( * )
  32 *       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
  33 *       COMPLEX*16         B( LDB, * ), D( * ), DF( * ), DL( * ),
  34 *      $                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
  35 *      $                   WORK( * ), X( LDX, * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> ZGTSVX uses the LU factorization to compute the solution to a complex
  45 *> system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
  46 *> where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
  47 *> matrices.
  48 *>
  49 *> Error bounds on the solution and a condition estimate are also
  50 *> provided.
  51 *> \endverbatim
  52 *
  53 *> \par Description:
  54 *  =================
  55 *>
  56 *> \verbatim
  57 *>
  58 *> The following steps are performed:
  59 *>
  60 *> 1. If FACT = 'N', the LU decomposition is used to factor the matrix A
  61 *>    as A = L * U, where L is a product of permutation and unit lower
  62 *>    bidiagonal matrices and U is upper triangular with nonzeros in
  63 *>    only the main diagonal and first two superdiagonals.
  64 *>
  65 *> 2. If some U(i,i)=0, so that U is exactly singular, then the routine
  66 *>    returns with INFO = i. Otherwise, the factored form of A is used
  67 *>    to estimate the condition number of the matrix A.  If the
  68 *>    reciprocal of the condition number is less than machine precision,
  69 *>    INFO = N+1 is returned as a warning, but the routine still goes on
  70 *>    to solve for X and compute error bounds as described below.
  71 *>
  72 *> 3. The system of equations is solved for X using the factored form
  73 *>    of A.
  74 *>
  75 *> 4. Iterative refinement is applied to improve the computed solution
  76 *>    matrix and calculate error bounds and backward error estimates
  77 *>    for it.
  78 *> \endverbatim
  79 *
  80 *  Arguments:
  81 *  ==========
  82 *
  83 *> \param[in] FACT
  84 *> \verbatim
  85 *>          FACT is CHARACTER*1
  86 *>          Specifies whether or not the factored form of A has been
  87 *>          supplied on entry.
  88 *>          = 'F':  DLF, DF, DUF, DU2, and IPIV contain the factored form
  89 *>                  of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not
  90 *>                  be modified.
  91 *>          = 'N':  The matrix will be copied to DLF, DF, and DUF
  92 *>                  and factored.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] TRANS
  96 *> \verbatim
  97 *>          TRANS is CHARACTER*1
  98 *>          Specifies the form of the system of equations:
  99 *>          = 'N':  A * X = B     (No transpose)
 100 *>          = 'T':  A**T * X = B  (Transpose)
 101 *>          = 'C':  A**H * X = B  (Conjugate transpose)
 102 *> \endverbatim
 103 *>
 104 *> \param[in] N
 105 *> \verbatim
 106 *>          N is INTEGER
 107 *>          The order of the matrix A.  N >= 0.
 108 *> \endverbatim
 109 *>
 110 *> \param[in] NRHS
 111 *> \verbatim
 112 *>          NRHS is INTEGER
 113 *>          The number of right hand sides, i.e., the number of columns
 114 *>          of the matrix B.  NRHS >= 0.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] DL
 118 *> \verbatim
 119 *>          DL is COMPLEX*16 array, dimension (N-1)
 120 *>          The (n-1) subdiagonal elements of A.
 121 *> \endverbatim
 122 *>
 123 *> \param[in] D
 124 *> \verbatim
 125 *>          D is COMPLEX*16 array, dimension (N)
 126 *>          The n diagonal elements of A.
 127 *> \endverbatim
 128 *>
 129 *> \param[in] DU
 130 *> \verbatim
 131 *>          DU is COMPLEX*16 array, dimension (N-1)
 132 *>          The (n-1) superdiagonal elements of A.
 133 *> \endverbatim
 134 *>
 135 *> \param[in,out] DLF
 136 *> \verbatim
 137 *>          DLF is COMPLEX*16 array, dimension (N-1)
 138 *>          If FACT = 'F', then DLF is an input argument and on entry
 139 *>          contains the (n-1) multipliers that define the matrix L from
 140 *>          the LU factorization of A as computed by ZGTTRF.
 141 *>
 142 *>          If FACT = 'N', then DLF is an output argument and on exit
 143 *>          contains the (n-1) multipliers that define the matrix L from
 144 *>          the LU factorization of A.
 145 *> \endverbatim
 146 *>
 147 *> \param[in,out] DF
 148 *> \verbatim
 149 *>          DF is COMPLEX*16 array, dimension (N)
 150 *>          If FACT = 'F', then DF is an input argument and on entry
 151 *>          contains the n diagonal elements of the upper triangular
 152 *>          matrix U from the LU factorization of A.
 153 *>
 154 *>          If FACT = 'N', then DF is an output argument and on exit
 155 *>          contains the n diagonal elements of the upper triangular
 156 *>          matrix U from the LU factorization of A.
 157 *> \endverbatim
 158 *>
 159 *> \param[in,out] DUF
 160 *> \verbatim
 161 *>          DUF is COMPLEX*16 array, dimension (N-1)
 162 *>          If FACT = 'F', then DUF is an input argument and on entry
 163 *>          contains the (n-1) elements of the first superdiagonal of U.
 164 *>
 165 *>          If FACT = 'N', then DUF is an output argument and on exit
 166 *>          contains the (n-1) elements of the first superdiagonal of U.
 167 *> \endverbatim
 168 *>
 169 *> \param[in,out] DU2
 170 *> \verbatim
 171 *>          DU2 is COMPLEX*16 array, dimension (N-2)
 172 *>          If FACT = 'F', then DU2 is an input argument and on entry
 173 *>          contains the (n-2) elements of the second superdiagonal of
 174 *>          U.
 175 *>
 176 *>          If FACT = 'N', then DU2 is an output argument and on exit
 177 *>          contains the (n-2) elements of the second superdiagonal of
 178 *>          U.
 179 *> \endverbatim
 180 *>
 181 *> \param[in,out] IPIV
 182 *> \verbatim
 183 *>          IPIV is INTEGER array, dimension (N)
 184 *>          If FACT = 'F', then IPIV is an input argument and on entry
 185 *>          contains the pivot indices from the LU factorization of A as
 186 *>          computed by ZGTTRF.
 187 *>
 188 *>          If FACT = 'N', then IPIV is an output argument and on exit
 189 *>          contains the pivot indices from the LU factorization of A;
 190 *>          row i of the matrix was interchanged with row IPIV(i).
 191 *>          IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
 192 *>          a row interchange was not required.
 193 *> \endverbatim
 194 *>
 195 *> \param[in] B
 196 *> \verbatim
 197 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 198 *>          The N-by-NRHS right hand side matrix B.
 199 *> \endverbatim
 200 *>
 201 *> \param[in] LDB
 202 *> \verbatim
 203 *>          LDB is INTEGER
 204 *>          The leading dimension of the array B.  LDB >= max(1,N).
 205 *> \endverbatim
 206 *>
 207 *> \param[out] X
 208 *> \verbatim
 209 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
 210 *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
 211 *> \endverbatim
 212 *>
 213 *> \param[in] LDX
 214 *> \verbatim
 215 *>          LDX is INTEGER
 216 *>          The leading dimension of the array X.  LDX >= max(1,N).
 217 *> \endverbatim
 218 *>
 219 *> \param[out] RCOND
 220 *> \verbatim
 221 *>          RCOND is DOUBLE PRECISION
 222 *>          The estimate of the reciprocal condition number of the matrix
 223 *>          A.  If RCOND is less than the machine precision (in
 224 *>          particular, if RCOND = 0), the matrix is singular to working
 225 *>          precision.  This condition is indicated by a return code of
 226 *>          INFO > 0.
 227 *> \endverbatim
 228 *>
 229 *> \param[out] FERR
 230 *> \verbatim
 231 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 232 *>          The estimated forward error bound for each solution vector
 233 *>          X(j) (the j-th column of the solution matrix X).
 234 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 235 *>          is an estimated upper bound for the magnitude of the largest
 236 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 237 *>          largest element in X(j).  The estimate is as reliable as
 238 *>          the estimate for RCOND, and is almost always a slight
 239 *>          overestimate of the true error.
 240 *> \endverbatim
 241 *>
 242 *> \param[out] BERR
 243 *> \verbatim
 244 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 245 *>          The componentwise relative backward error of each solution
 246 *>          vector X(j) (i.e., the smallest relative change in
 247 *>          any element of A or B that makes X(j) an exact solution).
 248 *> \endverbatim
 249 *>
 250 *> \param[out] WORK
 251 *> \verbatim
 252 *>          WORK is COMPLEX*16 array, dimension (2*N)
 253 *> \endverbatim
 254 *>
 255 *> \param[out] RWORK
 256 *> \verbatim
 257 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 258 *> \endverbatim
 259 *>
 260 *> \param[out] INFO
 261 *> \verbatim
 262 *>          INFO is INTEGER
 263 *>          = 0:  successful exit
 264 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 265 *>          > 0:  if INFO = i, and i is
 266 *>                <= N:  U(i,i) is exactly zero.  The factorization
 267 *>                       has not been completed unless i = N, but the
 268 *>                       factor U is exactly singular, so the solution
 269 *>                       and error bounds could not be computed.
 270 *>                       RCOND = 0 is returned.
 271 *>                = N+1: U is nonsingular, but RCOND is less than machine
 272 *>                       precision, meaning that the matrix is singular
 273 *>                       to working precision.  Nevertheless, the
 274 *>                       solution and error bounds are computed because
 275 *>                       there are a number of situations where the
 276 *>                       computed solution can be more accurate than the
 277 *>                       value of RCOND would suggest.
 278 *> \endverbatim
 279 *
 280 *  Authors:
 281 *  ========
 282 *
 283 *> \author Univ. of Tennessee
 284 *> \author Univ. of California Berkeley
 285 *> \author Univ. of Colorado Denver
 286 *> \author NAG Ltd.
 287 *
 288 *> \date September 2012
 289 *
 290 *> \ingroup complex16GTsolve
 291 *
 292 *  =====================================================================
 293       SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
 294      $                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
 295      $                   WORK, RWORK, INFO )
 296 *
 297 *  -- LAPACK driver routine (version 3.4.2) --
 298 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 299 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 300 *     September 2012
 301 *
 302 *     .. Scalar Arguments ..
 303       CHARACTER          FACT, TRANS
 304       INTEGER            INFO, LDB, LDX, N, NRHS
 305       DOUBLE PRECISION   RCOND
 306 *     ..
 307 *     .. Array Arguments ..
 308       INTEGER            IPIV( * )
 309       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
 310       COMPLEX*16         B( LDB, * ), D( * ), DF( * ), DL( * ),
 311      $                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
 312      $                   WORK( * ), X( LDX, * )
 313 *     ..
 314 *
 315 *  =====================================================================
 316 *
 317 *     .. Parameters ..
 318       DOUBLE PRECISION   ZERO
 319       PARAMETER          ( ZERO = 0.0D+0 )
 320 *     ..
 321 *     .. Local Scalars ..
 322       LOGICAL            NOFACT, NOTRAN
 323       CHARACTER          NORM
 324       DOUBLE PRECISION   ANORM
 325 *     ..
 326 *     .. External Functions ..
 327       LOGICAL            LSAME
 328       DOUBLE PRECISION   DLAMCH, ZLANGT
 329       EXTERNAL           LSAME, DLAMCH, ZLANGT
 330 *     ..
 331 *     .. External Subroutines ..
 332       EXTERNAL           XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS,
 333      $                   ZLACPY
 334 *     ..
 335 *     .. Intrinsic Functions ..
 336       INTRINSIC          MAX
 337 *     ..
 338 *     .. Executable Statements ..
 339 *
 340       INFO = 0
 341       NOFACT = LSAME( FACT, 'N' )
 342       NOTRAN = LSAME( TRANS, 'N' )
 343       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
 344          INFO = -1
 345       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 346      $         LSAME( TRANS, 'C' ) ) THEN
 347          INFO = -2
 348       ELSE IF( N.LT.0 ) THEN
 349          INFO = -3
 350       ELSE IF( NRHS.LT.0 ) THEN
 351          INFO = -4
 352       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 353          INFO = -14
 354       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 355          INFO = -16
 356       END IF
 357       IF( INFO.NE.0 ) THEN
 358          CALL XERBLA( 'ZGTSVX', -INFO )
 359          RETURN
 360       END IF
 361 *
 362       IF( NOFACT ) THEN
 363 *
 364 *        Compute the LU factorization of A.
 365 *
 366          CALL ZCOPY( N, D, 1, DF, 1 )
 367          IF( N.GT.1 ) THEN
 368             CALL ZCOPY( N-1, DL, 1, DLF, 1 )
 369             CALL ZCOPY( N-1, DU, 1, DUF, 1 )
 370          END IF
 371          CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
 372 *
 373 *        Return if INFO is non-zero.
 374 *
 375          IF( INFO.GT.0 )THEN
 376             RCOND = ZERO
 377             RETURN
 378          END IF
 379       END IF
 380 *
 381 *     Compute the norm of the matrix A.
 382 *
 383       IF( NOTRAN ) THEN
 384          NORM = '1'
 385       ELSE
 386          NORM = 'I'
 387       END IF
 388       ANORM = ZLANGT( NORM, N, DL, D, DU )
 389 *
 390 *     Compute the reciprocal of the condition number of A.
 391 *
 392       CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
 393      $             INFO )
 394 *
 395 *     Compute the solution vectors X.
 396 *
 397       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
 398       CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
 399      $             INFO )
 400 *
 401 *     Use iterative refinement to improve the computed solutions and
 402 *     compute error bounds and backward error estimates for them.
 403 *
 404       CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
 405      $             B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
 406 *
 407 *     Set INFO = N+1 if the matrix is singular to working precision.
 408 *
 409       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
 410      $   INFO = N + 1
 411 *
 412       RETURN
 413 *
 414 *     End of ZGTSVX
 415 *
 416       END