 LAPACK 3.3.1 Linear Algebra PACKage

# cgtsvx.f

Go to the documentation of this file.
```00001       SUBROUTINE CGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
00002      \$                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
00003      \$                   WORK, RWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          FACT, TRANS
00012       INTEGER            INFO, LDB, LDX, N, NRHS
00013       REAL               RCOND
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IPIV( * )
00017       REAL               BERR( * ), FERR( * ), RWORK( * )
00018       COMPLEX            B( LDB, * ), D( * ), DF( * ), DL( * ),
00019      \$                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
00020      \$                   WORK( * ), X( LDX, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CGTSVX uses the LU factorization to compute the solution to a complex
00027 *  system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
00028 *  where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
00029 *  matrices.
00030 *
00031 *  Error bounds on the solution and a condition estimate are also
00032 *  provided.
00033 *
00034 *  Description
00035 *  ===========
00036 *
00037 *  The following steps are performed:
00038 *
00039 *  1. If FACT = 'N', the LU decomposition is used to factor the matrix A
00040 *     as A = L * U, where L is a product of permutation and unit lower
00041 *     bidiagonal matrices and U is upper triangular with nonzeros in
00042 *     only the main diagonal and first two superdiagonals.
00043 *
00044 *  2. If some U(i,i)=0, so that U is exactly singular, then the routine
00045 *     returns with INFO = i. Otherwise, the factored form of A is used
00046 *     to estimate the condition number of the matrix A.  If the
00047 *     reciprocal of the condition number is less than machine precision,
00048 *     INFO = N+1 is returned as a warning, but the routine still goes on
00049 *     to solve for X and compute error bounds as described below.
00050 *
00051 *  3. The system of equations is solved for X using the factored form
00052 *     of A.
00053 *
00054 *  4. Iterative refinement is applied to improve the computed solution
00055 *     matrix and calculate error bounds and backward error estimates
00056 *     for it.
00057 *
00058 *  Arguments
00059 *  =========
00060 *
00061 *  FACT    (input) CHARACTER*1
00062 *          Specifies whether or not the factored form of A has been
00063 *          supplied on entry.
00064 *          = 'F':  DLF, DF, DUF, DU2, and IPIV contain the factored form
00065 *                  of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not
00066 *                  be modified.
00067 *          = 'N':  The matrix will be copied to DLF, DF, and DUF
00068 *                  and factored.
00069 *
00070 *  TRANS   (input) CHARACTER*1
00071 *          Specifies the form of the system of equations:
00072 *          = 'N':  A * X = B     (No transpose)
00073 *          = 'T':  A**T * X = B  (Transpose)
00074 *          = 'C':  A**H * X = B  (Conjugate transpose)
00075 *
00076 *  N       (input) INTEGER
00077 *          The order of the matrix A.  N >= 0.
00078 *
00079 *  NRHS    (input) INTEGER
00080 *          The number of right hand sides, i.e., the number of columns
00081 *          of the matrix B.  NRHS >= 0.
00082 *
00083 *  DL      (input) COMPLEX array, dimension (N-1)
00084 *          The (n-1) subdiagonal elements of A.
00085 *
00086 *  D       (input) COMPLEX array, dimension (N)
00087 *          The n diagonal elements of A.
00088 *
00089 *  DU      (input) COMPLEX array, dimension (N-1)
00090 *          The (n-1) superdiagonal elements of A.
00091 *
00092 *  DLF     (input or output) COMPLEX array, dimension (N-1)
00093 *          If FACT = 'F', then DLF is an input argument and on entry
00094 *          contains the (n-1) multipliers that define the matrix L from
00095 *          the LU factorization of A as computed by CGTTRF.
00096 *
00097 *          If FACT = 'N', then DLF is an output argument and on exit
00098 *          contains the (n-1) multipliers that define the matrix L from
00099 *          the LU factorization of A.
00100 *
00101 *  DF      (input or output) COMPLEX array, dimension (N)
00102 *          If FACT = 'F', then DF is an input argument and on entry
00103 *          contains the n diagonal elements of the upper triangular
00104 *          matrix U from the LU factorization of A.
00105 *
00106 *          If FACT = 'N', then DF is an output argument and on exit
00107 *          contains the n diagonal elements of the upper triangular
00108 *          matrix U from the LU factorization of A.
00109 *
00110 *  DUF     (input or output) COMPLEX array, dimension (N-1)
00111 *          If FACT = 'F', then DUF is an input argument and on entry
00112 *          contains the (n-1) elements of the first superdiagonal of U.
00113 *
00114 *          If FACT = 'N', then DUF is an output argument and on exit
00115 *          contains the (n-1) elements of the first superdiagonal of U.
00116 *
00117 *  DU2     (input or output) COMPLEX array, dimension (N-2)
00118 *          If FACT = 'F', then DU2 is an input argument and on entry
00119 *          contains the (n-2) elements of the second superdiagonal of
00120 *          U.
00121 *
00122 *          If FACT = 'N', then DU2 is an output argument and on exit
00123 *          contains the (n-2) elements of the second superdiagonal of
00124 *          U.
00125 *
00126 *  IPIV    (input or output) INTEGER array, dimension (N)
00127 *          If FACT = 'F', then IPIV is an input argument and on entry
00128 *          contains the pivot indices from the LU factorization of A as
00129 *          computed by CGTTRF.
00130 *
00131 *          If FACT = 'N', then IPIV is an output argument and on exit
00132 *          contains the pivot indices from the LU factorization of A;
00133 *          row i of the matrix was interchanged with row IPIV(i).
00134 *          IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
00135 *          a row interchange was not required.
00136 *
00137 *  B       (input) COMPLEX array, dimension (LDB,NRHS)
00138 *          The N-by-NRHS right hand side matrix B.
00139 *
00140 *  LDB     (input) INTEGER
00141 *          The leading dimension of the array B.  LDB >= max(1,N).
00142 *
00143 *  X       (output) COMPLEX array, dimension (LDX,NRHS)
00144 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
00145 *
00146 *  LDX     (input) INTEGER
00147 *          The leading dimension of the array X.  LDX >= max(1,N).
00148 *
00149 *  RCOND   (output) REAL
00150 *          The estimate of the reciprocal condition number of the matrix
00151 *          A.  If RCOND is less than the machine precision (in
00152 *          particular, if RCOND = 0), the matrix is singular to working
00153 *          precision.  This condition is indicated by a return code of
00154 *          INFO > 0.
00155 *
00156 *  FERR    (output) REAL array, dimension (NRHS)
00157 *          The estimated forward error bound for each solution vector
00158 *          X(j) (the j-th column of the solution matrix X).
00159 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00160 *          is an estimated upper bound for the magnitude of the largest
00161 *          element in (X(j) - XTRUE) divided by the magnitude of the
00162 *          largest element in X(j).  The estimate is as reliable as
00163 *          the estimate for RCOND, and is almost always a slight
00164 *          overestimate of the true error.
00165 *
00166 *  BERR    (output) REAL array, dimension (NRHS)
00167 *          The componentwise relative backward error of each solution
00168 *          vector X(j) (i.e., the smallest relative change in
00169 *          any element of A or B that makes X(j) an exact solution).
00170 *
00171 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00172 *
00173 *  RWORK   (workspace) REAL array, dimension (N)
00174 *
00175 *  INFO    (output) INTEGER
00176 *          = 0:  successful exit
00177 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00178 *          > 0:  if INFO = i, and i is
00179 *                <= N:  U(i,i) is exactly zero.  The factorization
00180 *                       has not been completed unless i = N, but the
00181 *                       factor U is exactly singular, so the solution
00182 *                       and error bounds could not be computed.
00183 *                       RCOND = 0 is returned.
00184 *                = N+1: U is nonsingular, but RCOND is less than machine
00185 *                       precision, meaning that the matrix is singular
00186 *                       to working precision.  Nevertheless, the
00187 *                       solution and error bounds are computed because
00188 *                       there are a number of situations where the
00189 *                       computed solution can be more accurate than the
00190 *                       value of RCOND would suggest.
00191 *
00192 *  =====================================================================
00193 *
00194 *     .. Parameters ..
00195       REAL               ZERO
00196       PARAMETER          ( ZERO = 0.0E+0 )
00197 *     ..
00198 *     .. Local Scalars ..
00199       LOGICAL            NOFACT, NOTRAN
00200       CHARACTER          NORM
00201       REAL               ANORM
00202 *     ..
00203 *     .. External Functions ..
00204       LOGICAL            LSAME
00205       REAL               CLANGT, SLAMCH
00206       EXTERNAL           LSAME, CLANGT, SLAMCH
00207 *     ..
00208 *     .. External Subroutines ..
00209       EXTERNAL           CCOPY, CGTCON, CGTRFS, CGTTRF, CGTTRS, CLACPY,
00210      \$                   XERBLA
00211 *     ..
00212 *     .. Intrinsic Functions ..
00213       INTRINSIC          MAX
00214 *     ..
00215 *     .. Executable Statements ..
00216 *
00217       INFO = 0
00218       NOFACT = LSAME( FACT, 'N' )
00219       NOTRAN = LSAME( TRANS, 'N' )
00220       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
00221          INFO = -1
00222       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00223      \$         LSAME( TRANS, 'C' ) ) THEN
00224          INFO = -2
00225       ELSE IF( N.LT.0 ) THEN
00226          INFO = -3
00227       ELSE IF( NRHS.LT.0 ) THEN
00228          INFO = -4
00229       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00230          INFO = -14
00231       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00232          INFO = -16
00233       END IF
00234       IF( INFO.NE.0 ) THEN
00235          CALL XERBLA( 'CGTSVX', -INFO )
00236          RETURN
00237       END IF
00238 *
00239       IF( NOFACT ) THEN
00240 *
00241 *        Compute the LU factorization of A.
00242 *
00243          CALL CCOPY( N, D, 1, DF, 1 )
00244          IF( N.GT.1 ) THEN
00245             CALL CCOPY( N-1, DL, 1, DLF, 1 )
00246             CALL CCOPY( N-1, DU, 1, DUF, 1 )
00247          END IF
00248          CALL CGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
00249 *
00250 *        Return if INFO is non-zero.
00251 *
00252          IF( INFO.GT.0 )THEN
00253             RCOND = ZERO
00254             RETURN
00255          END IF
00256       END IF
00257 *
00258 *     Compute the norm of the matrix A.
00259 *
00260       IF( NOTRAN ) THEN
00261          NORM = '1'
00262       ELSE
00263          NORM = 'I'
00264       END IF
00265       ANORM = CLANGT( NORM, N, DL, D, DU )
00266 *
00267 *     Compute the reciprocal of the condition number of A.
00268 *
00269       CALL CGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
00270      \$             INFO )
00271 *
00272 *     Compute the solution vectors X.
00273 *
00274       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00275       CALL CGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
00276      \$             INFO )
00277 *
00278 *     Use iterative refinement to improve the computed solutions and
00279 *     compute error bounds and backward error estimates for them.
00280 *
00281       CALL CGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
00282      \$             B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
00283 *
00284 *     Set INFO = N+1 if the matrix is singular to working precision.
00285 *
00286       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00287      \$   INFO = N + 1
00288 *
00289       RETURN
00290 *
00291 *     End of CGTSVX
00292 *
00293       END
```